通过 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;
}

最后输出的结果，获得了 3 倍多的加速（耗时不到原来串行版本的 1/4 ）