AVX计算精度

Posted 2023-02-16

【中文标题】AVX计算精度

【英文标题】：AVX calculation precision

【发布时间】：2018-08-19 12:24:01

【问题描述】：

我编写了一个程序来显示 mandelbrot 集。为了加快速度，我通过 <immintrin.h> 标头使用了 AVX（实际上是 AVX2）指令。 问题是：AVX 计算（双精度）的结果有伪影，并且与使用“正常”双精度计算时的结果不同。 详细来说，有一个函数getIterationCount 计算迭代次数，直到 mandelbrot 序列超过 4，或者如果在前 N 步中序列不超过 4，则假设该点包含在集合中。 代码如下所示：

#include "stdafx.h"
#include <iostream>
#include <complex>
#include <immintrin.h>

class MandelbrotSet 
public:
    int getIterationCount(const std::complex<double>, const int) const noexcept;
    __m256i getIterationCount(__m256d cReal, __m256d cIm, unsigned maxIterations) const noexcept;
;

inline int MandelbrotSet::getIterationCount(const std::complex<double> c, const int maxIterations) const noexcept

    double currentReal = 0;
    double currentIm = 0;
    double realSquare;
    double imSquare;
    for (int i = 0; i < maxIterations; ++i) 
        realSquare = currentReal * currentReal;
        imSquare = currentIm * currentIm;
        currentIm = 2 * currentReal * currentIm + c.imag();
        currentReal = realSquare - imSquare + c.real();
        if (realSquare + imSquare >= 4) 
            return i;
        
    
    return -1;


const __m256i negone = _mm256_set_epi64x(-1, -1, -1, -1);
const __m256i one = _mm256_set_epi64x(1, 1, 1, 1);
const __m256d two = _mm256_set_pd(2, 2, 2, 2);
const __m256d four = _mm256_set_pd(4, 4, 4, 4);

//calculates for i = 0,1,2,3
//output[i] = if ctrl[i] == 0b11...1 then onTrue[i] else onFalse[i]
inline __m256i _mm256_select_si256(__m256i onTrue, __m256i onFalse, __m256i ctrl) 
    return _mm256_or_si256(_mm256_and_si256(onTrue, ctrl), _mm256_and_si256(onFalse, _mm256_xor_si256(negone, ctrl)));


inline __m256i MandelbrotSet::getIterationCount(__m256d cReal, __m256d cIm, unsigned maxIterations) const noexcept 
    __m256i result = _mm256_set_epi64x(0, 0, 0, 0);
    __m256d currentReal = _mm256_set_pd(0, 0, 0, 0);
    __m256d currentIm = _mm256_set_pd(0, 0, 0, 0);
    __m256d realSquare;
    __m256d imSquare;
    for (unsigned i = 0; i <= maxIterations; ++i)
    
        realSquare = _mm256_mul_pd(currentReal, currentReal);
        imSquare = _mm256_mul_pd(currentIm, currentIm);

        currentIm = _mm256_mul_pd(currentIm, two);
        currentIm = _mm256_fmadd_pd(currentIm, currentReal, cIm);

        currentReal = _mm256_sub_pd(realSquare, imSquare);
        currentReal = _mm256_add_pd(currentReal, cReal);

        __m256i isSmaller = _mm256_castpd_si256(_mm256_cmp_pd(_mm256_add_pd(realSquare, imSquare), four, _CMP_LE_OS));
        result = _mm256_select_si256(_mm256_add_epi64(one, result), result, isSmaller);

        //if (i % 10 == 0 && !isSmaller.m256i_i64[0] && !isSmaller.m256i_i64[1] && !isSmaller.m256i_i64[2] && !isSmaller.m256i_i64[3]) return result;
    
    return result;


using namespace std;

int main() 
    MandelbrotSet m;
    std::complex<double> point(-0.14203954214360026, 1);

    __m256i result_avx = m.getIterationCount(_mm256_set_pd(-0.14203954214360026, -0.13995837669094691, -0.13787721123829355, -0.13579604578563975),
        _mm256_set_pd(1, 1, 1, 1), 2681);

    int result_normal = m.getIterationCount(point, 2681);
    cout << "Normal: " << result_normal << ", AVX: " << result_avx.m256i_i64[0] << ", at point " << point << endl;
    return 0;

当我运行这段代码时，我得到以下结果： （点 -0.14203954214360026 + i 是有意选择的，因为这两种方法在大多数点上返回相同/几乎相同的值）

Normal: 13, AVX: 20, at point (-0.14204,1)

1 的差异可能是可以接受的，但 7 的差异似乎很大，因为这两种方法都使用双精度。 AVX 指令的精度是否低于“普通”指令？如果不是，为什么两个结果差异如此之大？ 我使用 MS Visual Studio 2017、MS Visual C++ 2017 15.6 v14.13 141，我的电脑有一个 i7-7700K 处理器。该项目是为 x64 编译的。如果是没有优化或完全优化的编译器，结果是一样的。 渲染结果如下所示： AVX: 普通的

循环过程中realSquare和imSquare的值如下：

0, 0, 0
1, 0.0201752, 1
2, 1.25858, 0.512543
3, 0.364813, 0.367639
4, 0.0209861, 0.0715851
5, 0.0371096, 0.850972
6, 0.913748, 0.415495
7, 0.126888, 0.0539759
8, 0.00477863, 0.696364
9, 0.69493, 0.782567
10, 0.0527514, 0.225526
11, 0.0991077, 1.48388
12, 2.33115, 0.0542994
13, 4.5574, 0.0831971

在 AVX 循环中，值是：

0, 0, 0
1, 0.0184406, 1
2, 1.24848, 0.530578
3, 0.338851, 0.394109
4, 0.0365017, 0.0724287
5, 0.0294888, 0.804905
6, 0.830307, 0.478687
7, 0.04658, 0.0680608
8, 0.024736, 0.78746
9, 0.807339, 0.519651
10, 0.0230712, 0.0872787
11, 0.0400014, 0.828561
12, 0.854433, 0.404359
13, 0.0987707, 0.0308286
14, 0.00460416, 0.791455
15, 0.851277, 0.773114
16, 0.00332154, 0.387519
17, 0.270393, 1.14866
18, 1.02832, 0.0131355
19, 0.773319, 1.51892
20, 0.776852, 10.0336

【问题讨论】：

您的非 AVX 数学运算是否使用 x87 以 32 位代码完成，内部精度为 80 位？如果不是，则使用您正在使用的相同 AVX 指令的标量版本完成，例如vmulsd 而不是 vmulpd，它使用 IEEE754 64 位 double，包括不正常的支持。 您使用什么编译器和选项？

您是否尝试在标量循环中使用fma() 来匹配您在向量循环中执行的 FP 操作？我忘记了 MSVC 是否将 x*y + z 融合到 FMA 中，或者这取决于它是否等同于 -ffast-math。 （我从"stdafx.h" 认为您使用的是 MSVC，但这并不能告诉我您是否正在使用 x87 标量数学制作 32 位可执行文件。）相关：randomascii.wordpress.com/2012/03/21/… 关于中间 FP 精度以及 MSVC 的旧版本做了一些奇怪的事情来降低 x87 的精度。

@PeterCordes：我编辑了问题以包含信息。此外，反汇编使用mulsd 进行乘法运算。使用 fma() 并没有什么不同，即使标量部分没有编译为 fma 指令。

不要将图像放在 Dropbox 上。上传到 stack.imgur.com 并内联显示它们

不相关的代码审查：不要将向量常量放在全局范围内或将它们设为静态，例如const __m256i one = _mm256_set1_epi64x(1);。编译器实际上比在函数中定义它们更糟糕。 （它们在运行时通过从另一个向量常量复制来初始化静态存储。）另外，您不需要_mm256_select_si256。而是使用vcmpps 结果作为 0 / -1 整数。即result = _mm256_sub_epi64(result, _mm256_castpd_si256(cmp_result))。 x - (-1) 是 x++，x - (0) 是 x。另外，使用set1(2) instead of set(2,2,2,2)`。

【参考方案1】：

颠倒传递给_mm256_set_pd的参数顺序可以解决问题。

如果您在调试器中检查cReal 的值，您会看到第一个元素设置为-0.13579604578563975 而不是-0.14203954214360026。

【讨论】：

或者使用_mm256_setr_pd(0,1,2,3)而不是_mm256_set_pd(3,2,1,0)，如果您更喜欢以与数组初始值设定项相同的顺序编写向量，而不是按照英特尔的正常顺序使左移向左。

谢谢，对不起。我经常调试它并忽略了这一点。