AVX 指令集并行技术优化积分计算圆周率 π

// AVX_PI.cpp : 定义控制台应用程序的入口点。
//
#include "stdafx.h"
#include <iostream>
#include <immintrin.h>
#include <time.h>
#include "timer.h"
double compute_pi_naive(size_t dt)
{
  double pi = 0.0;
  double delta = 1.0 / dt;
  for (size_t i = 0; i < dt; i++)
  {
    double x = (double)i / dt;
    pi += delta / (1.0 + x*x);
  }
  return pi*4.0;
}
double compute_pi_avx(size_t dt)
{
  double pi = 0.0;
  double delta = 1.0 / dt;
  __m256d ymm0, ymm1, ymm2, ymm3, ymm4;
  ymm0 = _mm256_set1_pd(1.0); // 注意并不是对应一条指令，而是多条指令混合而成
  ymm1 = _mm256_set1_pd(delta);
  ymm2 = _mm256_set_pd(delta * 3, delta * 2, delta, 0.0); // 0.0， delta, 2*delta, 3*delta
  ymm4 = _mm256_setzero_pd();
  for (int i = 0; i <= dt - 4; i += 4)
  {
    ymm3 = _mm256_set1_pd(i * delta); // 归一化构造 x（积分区间为 0-1）
    ymm3 = _mm256_add_pd(ymm3, ymm2); // 分别累加四项（分别相差一个delta）
    ymm3 = _mm256_mul_pd(ymm3, ymm3); // 平方
    ymm3 = _mm256_add_pd(ymm0, ymm3); // + 1
    ymm3 = _mm256_div_pd(ymm1, ymm3);   // delta 向量除
    ymm4 = _mm256_add_pd(ymm4, ymm3); // 累加
  }
#ifndef _WIN32
  double tmp[4] __attribute__((aligned(32)));
#else
//#pragma pack(8)
  __declspec(align(32)) // arm 下的内存边界对齐使用 __align
  double tmp[4];
//#pragma pack()
#endif
  _mm256_store_pd(tmp, ymm4);
  pi = tmp[0] + tmp[1] + tmp[2] + tmp[3];
  return pi * 4.0;
}
int _tmain(int argc, _TCHAR* argv[])
{
  size_t dt = 1e10;
  Timer timer;
  double pi_naive = compute_pi_naive(dt);
  printf("Pi = %lf, %lfms\n", pi_naive, timer.elapsed());
  timer.restart();
  double pi_avx = compute_pi_avx(dt);
  printf("Pi = %lf, %lfms\n", pi_avx, timer.elapsed());
  return 0;
}