以”计算圆周率“为例说明AVX相关指令

圆周率计算公式如下：

AVX相关指令介绍：https://software.intel.com/sites/landingpage/IntrinsicsGuide/

AVX指令中有计算反正切函数的指令，但本文的目的是为了对比验证AVX相关指令实现SIMD并行计算对于计算效率的提升。因此，我们选取定积分部分作为计算圆周率的衡量方式，即将计算圆周率视作大量的可迭代的加法运算。

本文相关代码：https://github.com/vistart/IntelIntrinsics

一、准备工作

1. 引用头文件

使用AVX相关协议需要引用头文件”immintrin.h”：

#include "immintrin.h"

2. 验证自己的处理器对AVX指令的支持程度

方法一：如果是英特尔处理器，且知道处理器的型号，可以到英特尔官方网站查询。

以英特尔酷睿i9-9900K为例，查询结果如下：

从图1中可以看出，英特尔酷睿i9-9900K支持AVX2指令集，但不支持AVX-512指令集。

再以英特尔酷睿i5-1135G7为例，查询结果如下：

从图2中可以看出，英特尔酷睿i5-1135G7支持AVX-512指令集。

方法二：如果不知道处理器的型号，或没有可查看处理器指令集支持情况的页面，可以下载CPU-Z后运行CPU-Z查询。CPU-Z会直接显示当前处理器支持的指令集情况。

以 AMD Ryzen 9 3900X 为例，指令集支持情况如下：

从图3中可以看出，AMD Ryzen 9 3900X支持AVX2指令集，不支持AVX-512指令集。

二、原始计算方式

原始计算方式是仅使用默认四则运算来计算圆周率。源代码如下：

double compute_pi_native(size_t dt)
{
	double pi = 0.0;
	double const delta = 1.0 / dt;
	for (size_t i = 0; i < dt; i++)
	{
		double x = i * delta;
		pi += delta / (1 + x * x);
	}
	return pi * 4;
}

三、使用AVX2相关指令计算

double compute_pi_avx256(size_t dt)
{
    double pi = 0.0;
    double const 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);
    ymm4 = _mm256_setzero_pd();
    for (int i = 0; i < dt - 4; i += 4) {
        ymm3 = _mm256_set1_pd(i * delta); // x = i * delta
        ymm3 = _mm256_add_pd(ymm3, ymm2); // x1 = x, x2 = x + delta, x3 = x + 2 * delta, x4 = x + 3 * delta
        ymm3 = _mm256_mul_pd(ymm3, ymm3); // x * x
        ymm3 = _mm256_add_pd(ymm0, ymm3); // 1 + x * x
        ymm3 = _mm256_div_pd(ymm1, ymm3); // delta / (1 + x * x)
        ymm4 = _mm256_add_pd(ymm4, ymm3); // pi += delta / (1 + x * x)
    }
    double tmp[4];
    _mm256_store_pd(tmp, ymm4);
    pi += tmp[0] + tmp[1] + tmp[2] + tmp[3];
    return pi * 4.0;
}

四、使用AVX-512相关指令计算

double compute_pi_avx512(size_t dt)
{
	double pi = 0.0;
	double const delta = 1.0 / dt;
	__m512d ymm0, ymm1, ymm2, ymm3, ymm4;
	ymm0 = _mm512_set1_pd(1.0);
	ymm1 = _mm512_set1_pd(delta);
	ymm2 = _mm512_set_pd(delta * 7, delta * 6, delta * 5, delta * 4, delta * 3, delta * 2, delta, 0.0);
	ymm4 = _mm512_setzero_pd();
	for (int i = 0; i < dt - 8; i += 8)
	{
		ymm3 = _mm512_set1_pd(i * delta);
		ymm3 = _mm512_add_pd(ymm3, ymm2);
		ymm3 = _mm512_mul_pd(ymm3, ymm3);
		ymm3 = _mm512_add_pd(ymm0, ymm3);
		ymm3 = _mm512_div_pd(ymm1, ymm3);
		ymm4 = _mm512_add_pd(ymm4, ymm3);
	}
	double tmp[8];
	_mm512_store_pd(tmp, ymm4);
	pi += tmp[0] + tmp[1] + tmp[2] + tmp[3] + tmp[4] + tmp[5] + tmp[6] + tmp[7];
	return pi * 4;
}

五、实验

以下实验中，考量到每次运行耗时会有细微差异，每次Debug环境实验将运行十次，并计算耗时均值和方差；而 Release 模式会存在累进优化现象，即相同变量多次赋值只记录最后一次，因而 Release 下循环运行多次并分次计时会产生失真，故 Release 模式只通过观察方式记录十次耗时的均值和方差。

1. 实验环境一

编译环境： Visual Studio 2019，C++ 17

处理器：英特尔酷睿i5-1135G7，内存：16GB 4266MHz 四通道

在 Debug, x64 模式下，计算结果如下：

原始：3.14159265405554100070162348856684，30454毫秒

AVX2：3.14159265126151776570395668386482，15716毫秒

AVX-512：3.14159264753617062382318181334995，8362毫秒

在 Release, x64 模式下，计算结果如下：

原始：3.14159265405554100070162348856684，12325毫秒

AVX2：3.14159265126151776570395668386482，1045毫秒

AVX-512：3.14159264753617062382318181334995，1062毫秒

从以上结果可得出：Debug 模式下，AVX2指令的运行耗时比原始缩短了一半，而AVX-512又比AVX2缩短了一半。在 Release 模式下，AVX2指令的运行耗时相对原始缩短到十二分之一，而AVX-512与AVX2的运行耗时相差不大。反复测试若干次后，基本耗时和相互比例相差极小，可以认为此结果相对稳定。

2. 实验环境二

编译环境：Visual Studio 2019，C++ 17

处理器：英特尔酷睿 i9-9900K，内存：96GB，3000MHz 双通道

在 Debug, x64 模式下，计算结果如下：

原始：3.14159265405554100070162348856684，耗时均值18070.70毫秒，方差32.15毫秒

AVX2：3.14159265126151776570395668386482，耗时均值12475.30毫秒，方差35.33毫秒

在 Release, x64 模式下，计算结果如下：

原始：3.14159265405554100070162348856684，8650毫秒

AVX2：3.14159265126151776570395668386482，919毫秒

从以上结果可得出：Debug 模式下，AVX2指令的运行耗时约为原始的三分之二。在 Release 模式下，AVX2指令的运行还是低于原始的九分之一。因该处理器不支持AVX-512指令集，故此测试不涉及AVX-512指令集。

3. 实验环境三

编译环境：Visual Studio 2019，C++ 17

处理器：AMD Ryzen 3900X，内存：64GB，3000MHz 双通道

在 Debug, x64 模式下，计算结果如下：

原始：3.14159265405554100070162348856684，5770毫秒

AVX2：3.14159265126151776570395668386482，51343毫秒

在 Release, x64 模式下，计算结果如下：

原始：3.14159265405554100070162348856684，2616毫秒

AVX2：3.14159265126151776570395668386482，651毫秒

从以上结果可得出：Debug 模式下，AVX2指令的运行耗时接近原始的九倍。在 Release 模式下，AVX2指令的运行还是低于原始的四分之一。因该处理器不支持AVX-512指令集，故此测试不涉及AVX-512指令集。

六、总结

从实验结果中可看出，AVX指令集对并行计算效率提升十分显著。但在 Release 模式下 AVX2 和 AVX-512 未显示出明显差距，甚至 AVX-512 比 AVX2 还慢一点。不论是哪种方式，计算得出的结果有效位按十进制算，在小数点后7位。这个精度足以满足大多数计算需求，但更高精度计算仍待改进算法。接下来将深入分析其中原理，也欢迎广大网友提供建议。