使用 SSE/SSE2 生成分形

• 下载源代码 - 39.2 Kb

引言

你可能以前听说过分形。它们是像上面展示的那样美丽的图像。任何分形都可以用迭代公式来描述。所以你可以通过评估这些公式并确定每个像素的颜色来生成分形。这是一个巨大的计算任务，绘制分形需要一台快速的 CPU 和一个精心优化的程序。

分形及其性质是在非线性动力学中研究的。除了科学家，还有一些人喜欢探索分形，仅仅因为这些图像优雅而奇特。

本文将跳过数学背景。你可以在这里找到必要的解释

• Julia 和 Mandelbrot 集合，作者：David E. Joyce。
• Mandelbrot 集合是什么以及 Fabio Cesari 关于分形的其他文章。
• 你的分形，先生 - 谢谢你，James，作者：q123456789。
• Julia 和 Mandelbrot 集合 - 理论，作者：Martin Pfingstl。
• 分形理论，作者：Cygnus Software。

最后三篇文章还提供了一些 C 语言编写的分形生成代码。如果你是分形新手，我建议你阅读 Fabio Cesari 的文章。他的解释简单明了，配有表格和图像，胜过千言万语。

我的目标和几种分形生成器的评测

我需要找到或开发教育软件来学习分形及其性质。该软件将用于向学生展示 Mandelbrot 和 Julia 集合是什么，并给他们一些练习。学生预计只使用该程序一个课时，因此它需要简单直观。主要**功能**如下

• 使用快速的 SSE 驱动算法生成 Julia 和 Mandelbrot 集合；
• 将**函数的轨道**（即迭代序列）显示为图表和表格；
• 显示任何点的轴和坐标；
• 轻松缩放（左键单击放大，右键单击缩小或类似操作）；
• 生成与 Mandelbrot 集合中某个点对应的 C 参数的 Julia 集合。

最后一点意味着你可以从 Mandelbrot 集合的图像中选择 Julia 集合的 C 值。你只需按住鼠标按钮，将光标移到 Mandelbrot 集合上，Julia 集合就会随着你的操作而改变。

简称，我将称之为**画中画模式**（Picture iN Picture）。

再谈谈**非目标**，即我将不实现的功能

• 花哨的颜色和调色板（如果我这样做，学生就会专注于着色方案而不是数学）；
• 除 Julia 和 Mandelbrot 以外的任何类型的分形；
• 支持 3DNow（旧 AMD 处理器中的 SIMD 扩展；现代 Athlon XP 使用 SSE）；
• 移植到其他操作系统和处理器（x86 和 Win9x/XP 对我来说已经足够了）；
• 多处理器支持。

我的第一个尝试是Chaos Pro。它是一个复杂的免费应用程序，可以生成 Julia 和 Mandelbrot 集合以及许多其他分形。你甚至可以使用内置的编程语言添加自己的分形公式。虽然该程序能够创建炫酷的 3D 图像和动画，使用不同的调色板等，但它仍然缺少一些基本功能。首先，你无法查看轨道和轴。其次，Mandelbrot 集合在我 1.5 GHz 的 Pentium M 笔记本电脑上大约需要 1.6 秒生成，而 C = -0.12+0.74i 的 Julia 集合需要 1.4 秒绘制。Chaos Pro 支持画中画模式，但使用起来不方便，因为渲染速度太慢。分形无法完全实时渲染（当你移动光标时，图像不会立即改变）。第三（这也是我的主要原因），我经常遇到“访问违例”错误或其他与 Chaos Pro 相关的 bug。当然，对于一个业余爱好者编写的免费应用程序来说，这是可以接受的，但我想要更可靠的东西。

然后我偶然发现了FracInt，一个很老牌的 DOS 分形绘图程序。它已经开发了好几年，现在有一个 20.0 版本可用。该程序相当复杂和强大，具有许多有趣的特性，例如能够查看数十种分形类型、播放分形音乐、使用画中画模式以及创建立体图像。它的主要缺点是过时的用户界面。缩放需要很长时间：你拖动鼠标设置所需的矩形大小，然后将其移动到目标点，最后点击两次（左键单击放大，右键单击缩小）。这对我来说是不可接受的，因为普通学生将花一半的时间来理解如何使用这个缩放界面。支持高达 1600x1200 的视频模式（1024x768 VESA 应该适用于大多数显卡），但当我在 Windows XP 中运行 FracInt 并按下 Alt+Tab 时，我的视频驱动程序崩溃了。输出到声卡不起作用，我不得不使用我的 PC 扬声器（还有人记得这个臭名昭著的尖叫声吗？;）。更不用说这个旧应用程序缺少 MMX 或 SSE 支持了。

经过几个小时的网上冲浪，我发现了 Fast Floating Fractal Fun，也叫FFFF。它是一个非常快速的 Mandelbrot 集合生成器，支持 SSE、SSE2、3DNow 和 GPU 顶点程序。如果我早知道它，就不会开始开发自己的分形渲染引擎了 :)。不过，除了缩放和在不同计算引擎之间切换外，FFFF 没有太多功能，所以如果你决定使用它，你需要添加一个 Julia 集合生成器并开发一个更精致的 GUI。

算法

以下是 Mandelbrot 集合生成算法的伪代码

complex P, Z;
for each point P on the complex plane {
   Z = 0;
   for i = 0 to ITER {
      if abs(Z) > 4.0
         set point color to palette[i];
         break the inner loop and go to the next point;
      Z = Z * Z + P;
   }
   set point color to black;
}

P 和 Z 是复数变量。你可以在上述文章中找到关于复数的通用解释。

对于图像的每个像素 P，循环会重复 ITER 次。因此，**执行时间**为 O(X*Y*ITER)，其中 ITER 是迭代次数，图像的尺寸为 X x Y 像素。

abs(Z) 是一个计算复数绝对值的函数。想象一个半径为 2.0、中心在原点的圆。如果点 Z 在圆外，abs(Z) 将大于 4.0。在这种情况下，内循环应该终止，因为你知道 Mandelbrot 集合的每个点都位于圆内。如果 Z 超出了圆，那么点 P **就不属于 Mandelbrot 集合**。我们用一个颜色来标记这个点，该颜色显示了它在迭代 Z = Z * Z + P 计算时逃离圆的速度。完整的迭代次数 i 显示了这个速度。因此，你从某个调色板中获取第 i 种颜色，然后继续处理下一个点。调色板就是一个颜色数组。例如，在第一次迭代后超出圆的点被标记为深绿色；第二次迭代后离开圆的点被标记为浅绿色，依此类推。

否则，如果循环执行了 ITER 次并且点 Z 从未超出圆，那么**该点属于 Mandelbrot 集合**。我们将其涂成黑色。

内循环重复计算 Z<SUB>N+1</SUB> = Z<SUB>N</SUB>^2 + P，这是 Mandelbrot 集合的迭代公式。Z 初始化为零，但很容易看出**可以将 Z 设置为 P**，因为 Z<SUB>1</SUB> = Z<SUB>0</SUB>^2 + P = 0^2 + P = P。

使用调色板是为 Mandelbrot 集合着色的最简单、最快的方法，但最终图像只有 ITER+1 种颜色。有一些方法可以渲染 24 位图像。例如，ChaosPro 使用这个公式

color_index = (i - log(log(abs(Z))) * 1.442695 + 2) / 252.0;

系数 252.0 和 2.0 将 color_index 值缩放到 0..1 的范围。然后，color_index 被转换为 24 位颜色。不幸的是，计算对数是一个缓慢的操作，它会使整体执行时间增加大约 30-50%。这就是为什么我为了速度牺牲了 24 位颜色，而选择了有限的调色板。

要将算法转化为 **C 代码**，你需要记住复数乘法和计算其绝对值的公式（我希望你没有错过大学的数学课！）。

void PaintMandelbrotFPU(int w, int h, long* b, double left, double top,
                        double right, double bottom) {
   int x, y, i; double x2, y2, px, py, zx, zy, dx, dy;
   dx = (right - left) / w; // Find horizontal and vertical increase for each pixel
   dy = (bottom - top) / h;
   py = top;
   for(y = 0; y < h; ++y) {
      px = left;
      for(x = 0; x < w; ++x) {
         zx = px, zy = py;
         for(i = 0; i < ITER; ++i) {
            x2 = zx * zx,   y2 = zy * zy;
            if(x2 + y2 > 4.0)
               break;
            zy = zx * zy * 2 + py; // Calculate z = z * z + p
            zx = x2 - y2 + px;
         }
         *b++ = a[i];
         px += dx; // Move to the next pixel
      }
      py += dy; // Move to the next line of pixels
   }
}

复数变量 p 变成了 px（实部）和 py（虚部）。对于 zx 和 zy 也是如此。zx 和 zy 的平方在两个地方被使用：在计算绝对值和在计算 z * z + p 时。这就是为什么我为它们创建了两个变量 x2 和 y2。变量 zx 和 zy 被初始化为 px 和 py 而不是零，因为这样你可以通过跳过 0*0 乘法来节省时间。

请注意整数变量 x 和 y，它们计算当前像素的坐标，以及一对 (px, py)，它们是复平面上的坐标。px 和 py 的增加独立于 x 和 y。这极大地提高了性能，因为你不必在每次迭代循环时都将屏幕坐标 (x, y) 转换为复平面坐标 (px, py)。这种性能提升的代价是**精度的损失**。浮点数精度有限，并且在多次加法 px += dx, py += dy 过程中会累积误差。这会导致一个令人不快的效应：当你放大图像时，它会“抖动”。另一种解决方案可能是每次迭代都计算 px = (right - left) * x / w, py = (bottom - top) * y / h，但这似乎不切实际，因为除法延迟很长。

最后，a 是我的调色板，b 是函数生成的**位图**。位图 b 稍后使用 DirectDraw 或 GDI 函数进行绘制。这个过程的细节对我们来说现在不重要。只需记住 b 的大小是 w * h 像素，你需要用调色板中的颜色填充它。

当内循环使用 break 终止时，i 的值小于 ITER。当它未终止时，循环后 i == ITER。利用这个属性，我做了一个巧妙的技巧，创建了数组 a[ITER+1]，并将这个数组的最后一个元素设置为黑色。循环结束后，a[i] 在 a[ITER] 中引用黑色。更严谨的做法是

// ...
for(i = 0; i < ITER; ++i) {
    x2 = zx * zx,   y2 = zy * zy;
    if(x2 + y2 > 4.0) {
       *b++ = a[i];
       px += dx;
       goto NEXTPIXEL;
    }
    zy = zx * zy * 2 + py; // Calculate z = z * z + p
    zx = x2 - y2 + px;
}
*b++ = RGB(0,0,0);
px += dx; // Move to the next pixel
NEXTPIXEL:
// ...

Julia 集合的生成算法与 Mandelbrot 集合的算法非常相似，所以我不会在这里展示。主要区别在于引入了复数参数 (PX, PY)。

MOV-ADD 方法

我做了一些性能剖析，发现 PaintMandelbrotFPU 函数占用了整个执行时间的 95%-99%。绘制只占 GDI 的 5% 或 DirectDraw 的 1%。因此，优化这个函数比程序中的任何其他部分都重要。

进行此类优化的最佳方法之一是使用 **SSE**（流式 SIMD 扩展）。这些指令允许你同时计算四个像素值。所以，至少在理论上，速度会提高 4 倍（实际数字约为 3.2-3.3）。你可以在 Intel 架构软件开发者手册第 10 章和 AMD64 架构程序员手册第 1 卷第 4 章中找到更多关于 SSE 的信息。在接下来的解释中，我将假设你了解 SSE 和 x86 汇编编程的基础知识。

有许多不同的方法可以为这个函数安排 SSE 指令，其中一些方法比其他方法更快。以下是针对 Pentium IV 处理器的一种非常快速的序列。我称之为 **MOV-ADD** 方法。

; xmm0 = zx  xmm1 = zy  xmm2 = tmp  xmm3 = tmp  
;            xmm4 = px  xmm5 = py  xmm6 = result  xmm7 = 4.0
; eax = tmp    ecx = loop counter
   MOVAPS xmm0,xmm4
   XORPS  xmm6,xmm6
   MOVAPS xmm1,xmm5
   mov ecx, ITER
ILOOP:
   ; xmm0 = zx             xmm1 = zy
   MOVAPS xmm2, xmm0
   MULPS  xmm0, xmm0
   MOVAPS xmm3, xmm1
   ADDPS  xmm1, xmm1
   ; xmm0 = zx^2           xmm1 = 2 * zy     xmm2 = zx           xmm3 = zy
   MULPS  xmm1, xmm2
   MOVAPS xmm2, xmm0
   MULPS  xmm3, xmm3
   ; xmm0 = zx^2           xmm1 = 2*zy*zx    xmm2 = zx^2         xmm3 = zy^2
   ADDPS  xmm1, xmm5
   SUBPS  xmm0, xmm3
   ADDPS  xmm2, xmm3
   ; xmm0 = zx^2 - zy^2    xmm1=2*zy*zx+py   xmm2 = zx^2 + zy^2  xmm3 = zy^2
   CMPLEPS xmm2, xmm7
   ADDPS   xmm0, xmm4
   MOVMSKPS eax, xmm2
   test eax, eax
   jz EXIT
   ANDPS   xmm2, xmm7      ; xmm6 += (xmm2 < 4.0) ? 4.0 : 0.0;
   ADDPS   xmm6, xmm2
   sub ecx, 1
   jnz ILOOP
EXIT:

这里只显示了内循环。循环开始时，xmm0 和 xmm1 分别包含 px 和 py，这是 zx 和 zy 的初始值。寄存器 xmm4 和 xmm5 存储 px 和 py（它们都是循环不变量）。寄存器 xmm7 在所有双字中都包含常量 4.0。

程序开始时，将 xmm0 和 xmm1 保存到临时寄存器 xmm2 和 xmm3。然后将 zx 平方，并将 zy 加倍（将 xmm1 加到自身）。其他计算步骤在上面的注释中显示。

寄存器 xmm6 存储**完整迭代次数**。请注意，xmm6 中的数字应该为每个 4 个像素独立计数。CMPLEPS 指令将一个**位掩码**放入寄存器 xmm2，而序列 MovMskPS-test-jz 检查是否所有 4 个像素都超出了半径为 2.0 的圆。如果发生这种情况，循环将终止。在其他情况下，程序应该为仍在圆内的像素增加完整迭代次数。然后将掩码与 4.0 进行 AND 运算，以在对应于这些像素的双字中得到 4.0，在其他双字中得到 0.0。这个值被加到 xmm6。简而言之，当循环结束时，xmm6 将包含每个 4 个像素的 4.0 * 完整迭代次数。该值稍后将用作我们调色板的偏移量。

现在，让我们计算 **Pentium IV 的执行时间**并找出代码中的瓶颈。Agner Fog 的手册描述了如何执行此类分析；如果你是 Pentium IV 优化新手，请阅读该手册。另请参阅 Intel 的官方文档：Intel Architecture Optimization Reference Manual 和 Microarchitecture of the Pentium 4 Processor。

首先，内存访问不是这里的瓶颈，因为计算比写入内存花费的时间多得多。对于通常的 1024x768 分辨率，程序写入 1024*768*4 字节 = 3 Mb。现代 DDR 内存的典型吞吐量可能约为 1-2 Gb/s。粗略估计表明，内存访问需要 1.5-3 毫秒，而代码的实际执行时间为 70-80 毫秒。对代码的分析也得出了相同的结果：内存只访问 1 个像素两次（见下一章代码），但内循环执行 ITER 次（通常约为 64 次），并且它包含了具有长延迟的 SSE 指令。

因此，用于优化缓存性能的各种技巧对这段代码毫无用处。指令解码也是如此：循环很小，完全适合跟踪缓存，并且只解码一次。真正的瓶颈是 SSE 指令的**长延迟和有限吞吐量**。

为简单起见，我们假设管道最初是空的，并且整个代码都在跟踪缓存中。寄存器名称 xmm1、xmm2 等是 x1、x2 的缩写，以适应图表。

前三个指令（MOVAPS x2, x0、MULPS x0, x0 和 MOVAPS x3, x1）从跟踪缓存发送到操作队列。FP/SSE 执行集群有两个端口：一个用于 add、mul、div 操作，在 Agner Fog 的手册中称为 **fp** 单元，另一个用于寄存器之间以及存储的操作（**mov** 执行单元）。这两个单元可以同时开始执行 MOVAPS 和 MULPS 指令。

在第二个时钟周期，接下来的三个指令进入队列。其中，MOVAPS x3, x1 和 ADDPS x1, x1 可以并行执行。MOVAPS 指令的倒数吞吐量为 1，因此第二个指令可以在这个周期开始。

MOVAPS 指令的延迟为 6 个时钟周期，而 MULPS x1, x2 依赖于 MOVAPS x2, x0 的结果（**依赖关系**用图表上的箭头表示）。所以 MULPS 可以在第六个周期开始（此时，队列应该已经充满了等待执行的操作）。

下一个指令 MOVAPS x2, x0 依赖于 MULPS x0, x0。MULPS 指令还有一个**额外延迟**为 1，这意味着如果下一个操作进入不同的执行单元，它需要等待 1 个额外的时钟周期。MULPS 和 MOVAPS 来自不同的执行单元 **fp** 和 **mov**，我必须在它们之间插入一个周期的延迟。

现在，指令 MULPS x3, x3。乘法子单元（Intel 称为 FP 乘法器）每两个周期可以执行一次乘法，因此第二个 MULPS 指令必须等待一个周期。它将在第八个周期开始。

接下来的五个指令进入 FP 加法子单元（**FP adder**）。它们都有 4 个时钟周期的延迟。ADDPS x1, x5 可以在第十二个周期开始，此时 MUL x1, x2 乘法的结果已经准备好。SUBPS x0, x3 也是如此；它在第十四个周期开始。指令 ADDPS x2, x3 可以在第十六个周期之前开始，但 FP 加法的吞吐量将同时运行的指令数量限制为 2。序列 SUBPS x0, x3; ADDPS x2, x3; ADDPS x0, x4; CMPPS x2, x7 完美地满足了 FP 加法的这一要求。两个相邻的指令来自不同的依赖链，因此它们是完全独立的，可以并行执行。

我们的目标是尽快完成 x0 = x0 ^ 2 - x1 ^ 2 + x4, x1 = 2 * x0 * x1 + x5 的计算。在图表中，计算**x0 和 x1 新值**的指令用不同深浅的灰色标记，非相关指令用白色标记。你需要比其他指令更早地调度“灰色”指令，因为循环的下一次迭代在没有 x0 和 x1 新值的情况下无法开始。前一次迭代中的 ANDPS 和 MOVMSKPS 指令可以与下一次迭代并行执行，前提是分支预测正确。因此，你可以将它们安排在灰色指令之后。

我的时间测量证实了这些考虑。现在 MOVAPS x2, x0 和 ADDPS x2, x3 之间有一个很长的延迟（见图表）。但如果你尝试将 ADDPS x2, x3 放在 SUBPS x0, x3 之前，总执行时间会增加。这意味着只要 x0 和 x1 尽早计算出来，延迟就不重要。

你看到代码的速度取决于指令延迟和吞吐量，因为有**三个相当独立的链**可以同时执行。**关键路径**用深灰色标记。MOVAPS 增加了不必要的延迟（也许，如果 x86 体系结构有像 RISC 那样的三地址指令，代码会执行得更快，谁知道呢）。子单元没有得到充分利用。在循环的第一阶段（周期 0-12），MOV 单元和 FP 乘法器承受了很大的负载；之后，在周期 14-22，FP 加法器过载，FP 乘法器根本没有使用。你可以从两个迭代中提取指令并重新排序以获得更均衡的负载，但这需要更多的 SSE 寄存器。使用 AMD64（或 Intel EM64T），可以使用 8 个额外的寄存器，因此代码**可以在 64 位模式下进一步优化**。

绘制图表后，**我按启动时间对指令进行了排序**。例如，MOVAPS x2, x0 和 MULPS x0, x0 将是最早启动的指令，所以我将它们放在代码的最前面。然后是 MOVAPS x3, x1 和 ADDPS x1, x1 开始，它们是源代码中的下一条指令。对指令进行排序可以减少操作队列中等待的操作数量。现在，你的处理器无需在运行时重新排序指令：它们已经根据其延迟、吞吐量和依赖关系进行了排序。

但这段代码中的**分支预测**怎么样？它会成为瓶颈吗？事实上，错误预测分支的惩罚会花费程序大量执行时间。但完整迭代次数是我们想要得到的结果。你无法准确预测循环会执行多少次，因为这个值本身就是程序的输出。如果你确切地知道循环重复了多少次，那你为什么还要运行这个程序呢？:)

我随机选择了两个集合的放大图像，然后计算了它们以及一张大尺度图像的**迭代次数分布**。

蓝色条表示 Mandelbrot 集合的大尺度图像；红色和黄色条表示两个放大图像。我的程序中 ITER 等于 64，所以图表显示了从 0 到 5 次迭代的点频率，从 6 到 11 次迭代，依此类推，直到 60 次以上迭代。显然，**最小和最大迭代次数**构成了图像的最大部分。对于蓝色条，最小值为 0-5 次迭代（占所有像素的 59%）；对于红色条，为 12-17 次（43%）；对于黄色条，为 30-35 次（8%）和 36-41 次（15%）。最小迭代次数对应于被一种颜色填充的巨大背景区域，而最大迭代次数则表示属于该集合的黑色点。因此，在大图像中 26% 的像素属于 Mandelbrot 集合，而在两个放大图像中，分别有 21% 和 62% 的点属于该集合。

从这个分析可以看出，如果循环终止，那么很可能它会在最小迭代次数时终止。但如果循环没有早期终止，它通常就不会终止。这个统计数据很有趣，但我无法在代码中使用它。也许你，读者，可以做到。任何建议都将受到赞赏。

让我们对这种 MOV-ADD 方法做一些**结论**

• 总体策略是尽可能早地计算下一迭代将需要的值（在本例中是 x0 和 x1）。这些计算构成了程序的关键路径，而总体执行时间取决于它们的完成时间。
• 通过调度独立的指令并行执行，我试图提高指令级并行性（**ILP**）。例如，前四条指令并行执行加法、乘法和两次移动。
• 指令按它们的启动时间排序，这有助于减小操作队列的长度。

FFFF 方法

上面提到的 Mandelbrot 集合生成器 **FFFF** v3.2.2 使用另一种方法。FFFF 是开源的，因此在此处发布其代码是完全合法的。

; xmm0 = zx  xmm1 = zy  xmm2 = tmp  xmm3 = tmp  xmm4 = px  
;            xmm5 = py  xmm6 = result  xmm7 = 4.0
; eax = tmp    ecx = loop counter
   MOVAPS xmm0,xmm4
   XORPS  xmm6,xmm6
   MOVAPS xmm1,xmm5
   mov ecx, ITER
ILOOP:
   ; xmm0=zx      xmm1=zy
   MOVAPS xmm2, xmm0
   MULPS  xmm2, xmm1
   MULPS  xmm0, xmm0
   MULPS  xmm1, xmm1
   ; xmm0=zx^2    xmm1=zy^2  xmm2=zx*zy 
   addps xmm2, xmm2
   movaps xmm3, xmm0
   addps xmm3, xmm1
   ; xmm0=zx^2    xmm1=zy^2  xmm2=2*zx*zy  xmm3=zx^2+zy^2
   cmpltps xmm3, xmm7
   movmskps eax, xmm3
   test eax, eax
   jz EXIT

   subps xmm0, xmm1
   movaps xmm1, xmm2
   ; xmm0=zx^2-zy^2  xmm1=2*zx*zy   xmm3=zx^2+zy^2
   ANDPS  xmm3, xmm7 ; xmm6 += (xmm3<radius) ? 4.0 : 0.0
   ADDPS  xmm6, xmm3
   ADDPS  xmm0, xmm4
   ADDPS  xmm1, xmm5
   dec ecx
   jnz ILOOP
EXIT:

FFFF 尝试将 x0^2 + x1^2 与 4.0 进行比较，只有在比较失败时，它才找到 x0 和 x1 的新值。如果比较返回 true，循环将终止，**无需计算 x0 和 x1**。虽然这似乎是一个好策略，但我的实验表明并非如此。看这个长链：MOVAPS xmm3, xmm0; ADDPS xmm3, xmm1; CMPLTPS xmm3, xmm7; MOVMSKPS eax, xmm3。当这条链运行时，有很多机会可以利用空闲的执行单元并并行运行第二条链。但 FFFF 作者没有在该链的指令之间插入任何指令。调度程序必须自己重新排序指令。有时它会失败，执行单元会空闲。例如，MULPS xmm1, xmm1 和 SUBPS xmm0, xmm1 之间有很大的距离，以至于第二条指令可能没有及时发送到队列。在这种情况下，应该在这些指令之间添加额外的延迟。同样对于 ADDPS xmm2, xmm2 和 MOVAPS xmm1, xmm2：它们有很小的几率在没有延迟的情况下执行。

图表显示了理想情况，没有任何延迟。即使如此，x0 的计算需要 20 个时钟周期，x1 需要 27 个。对于我的 MOV-ADD 方法，分别为 16 和 22 个时钟周期。关键路径（用深灰色标记）包括**两个具有长延迟的 MOVAPS 指令**，而 MOV-ADD 方法在其关键路径上只使用了一个 MOVAPS。

在 MOVAPS x1, x2 和 MOVAPS x3, x0 之前有**额外的延迟**，因为前一条指令在不同的单元中执行。由于 FP 加法的吞吐量有限，ADDPS x3, x1 被推迟了一个周期。

在其他方面，MOV-ADD 与 FFFF 相似。两种方法都使用 SSE 同时计算 4 个像素值，都具有相同的精度限制和相同的缩放抖动效应。我使用了类似的算法，但以不同的方式安排了指令以提高性能。FFFF 与我的程序之间的其他区别将在后面展示。

在上面的 FFFF 代码中，寄存器名称已更改为 MOV-ADD 中使用的名称（xmm0 用于 zx，xmm1 用于 zy，依此类推）。此更改不影响性能，但使方法的代码更易于比较。

Vodnaya 方法

我还开发了 **Vodnaya**，一种针对 Pentium M 优化的方法。这个名字完全没有意义，所以如果你觉得说俄语单词困难，可以称它为 Vod-method 或 V-method ;)。它在 Pentium M 和 Pentium III 上都优于 MOV-ADD 和 FFFF 方法。关于 Pentium M 的微体系结构只知道很少，所以我无法分析延迟和吞吐量。经过大量的试错，我发现了一系列指令，它们在这些处理器上似乎是最快的。该方法在指令数量上也名列前茅：比 MOV-ADD 或 FFFF 少一条指令。

   MOVAPS xmm0,xmm4
   XORPS  xmm6,xmm6
   MOVAPS xmm1,xmm5
   mov ecx, ITER
ILOOP:
   ; xmm0 = zx             xmm1 = zy
   MOVAPS xmm2, xmm1
   MULPS  xmm2, xmm2
   MULPS  xmm1, xmm0
   ; xmm0 = zx             xmm1 = zy*zx           xmm2 = zy^2
   MOVAPS xmm3, xmm2 ; save zy^2
   MULPS  xmm0, xmm0
   ADDPS  xmm2, xmm0
   ; xmm0 = zx^2           xmm1 = zy*zx           xmm2 = zx^2+zy^2     xmm3 = zy^2
   CMPLEPS  xmm2, xmm7
   ADDPS  xmm1, xmm1
   SUBPS  xmm0, xmm3
   ; xmm0 = zx^2-zy^2      xmm1 = zy*zx*2         xmm2 = zx^2+zy^2
   ADDPS  xmm1, xmm5
   ADDPS  xmm0, xmm4
   MOVMSKPS eax, xmm2
   test eax, eax
   jz EXIT
   ANDPS  xmm2, xmm7
   ADDPS  xmm6, xmm2
   sub ecx, 1
   jnz ILOOP
EXIT:

zx^2+zy^2 的和在临时寄存器 xmm2 中计算。zy^2 被复制到这里，以防止在程序计算 zx^2-zy^2 时被覆盖。zx 和 zy 的新值在原地计算，在 xmm0 和 xmm1 寄存器中。我猜整个方法都相当优雅。不幸的是，它不是 Pentium-IV 和 Athlon XP 处理器中最快的方法。

大战：MOV-ADD vs. FFFF vs. Vodnaya

现在，这三种方法的速度对比

在测试集中，图像尺寸为 1024x768 像素；使用 RDTSC 指令测量时间（参见 _SSEroutines_test.asm_ 文件中的源代码）。我重复了五次测试，并取了每种方法的最小值。图表中的时间以百万个时钟周期表示。较低的时间更好。处理器名称后面显示了它们的 CPUID 签名。例如，**Celeron F13** 是一个家族 15、型号 1、步进 3 的处理器，实际上是精简版的 Pentium IV。它的缓存比 Pentium IV 小，但核心是相同的，因此这款 Celeron 的结果应该与 Pentium IV 的结果几乎相等（抱歉，我无法获得 Pentium IV 进行测试）。“兔子”表示 Julia 集合 C = -0.12+0.74i。

显然，MOV-ADD 是 Pentium IV 和 Athlon XP 处理器的最佳方法。Vodnaya 在 Intel 的 Family 6 处理器（尤其是 Pentium III 和 Pentium M）上表现良好。FFFF 在任何类型的处理器上都绝对是输家。MOV-ADD 相对于 FFFF 的速度提升在 Pentium IV 上约为 12%，Vodnaya 相对于 FFFF 在 Pentium M 和 Pentium III 上为 7-11%。

现成的程序使用**Intel 的 Family 6 处理器上的 Vodnaya** 方法，在其他处理器上使用 MOV-ADD 方法。该方法使用 CPUID 指令检索的处理器签名在运行时进行选择。

CvtTss2si

计算完成后，程序应将结果转换为颜色值。C 语言中的单行代码（*b++ = a[i];）可以完成转换，但在汇编语言中，你需要做更多的工作。

程序将 xmm6 初始化为零，并在每次迭代时向其添加 4.0（参见上面任何方法的代码）。如果内循环未终止，它将重复 ITER 次。循环结束后，如果点不属于 Mandelbrot 集合，则完整迭代次数乘以 4 存储在 xmm6 寄存器中，否则其中包含 ITER * 4。还可以进行**相同的数组最后一个元素的技巧**：向数组添加一个元素，并将其值设置为黑色。这样，你就可以使用一种着色代码来处理属于 Mandelbrot 集合的点和不属于该集合的点。FFFF v3.2.2 使用条件移动来处理这两种类型的点（这是一种较慢的方式，并且 CMOV 指令不受某些旧 AMD 处理器的支持）。

请注意，xmm6 寄存器保存**乘以 4 的完整迭代次数**。FFFF 会对其进行除法以获取调色板中的索引。对于除法，它使用 shr 指令，该指令在模型号低于 3 的 Pentium IV 处理器上很慢。但是，一个调色板元素在内存中占 4 个字节。因此 FFFF 将数字除以 4，然后乘以 4 来获取内存偏移量。多么浪费时间！:) 我做得更简单

CVTTSS2SI eax, xmm6
; xmm6 contains an offset into palette, i.e. 4*number-of-iterations

mov eax, [esi + eax]  ; palette is started at [esi]
mov [edi], eax        ; the current pixel of the bitmap is at [edi]

CvtTss2si 指令**转换**（截断）一个单精度浮点值到 eax 寄存器中的整数值。迭代次数显然是整数。它将无损地转换，因此舍入模式在此处无关紧要。下一条指令从调色板中提取颜色，最后一条指令将颜色写入位图。位图可以位于视频内存或 RAM 中（现在不用管它）。

还有**其他指令**用于转换浮点数（即 CvtPs2Pi 和 CvtPs2Dq），但没有一个将结果写入通用寄存器。你必须将数字存储到内存中并立即从同一位置读取，这可能会在 Pentium IV 上导致额外的延迟。我用 CvtPs2Dq 测试了代码，发现它并不比使用 CvtTss2si 的代码快。

FASM 宏的强大之处

为 Mandelbrot 集合和 Julia 集合、黑白图像和彩色图像、SSE 和 SSE2 编写几种相同代码的变体需要大量例行工作。值得庆幸的是，大多数汇编器都有一些特殊的工具可以完成这项工作。这些工具是**宏指令**和条件汇编。宏允许你使用不同的参数重复指令序列，而条件汇编则允许根据这些参数插入指令。

我在几个汇编器中选择了 FASM 1.62，因为它体积小、速度快，并且具有丰富的宏语法。FASM 手册详细描述了语法。以下是 Vodnaya 方法中的一小段，用于说明我的程序中宏的使用。

macro JuliaMandelPaint color, type {
; ...
   ADDPS  xmm1, xmm1
   SUBPS  xmm0, xmm3
if type = Julia
   ADDPS  xmm1, dqword[cy1]
   ADDPS  xmm0, dqword[cx1]
else if type = Mandel
   ADDPS  xmm1, xmm5
   ADDPS  xmm0, xmm4
end if
; ...
}

; later in the code:
   mov eax, [fractaltype]
   test eax, eax
   jz MANDELBROT
JULIA:
   JuliaMandelPaint 1, Julia
   ret
MANDELBROT:
   JuliaMandelPaint 1, Mandel
   ret

如果参数 type 等于 Julia，FASM 将生成将 cx1 加到 xmm0 和 cy1 加到 xmm1 的代码，否则它将使用寄存器 xmm4-xmm5 中的值。这两种代码变体都将生成：第一种在标签 JULIA 处，第二种在 MANDELBROT 处。你只需检查变量 fractaltype 并决定使用哪种变体。

你可能会说有一种更简单的方法

   mov eax, [fractaltype]
   test eax, eax
   jz MANDELBROT
   ADDPS  xmm1, dqword[cy1]
   ADDPS  xmm0, dqword[cx1]
   jmp NEXT
MANDELBROT:
   ADDPS  xmm1, xmm5
   ADDPS  xmm0, xmm4
NEXT:

为什么不使用这段代码？因为检查位于深度嵌套的循环中，比较会重复数百万次，而只需要进行一次。更不用说如果条件跳转被错误预测可能发生的延迟了……因此，将比较移出循环是有利的。宏可以帮助简短地编写。

下一个有趣的任务是将方法改编为 **SSE2**。SSE2 指令并行处理两个双精度值。在其他方面，它们与 SSE 没有区别。即使是 SSE2 的延迟和吞吐量也与 SSE 相同。

所以你必须进行大量的搜索和替换，将 ADDPS 更改为 ADDPD，MULPS 更改为 MULPD，MOVAPS 更改为 MOVAPD 等。然后调整指针增量，得到两个函数：一个用于 SSE，一个用于 SSE2。这种方法用于 FFFF 和其他一些程序。它有一个严重的缺点。当你想更新程序时，你必须在两个函数中更改相同的代码。这很乏味且令人沮丧。有时你纠正了一个函数，但忘记更改另一个。

这个想法是创建一个名为 MOVAP 的宏，它**将被替换为 SSE 的 MOVAPS 和 SSE2 的 MOVAPD**。现在，你只需编写一个函数而不是两个函数，并在其中使用 MOVAP、MULP 和 ADDP 宏。在编译过程中，实际的 SSE 或 SSE2 指令将被替换为这些宏。

macro MOVAP a, b {
   if ssewidth = 4
      MOVAPS a, b
   else
      MOVAPD a, b
   end if
}
macro ADDP a, b {
   if ssewidth = 4
      ADDPS a, b
   else
      ADDPD a, b
   end if
}
macro MULP a, b {
   if ssewidth = 4
      MULPS a, b
   else
      MULPD a, b
   end if
}
; another 2 pages of such clumsy code
; for every SSE/SSE2 instruction

别急着写这段代码！有一个更优雅的解决方案。

macro defineinstr instr, [args] {
   common
   if ssewidth = 4
      instr#S args
   else
      instr#D args
   end if
}

macro MOVAP a, b   {defineinstr MOVAP, a, b}
macro ADDP a, b    {defineinstr ADDP, a, b}
macro MULP a, b    {defineinstr MULP, a, b}

宏 defineinstr 接受**可变数量的参数**。例如，它可以用于带有三个参数的 SHUFPS，用于带有两个参数的 ADDPS。宏发出一个带有 S 后缀的指令（如果 ssewidth = 4），或者否则是带有 D 后缀的相同指令。# 运算符连接名称；common 表示指令只重复一次。

你可以让这段代码更短。

macro defineinstr instr {
   macro instr [args] \{
   \common
      if ssewidth = 4
    instr#S args
      else
    instr#D args
      end if
   \}
}
defineinstr MOVAP
defineinstr ADDP
defineinstr MULP

这里有两个**嵌套宏**。外部宏 defineinstr 定义了内部宏 instr。请注意，instr 也是 defineinstr 的参数，所以当你写 defineinstr MOVAP 时，一个名为 MOVAP 的宏将被定义。这正是我想要的。花括号和 common 指令需要**用反斜杠转义**，因为 FASM 内部存在一些奇怪的解析问题。现在，新 SSE/SSE2 指令的定义非常简短易懂：只需写 defineinstr 和指令名称。参数的数量无关紧要。

我希望你现在相信 FASM 宏非常强大和优雅。如果不是，只需尝试创建一个带有可变数量参数的 C 预处理器宏，并将其与 FASM 等效项进行比较！

使用 GDI 和 DirectDraw 绘制位图

SetPixel 函数似乎是输出像素到屏幕的最慢方式。

case WM_PAINT: {
   PAINTSTRUCT ps;
   HDC hdc = BeginPaint(hWnd, &ps);
   for(x = 0; x < w; ++x)
      for(y = 0; y < h; ++y)
          SetPixel(hdc, x,y, b[i]);
   EndPaint(hWnd, &ps);
}

这个函数只适用于设置 1 或 2 个像素；对于大型图像来说太慢了。这是一个**我的程序**的简化版本，速度快得多。

int w, h; HBITMAP hBMP = 0, holdBMP; HDC hmem;
BITMAPINFOHEADER bi = {sizeof(BITMAPINFOHEADER), 0, 0, 1, 32, BI_RGB};
void ReDraw() {
   PaintJuliaSSE(b, ceil(w, 4), h);
   SetDIBits(hmem, hBMP, 0, h, b, 
            (BITMAPINFO*)&bi, DIB_RGB_COLORS);
            // Copy the bitmap to hmem
}
LRESULT CALLBACK WndProc(HWND hWnd,UINT msg,WPARAM wParam,LPARAM lParam) {
   switch(msg) {
      case WM_SIZE: {
         HDC hdc = GetDC(hWnd);
         if(hBMP) { DeleteObject(SelectObject(hmem, holdBMP)); DeleteDC(hmem); }
         w = LOWORD(lParam), bi.biWidth = ceil(w, 4),
         h = HIWORD(lParam), bi.biHeight = h;
         bi.biSizeImage = bi.biWidth * h * 4;
         b = (long*)realloc(b, bi.biSizeImage);
         hBMP = CreateCompatibleBitmap(hdc, ceil(w, 4), h);
         hmem = CreateCompatibleDC(hdc);
         holdBMP = (HBITMAP)SelectObject(hmem, hBMP);
         ReleaseDC(hWnd, hdc);
         ReDraw();
         return 0;  }
      case WM_PAINT: {
         PAINTSTRUCT ps; HDC hdc = BeginPaint(hWnd, &ps);
         BitBlt(hdc, 0, 0, w, h, hmem, 0, 0, SRCCOPY);  // Copy hmem to the screen
         DrawPolyline(hdc);
         EndPaint(hWnd, &ps);
         return 0;   }
      case WM_DESTROY:
         DeleteObject(SelectObject(hmem, holdBMP)); DeleteDC(hmem);
         PostQuitMessage(0); return 0;
      default:
         return DefWindowProc(hWnd,msg,wParam,lParam);
   }
}

每次图像完全更新时（例如，放大或缩小），你都会调用 ReDraw 函数。它使用 PaintJuliaSSE 在 b 中生成分形位图。然后调用 SetDIBits 函数。它将整个 b 数组复制到 hmem 内存 DC。

响应 WM_SIZE，程序设置新的位图大小并重新创建它。此外，应该重新分配位图的内存，并重新绘制整个分形。WM_PAINT 的处理程序只是将内存 DC 绘制（复制）到屏幕上，并添加一条折线来显示函数轨道。

分形图像存储在**内存 DC** 中，该 DC 在程序整个生命周期内都存在。有了这个图像，你可以通过绘制存储的图像并在其上方绘制新的折线位置来快速绘制移动的折线。这会导致折线**稍微闪烁**，但这是可以容忍的。添加第二个内存 DC 和第二个位图可以完全消除闪烁，但那样程序将变慢且内存消耗大。

SSE 指令并行处理 4 个值，所有上述 SSE 方法都需要每行像素数能被 4 整除。这就是为什么宽度（w）**向上取整**使用 ceil(w) 函数。位图会比需要的大一点，而在 WM_PAINT 处理程序中，你只从其中绘制 w * h 的矩形。

除了 GDI，该程序还支持 **DirectDraw**。我不会在这里展示代码，因为它很长而且很无聊（错误处理占用了大量空间）。本质上，它与 GDI 代码非常接近。创建一个屏幕外表面；将分形渲染到其中。在 WM_PAINT 中，将屏幕外表面绘制到主表面（屏幕）上，并绘制折线。现代显卡将两个表面都存储在视频内存中，因此 DirectDraw 的绘制速度可能比 GDI 绘制快 50 倍（我在我的笔记本电脑上使用简单的集成显卡芯片 Intel 915 GM 得到了这个结果）。与 GDI 相比，DirectDraw 将整体速度提高了 **5-6%**。

DirectDraw 支持 Windows 98 及更高版本。对于 Windows 95 用户，我编写了几行代码，以便在 DirectDraw 不可用时，程序可以优雅地降级到 GDI :).

结论

我的程序是一个非常快速的 Mandelbrot 和 Julia 集合生成器，使用 SSE/SSE2 和 DirectDraw。我不会声称它是“Win32 上最快的图形生成器”，就像 FFFF 的作者所说的那样，但你看到它相当快速且易于使用。有些方面仍需进一步完善，例如

• 存在更快的算法，它们可以**“猜测”**相邻像素的值；
• 你可以进一步优化 MOV-ADD 和 Vodnaya 方法，以在 **64 位处理器上**获得最大性能；
• 也许存在一种优化分支预测的技巧，但我没有找到；
• 程序无法保存当前位置和 C 参数，以便将来返回（下一版本将实现此功能）；
• 不支持更改调色板，但这不是我的目标。

所以代码是好的，但并不完美。请添加你的评论和改进建议。

你可以在许多不同的图形应用程序中使用本文提出的技术。**延迟和吞吐量分析**可以帮助你找到瓶颈。然后，你可以从关键路径中排除长延迟操作，在代码中创建几个独立的链，并通过重新排序指令来增加指令级并行性。最后，**汇编器宏**使你免于进行 SSE2 重复编写相同代码等枯燥的工作。

鸣谢

感谢 FFFF 的作者 Daniele Paccaloni。查看他的程序源代码帮助我发现了自己代码中的次优解决方案，总的来说，使我的程序更好。谢谢你，Daniele！

感谢我的老师 Tatyana Valenteenovna Korablina，她提供了关于分形和非线性动力学的有用信息。我也感谢帮助我测试的人：Anton Orlov、Larisa Alekseevna Morozova 和 Victor "vito"。特别感谢 Stas，他在那个星期五晚上留下来，让我使用了很多电脑进行测试。