9

考虑以下程序：

for i=1 to 10000000 do
  z <- z*z + c

其中z和c是复数。

使用 x87 与 SSE 以及单精度与双精度算术的该程序的有效 x86 汇编器实现是什么？

编辑我知道我可以用另一种语言编写它并相信编译器会为我生成最佳机器代码，但我这样做是为了学习如何自己编写最佳 x86 汇编程序。我已经看过由 生成的代码gcc -O2，我的猜测是还有很大的改进空间，但我不够熟练，无法自己手动编写最佳的 x86 汇编程序，所以我在这里寻求帮助。

4

3 回答 3

7

您不需要在汇编程序本身中执行此操作- 您可以通过内在函数使用 SSE 来实现有效的实现，特别是如果您可以使用单精度。

temp.re = z.re * z.re - z.im * z.im;
temp.im = 2.0 * z.re * z.im;
z.re = temp.re + c.re;
z.im = temp.im + c.im;

如果您适当地打乱您的输入向量，那么您可以在一条指令中获得所有乘法 ( _mm_mul_ps)，并在第二条指令 ( _mm_hadd_ps) 中获得加法。

如果您需要双精度，则适用相同的一般原则，但您需要两次乘法和两次水平加法。

请注意，大多数现代 x86 CPU 都有两个标量 FPU，因此 SSE 中双精度的好处可能不值得 - 但是单精度绝对应该是一个胜利。

这是使用 SSE 的初始工作实现 - 我认为它现在或多或少已调试 - 性能并不比使用 gcc -O3 编译的标量代码好多少，因为 gcc 在为此生成 SSE 代码方面做得很好：

static Complex loop_simd(const Complex z0, const Complex c, const int n)
{
    __m128 vz = _mm_set_ps(z0.im, z0.re, z0.im, z0.re);
    const __m128 vc = _mm_set_ps(0.0f, 0.0f, c.im, c.re);
    const __m128 vs = _mm_set_ps(0.0f, 0.0f, -0.0f, 0.0f);
    Complex z[2];
    int i;

    for (i = 0; i < n; ++i)
    {
        __m128 vtemp;

        vtemp = _mm_shuffle_ps(vz, vz, 0x16); // temp = { z.re, z.im, z.im, z.re }
        vtemp = _mm_xor_ps(vtemp, vs);        // temp = { z.re, -z.im, z.im, z.re }
        vtemp = _mm_mul_ps(vtemp, vz);        // temp = { z.re * z.re, - z.im * z.im, z.re * z.im, z.im * z.re }
        vtemp = _mm_hadd_ps(vtemp, vtemp);    // temp = { z.re * z.re - z.im * z.im, 2 * z.re * z.im, ... }
        vz = _mm_add_ps(vtemp, vc);           // temp = { z.re * z.re - z.im * z.im + c.re, 2 * z.re * z.im + c.im, ... }
    }
    _mm_storeu_ps(&z[0].re, vz);
    return z[0];
}

请注意，内部循环只有 6 条 SSE 指令（实际上应该是 5 条）+ 循环本身的一些内务管理：

L4:
    movaps  %xmm0, %xmm1
    shufps  $22, %xmm0, %xmm1
    xorps   %xmm3, %xmm1
    mulps   %xmm1, %xmm0
    haddps  %xmm0, %xmm0
    addps   %xmm2, %xmm0
    incl    %eax
    cmpl    %edi, %eax
    jne L4
L2:
于 2012-04-26T08:43:08.707 回答
6

查看您最喜欢的编译器的反汇编。如果您希望对z和的多个值执行此计算c（例如计算 mandelbrot 图像），我建议您一次处理四个值并将它们放入 SSE 寄存器中。如果您查看 Paul R 答案中的代码，您可以一次对四个值进行所有这些计算：

__m128 z_im, z_re, c_im, c_re; //Four z and c values packed
__m128 re = _mm_sub_ps(_mm_mul_ps(z_re, z_re), _mm_mul_ps(z_im, z_im));
__m128 im = _mm_mul_ps(z_re, z_im);
im = _mm_add_ps(im, im); // Multiply by two
z_re = _mm_add_ps(re, c_re);
z_im = _mm_add_ps(im, c_im);
于 2012-04-26T08:43:23.613 回答
3

Z = Z*Z + C

这就是曼德布洛特分形迭代。

我相信你会在网上找到高度优化的代码。我将从 Xaos 和 Fractint 的源代码开始。

于 2012-04-26T08:52:54.653 回答