我写了一个程序来显示 mandelbrot 集。为了加快速度，我通过 <immintrin.h> 使用了 AVX(实际上是 AVX2)指令。 header 。
问题是:AVX 计算( double )的结果有伪影，它与使用“正常” double 计算的结果不同。
详细来说，有一个函数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 的差异似乎很大，因为这两种方法都使用 double 。
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

最佳答案

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

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

关于c++ - AVX计算精度，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/51917634/