|
SSE（为 Streaming SIMD Extensions 的缩写）是由 Intel 公司，在 1999 年推出 Pentium III 处理器时，同时推出的新指令集。如同其名称所表示的，SSE 是一种 SIMD 指令集。所谓的 SIMD 是指 single instruction, multiple data，也就是一个指令同时对多个数据进行相同的动作。较早的 MMX 和 AMD 的 3DNow! 也都是 SIMD 指令集。因此，SSE 本质上是非常类似一个向量处理器的。SSE 指令包括了四个主要的部份：单精度浮点数运算指令、整数运算指令（此为 MMX 之延伸，并和 MMX 使用同样的缓存器）、Cache 控制指令、和状态控制指令。这里主要是介绍浮点数运算指令和 Cache 控制指令。
目录：
[Part 1]
[Part 2]
[Part 3]
[Part 4]
[Part 5]
[Part 6]
[Part 7]
[Part 8]
文档内容：
[Part 1]
SSE 新增了八个新的 128 位缓存器，xmm0 ~ xmm7。这些 128 位的缓存器，可以用来存放四个 32 位的单精度浮点数。SSE 的浮点数运算指令就是使用这些缓存器。和之前的 MMX 或 3DNow! 不同，这些缓存器并不是原来己有的缓存器（MMX 和 3DNow! 均是使用 x87 浮点数缓存器），所以不需要像 MMX 或 3DNow! 一样，要使用 x87 指令之前，需要利用一个 EMMS 指令来清除缓存器的状态。因此，不像 MMX 或 3DNow! 指令，SSE 的浮点数运算指令，可以很自由地和 x87 指令，或是 MMX 指令共享。但是，这样做的主要缺点是，因为多任务操作系统在进行 context switch 时，需要储存所有缓存器的内容。而这些多出来的新缓存器，也是需要储存的。因此，既存的操作系统需要修改，在 context switch 时，储存这八个新缓存器的内容，才能正确支持 SSE 浮点运算指令。下图是 SSE 新增的缓存器的示意图：
目前 Intel 新推出的 x86 CPU 均已支持 SSE 指令集，包括 Pentium III、Pentium 4、Celeron 533A 及之后的处理器。AMD 方面，在其 Athlon 处理器上，新增的 Enhanced 3DNow! 指令集中，即包含了 SSE 中的整数运算指令和 Cache 控制指令，但仍不包含浮点数指令和状态控制指令。但在最新的 Athlon XP 处理器上，已将 SSE 中的浮点运算指令和状态控制指令加入。因此 Athlon XP 已经是完全支持 SSE 指令集了。
[Part 2]
SSE 的浮点数运算指令，大致上可以分成两种：packed 和 scalar。Packed 指令是一次对 XMM 缓存器中的四个浮点数（即 DATA0 ~ DATA3）均进行计算，而 scalar 则只对 XMM 缓存器中的 DATA0 进行计算。如下图所示：
SSE 指令和一般的 x86 指令很类似，基本上包括两种寻址方式：reg-reg 和 reg-mem。下面是两个例子：
addps xmm0, xmm1 ; reg-reg
addps xmm0, [ebx] ; reg-mem
指令的运算结果会覆盖到第一个参数中。例如，以上面的第一个例子来说，xmm0 缓存器会存放最后计算的结果。
另外，绝大部份需要存取内存的 SSE 指令，都要求地址是 16 的倍数（也就是对齐在 16 bytes 的边上）。如果不是的话，就会导致 exception。这是非常重要的。因为，一般的 32 位浮点数只会对齐在 4 bytes 或 8 bytes 的边上（根据 compiler 的设定而不同）。另外，若是处理数组中的数字，也需要特别注意这个问题。
[Part 3]
目前要在 C 或 C++ 程序中，利用 SSE 指令，有很多方法。最一般的方式，是利用内嵌式汇编语言。现在包括 Intel C++ Compiler 5.0 和 Microsoft Visual C++ 6.0 Processor Pack 的内嵌式汇编语言都支持 SSE 指令集。另外，它们也都支持 SSE 指令集的 intrinsics。Intrinsics 是一种外表类似一般的函数，但是实际上会被 compiler 直接编译成汇编语言的东西。以下是一些例子：
内嵌式汇编语言：
__asm addps xmm0, xmm1
__asm movaps [ebx], xmm0
...
__m128 data;
...
__asm
{
lea ebx, data
addps xmm0, xmm1
movaps [ebx], xmm0
}
使用 intrinsics：
__m128 data1, data2;
...
__m128 out = _mm_add_ps(data1, data2);
...
使用 intrinsics 可以增加程序的可读性，也比较容易使用。不过，在某些情形下，compiler 可能没办法产生最好的程序代码，而且，其产生的程序代码的效率，也会随着 compiler 的不同而有改变。但是，对于大部份的应用来说，使用 intrinsics 的好处通常是很明显的。因此，在这里我们也以 intrinsics 为主，而不讨论内嵌式汇编语言。
SSE 的 intrinsics 的名称看起来虽然怪异，但是它是有一定的规则的。基本上一个 SSE 的 intrinsics 的型式为：
_mm_<opcode>_<suffix>
其中，<opcode> 是指令的类别，像是 add、sub 等等。而 <suffix> 则是资料的种类。在 SSE 浮点运算指令中，只有两个种类：ps 和 ss。ps 是指 Packed Single-precision，也就是这个指令对缓存器中的四个单精度浮点数进行运算。而 ss 则是指 Scalar Single-precision，也就是这个指令只对缓存器中的 DATA0 进行运算。所以，像上面的 _mm_add_ps 指令，就是把两个四维向量相加的指令。
SSE 也需要一个新的数据型态，也就是上面所看到的 __m128。__m128 是一个 16 bytes（128 bits）的数据型态，对应 SSE 的 128 位缓存器。几乎所有的 SSE 浮点运算的 intrinsics 都是使用这个数据型态。比如说，_mm_add_ps 这个 intrinsic 的函数宣告为：
__m128 _mm_add_ps(__m128 a, __m128 b);
可以看到，它的参数和传回值的型态都是 __m128。基本上，这个 intrinsic 的动作就是把两个参数相加，并把结果以传回值的方式传回。
另外，要使用 SSE 的 intrinsics 之前，要记得先包含 xmmintrin.h 这个 header 档。如下：
#include <xmmintrin.h>
这样才能使用这些 instrinsics。另外，这个档案也包含了一些方便的宏定义，用来配合一些 SSE 的 intrinsics 使用。
[Part 4]
现在我们就先以一个很简单的例子，来说明 SSE 浮点运算指令的使用。下面是一个简单的程序片断：
- float input1[4] = { 1.2f, 3.5f, 1.7f, 2.8f };
float input2[4] = { -0.7f, 2.6f, 3.3f, -0.8f };
float output[4];
for(int i = 0; i < 4; i++) {
- output[i] = input1[i] + input2[i];
}
这个程序片断的动作非常简单，只是把 input1 和 input2 中的四组数字两两相加，并把结果写到 output 中。在执行完后，output 的内容应该是 { 0.5f, 6.1f, 5.0f, 2.0f }。
上面的程序很适合用 SSE 浮点运算指令来做。它所进行的四个加法运算，刚好可以用一个 SSE 浮点运算指令，也就是 addps 来完成。以下是修改后的程序：
- float input1[4] = { 1.2f, 3.5f, 1.7f, 2.8f };
float input2[4] = { -0.7f, 2.6f, 3.3f, -0.8f };
float output[4];
__m128 a = _mm_load_ps(input1);
__m128 b = _mm_load_ps(input2);
__m128 t = _mm_add_ps(a, b);
_mm_store_ps(output, t);
上面的程序中，先用两个 _mm_load_ps 把浮点数数据加载 __m128 的变数中。要注意的是，这里是假设这两个浮点数数组都是对齐在 16 bytes 的边上。如果不是的话，就不能用 _mm_load_ps 这个 intrinsic 来载入，而要改用 _mm_loadu_ps 这个 intrinsic。它是专门用来处理没有对齐在 16 bytes 边上的资料的。但是，它的速度会比较慢。
另外，因为 x86 的 little endian 特性，地址较低的 byte 会放在缓存器的右边。也就是说，若以上面的 input1 为例，在加载到 XMM 缓存器后，缓存器中的 DATA0 会是 1.2，而 DATA1 是 3.5，DATA2 是 1.7，DATA3 是 2.8。如果需要以相反的顺序加载的话，可以用 _mm_loadr_ps 这个 intrinsic。当然，在这个例子中，顺序并不影响结果，所以不需要利用这个 intrinsic。
一般来说，宣告一个 float 的数组，并不会对齐在 16 bytes 的边上。如果希望它会对齐在 16 bytes 的边上，以便利用 SSE 指令的话，Visual C++ 6.0 Processor Pack 和 Intel C++ compiler 都可以指定对齐方式。指定的方式如下：
- __declspec(align(16)) float input1[4];
在上例中，input1 会被对齐在 16 bytes 的边上。这样就可以直接用较快的 _mm_load_ps 来加载数据了。因为 SSE 浮点数指令常常需要数据对齐在 16 bytes 的边上，所以在 xmmintrin.h 也定义了一个宏 _MM_ALIGN16，是同样的意义。因此，上面的程序也可以写成：
- _MM_ALIGN16 float input1[4];
利用 SSE 浮点数指令会带来多大的差别呢？对 1,024 个浮点数进行测试的结果，结果如下（使用 Intel C++ Compiler 5.0 for Windows 编译，在 Intel Pentium III 1.0B Ghz 上执行）：
可以看到利用 SSE 浮点运算可以得到大约两倍于 x87 浮点运算的效率。事实上，因为 Pentium III 架构上的因素，所以，如果同时使用加法和乘法运算，最多可以得到四倍的效率。在后面会有一个例子。
[Part 5]
看了简单的例子之后，我们现在列出一些 SSE 浮点运算指令所支持的一些基本的运算：
| 指令 |
Intrinsic |
功能 |
addps/addss |
_mm_add_ps/_mm_add_ss |
加法 |
subps/subss |
_mm_sub_ps/_mm_sub_ss |
减法 |
mulps/mulss |
_mm_mul_ps/_mm_mul_ss |
乘法 |
divps/divss |
_mm_div_ps/_mm_div_ss |
除法 |
sqrtps/sqrtss |
_mm_sqrt_ps/_mm_sqrt_ss |
平方根 |
maxps/maxss |
_mm_max_ps/_mm_max_ss |
逐项取最大值 |
minps/minss |
_mm_min_ps/_mm_min_ss |
逐项取最小值 |
这些运算基本上都符合 IEEE 754 中的规范，根据其计算结果，设定适当的旗标（divide-by-zero、invalid 等等），或是产生 exception。这可以由一个 MXCSR 缓存器来设定。MXCSR 缓存器是一个 32 位的旗标缓存器，可以设定是否要产生各种 exception，并会记录上次的计算中，发生了哪些情况。下面是 MXCSR 缓存器的说明图：
以下是 MXCSR 各字段的说明：
| 字段 |
名称 |
Intrinsic 宏 |
说明 |
| IE |
Invalid Operation Flag |
_MM_EXCEPT_INVALID |
运算中发生 invalid operation |
| DE |
Denormal Flag |
_MM_EXCEPT_DENORMA |
试图加载 denormal value 或对 denormal value 进行运算 |
| ZE |
Divide-by-Zero Flag |
_MM_EXCEPT_DIV_ZERO |
运算中发生 divide by zero |
| OE |
Overflow Flag |
_MM_EXCEPT_OVERFLOW |
运算中发生 overflow |
| UE |
Underflow Flag |
_MM_EXCEPT_UNDERFLOW |
运算中发生 underflow |
| PE |
Precision Flag |
_MM_EXCEPT_INEXACT |
运算结果不精确（inexact） |
| DAZ |
Denormals Are Zeros |
N/A |
强制将 denormals 视为 0 |
| IM |
Invalid Operation Mask |
_MM_MASK_INVALID |
关闭 invalid operation exception |
| DM |
Denormal Operation Mask |
_MM_MASK_DENORM |
关闭 denormal operation exception |
| ZM |
Divide-by-Zero Mask |
_MM_MASK_DIV_ZERO |
关闭 divide by zero exception |
| OM |
Overflow Mask |
_MM_MASK_OVERFLOW |
关闭 overflow exception |
| UM |
Underflow Mask |
_MM_MASK_UNDERFLOW |
关闭 underflow exception |
| PM |
Precision Mask |
_MM_MASK_INEXACT |
关闭 inexact exception |
| RC |
Rounding Control |
N/A |
设定舍去方式 |
| FZ |
Flush to Zero |
N/A |
打开 Flush-to-zero 模式 |
上表中有几个地方是需要特别注意的：
- DAZ 若设为 1，则会开启 Denormals Are Zeros 模式。当此模式开启时，若来源参数中有任何 denormals，都会被视为 0，而且不会设定 DE 或产生 denormal operation exception。这个模式和 IEEE 754 规范不兼容，但是通常会得到更好的效率。另外 DAZ 目前只在 Intel Pentium 4 和 Intel Xeon 处理器上有支持。有关判断 DAZ 是否支持的方式，请参考 IA-32 Intel Architecture Software Developer's Manual, Volume 1: Basic Architecture, Section 11.6.3。
- RC 可以设定成四种方式：00 为 Round to nearest (even)，即将运算结果舍去至最接近数值。如果两个数值同样接近，则舍去至偶数。01 为 Round down，即往负无限大方向舍去。10 为 Round up，即往正无限大方向舍去。11 则是 Round toward zero (truncate)，即往零的方向舍去。
- FZ 若设为 1，则会开启 Flush-to-zero 模式。当此模式开启，且 underflow exception 被关闭时，若在运算时发生 underflow，则结果为 0（符号同真正结果），且会设定 PE 和 UE。若 underflow exception 未关闭，则此模式无任何影响。注意此模式和 IEEE 754 规范不兼容（规范要求 underflow 需产生 denormal 结果），但通常会得到更好的效率。
设定和读取 MXCSR 同样有 intrinsic 可以用，即 _mm_getcsr 和 _mm_setcsr。不过，因为 MXCSR 基本上是由 bit field 组成，因此在 xmmintrin.h 里面也定义了一些方便的宏。这些宏包括：
| 宏 |
功用 |
_MM_SET_EXCEPTION_STATE |
设定 exception 状态（IE ~ PE） |
_MM_GET_EXCEPTION_STATE |
取得 exception 状态（IE ~ PE） |
_MM_SET_EXCEPTION_MASK |
设定 exception mask（IM ~ PM） |
_MM_GET_EXCEPTION_MASK |
取得 exception mask（IM ~ PM） |
_MM_SET_ROUNDING_MODE |
设定舍去方式 |
_MM_GET_ROUNDING_MODE |
取得目前的舍去方式 |
_MM_SET_FLUSH_ZERO_MODE |
设定 Flush-to-zero 模式（_MM_FLUSH_ZERO_ON 和 _MM_FLUSH_ZERO_OFF） |
_MM_GET_FLUSH_ZERO_MODE |
取得 Flush-to-zero 模式设定 |
舍去方式的设定包括：
_MM_ROUND_NEAREST |
Round to nearest (even) |
_MM_ROUND_DOWN |
Round down |
_MM_ROUND_UP |
Round up |
_MM_ROUND_TOWARD_ZERO |
Round toward zero (truncate) |
在开机时，所有的 exception mask 都是设定为 1，也就是所有的 exception 都是关闭的。如果有需要的话，可以将其开启。下面的程序片断是一个例子：
- _MM_ALIGN16 float test1[4] = { 0, 0, 0, 1 };
_MM_ALIGN16 float test2[4] = { 1, 2, 3, 0 };
_MM_ALIGN16 float out[4];
_MM_SET_EXCEPTION_MASK(0); // 打开所有的 exception
__try {
- __m128 a = _mm_load_ps(test1);
- __m128 b = _mm_load_ps(test2);
- a = _mm_div_ps(a, b);
- _mm_store_ps(out, a);
}
__except(EXCEPTION_EXECUTE_HANDLER) { // exception handler
- if(_mm_getcsr() & _MM_EXCEPT_DIV_ZERO) {
- cout << "Divide by zero" << endl;
- }
-
- return;
}
cout << out[3] << endl;
上面的程序在执行时，会产生 divide by zero 的 exception，并进入 exception handler。在 exception handler 中，因为 ZE 会被设定，因此会显示出 "Divide by zero" 的错误讯息。如果把 _MM_SET_EXCEPTION_MASK(0); 这一行去掉，则不会产生 exception，而会直接输出正无限大（1.#INF）的结果。
[Page 6]
现在我们来看看一个实际的例子，是一个 3D 绘图运算中很常需要进行的动作，也就是三维向量的内积。
两个三维向量 <x1, y1, z1> 和 <x2, y2, z2> 的定义为：
x1x2 + y1y2 + z1z2
但是 SSE 浮点运算指令是以四维向量为单位的，而且也没有计算内积的指令。要直接对一个三维向量计算出内积，会需要相当多的数据重排动作，反而会浪费很多时间。
而且，如果三维向量在内存中是以 x1, y1, z1, x2, y2, z2, ... 的方式排列的话，还会有对齐的问题。因为一个三维向量只需要 12 bytes，所以无法让每个向量都对齐在 16 bytes 的边上。这也会降低效率。
一个简单的方法是改变数据的排列方式。如果我们把数据改成以下面的方式排列：
| x1, x2, x3, x4, ... |
| y1, y2, y3, y4, ... |
| z1, z2, z3, z4, ... |
那我们就同时解决了两个问题。首先，因为一次读入四个 x 值（及 y 值和 z 值），所以不会有对齐的问题。另外，运算也变得容易。下面的程序片断一次可以算出四个向量内积：
- __m128 x1 = _mm_load_ps(vec1_x);
- __m128 y1 = _mm_load_ps(vec1_y);
- __m128 z1 = _mm_load_ps(vec1_z);
-
- __m128 x2 = _mm_load_ps(vec2_x);
- __m128 y2 = _mm_load_ps(vec2_y);
- __m128 z2 = _mm_load_ps(vec2_z);
-
- __m128 t1 = _mm_mul_ps(x1, x2);
- __m128 t2 = _mm_mul_ps(y1, y2);
- t1 = _mm_add_ps(t1, t2);
- t2 = _mm_mul_ps(z1, z2);
- t1 = _mm_add_ps(t1, t2);
-
- _mm_store_ps(output, t1);
由于 SIMD 指令的特性，很多时候，不要把工作当成向量来考虑。而应该是把多个工作，结合到一个向量中来进行。以上面的例子来说，一个三维向量的内积，是很难有效率地用 SSE 指令来完成的。事实上，配合 SSE 的数据重排指令（后面会介绍），一个四维向量的内积，需要八个 SSE 浮点运算指令才能完成。这并不是很有效率。但是，如果我们不要以向量为单位，而改成考虑一次对四个向量进行运算的话（如同上面的例子），四个三维向量的内积，由五个指令就可以完成了。四个四维向量的内积也只需要七个指令。
以下是以 1,024 个三维向量，以上述的数据排列方式，在 Intel Pentium III 1.0B Ghz 上测试的结果（使用 Intel C++ Compiler 5.0 for Windows 编译）：
可以看到使用 SSE 浮点运算指令，得到的效率比使用 x87 浮点数要高出甚多，几乎达到三倍。
[Page 7]
前面的向量内积的测试，是以 1,024 个三维向量，对一个固定的三维向量进行内积得到的结果。它的数据量算是相当小的，甚至可以整个放到 L1 cache 里面。但是，实际情形常常不是这样。我们常常需要对一大堆数据进行运算，因此数据是无法放在 cache 里面的，而需要从主存储器中取得。这样一来，速度会变得相当的慢。下面是对 102,400 个三维向量进行同样动作的结果：
可以看到效率差了很多，这是因为数据都需要从主存储器中读取，因而读取的时间成为瓶颈，而不是计算的时间。在这种情形下，SSE 浮点运算指令速度虽快，也没有什么帮助。而且，因为内积的计算是相当简单的，所以这也让 SSE 和 x87 的结果差距变得很小。
不过，现在我们会讨论到 SSE 的 "streaming" 部份。SSE 除了运算的指令之外，还支持了一些 cache 控制指令。我们在这里先介绍两个较简单的指令：prefetch 和 movntps。prefetch 指令实际上有四个不同的指令，包括 prefetch0、prefetch1、prefetch2、和 prefetchnta。不过，它们都是用同一个 intrinsic 表示的，也就是 _mm_prefetch。
prefetch 指令的主要目的，是提前让 CPU 加载稍后运算所需要的数据。通常是在对目前的数据进行运算之前，告诉 CPU 加载下一笔数据。这样就可以让目前的运算，和加载下一笔资料的动作，可以同时进行。如果运算的复杂度够高的话，这样可以完全消除读取主存储器的 latency。不同的 prefetch 指令则是告诉 CPU 将数据加载不同层次的 cache。不过，最常用的还是 prefetchnta，这个指令会把数据加载到离 CPU 最近的 cache 中（通常是 L1 cache 或 L2 cache），适用于数据在短时间内就会用到的情形。
另外 prefetch 指令不会产生任何 exception。它本质上只是一个 hint，CPU 并不一定会真的进行加载的动作。所以，即使 prefetch 一个不合法的内存地址，也不会产生错误。这让程序可以不用处理讨厌的边界问题。
除了 prefetch 之外，另一个指令是 movntps，它的 intrinsics 是 _mm_stream_ps。这个指令的用途，是要求 CPU 在写入数据的时候，不要把数据写到 cache 中，而是直接将数据写到主存储器中。实际上它以 write combining 的方式写入的。为什么要这样做呢？这是因为，很多时候计算的结果并不是立刻需要用到的，通常是很久以后才会用到。所以，这些数据如果被放在 cache 中，完全是浪费空间。而且，更糟的是，它们可能会把 cache 中有用的数据挤掉，而使得这些数据常常需要重新从主存储器中加载。因此，如果让这些资料不要被放在 cache 中，就可以避免这种问题。
对上面的程序，加上适当的 prefetch，并利用 movntps 指令，可以修改成类似下面的程序片断：
- __m128 x1 = _mm_load_ps(vec1_x);
- __m128 y1 = _mm_load_ps(vec1_y);
- __m128 z1 = _mm_load_ps(vec1_z);
-
- __m128 x2 = _mm_load_ps(vec2_x);
- __m128 y2 = _mm_load_ps(vec2_y);
- __m128 z2 = _mm_load_ps(vec2_z);
-
- _mm_prefetch((const char*)(vec1_x + next), _MM_HINT_NTA);
- _mm_prefetch((const char*)(vec1_y + next), _MM_HINT_NTA);
- _mm_prefetch((const char*)(vec1_z + next), _MM_HINT_NTA);
-
- _mm_prefetch((const char*)(vec2_x + next), _MM_HINT_NTA);
- _mm_prefetch((const char*)(vec2_y + next), _MM_HINT_NTA);
- _mm_prefetch((const char*)(vec2_z + next), _MM_HINT_NTA);
-
- __m128 t1 = _mm_mul_ps(x1, x2);
- __m128 t2 = _mm_mul_ps(y1, y2);
- t1 = _mm_add_ps(t1, t2);
- t2 = _mm_mul_ps(z1, z2);
- t1 = _mm_add_ps(t1, t2);
-
- _mm_stream_ps(output, t1);
同样对 102,400 个三维向量进行测试，结果为：
由于 x87 无法利用 movntps 指令，所以效率较差。可以看到，加上 prefetch 和 movntps 指令后，SSE 浮点运算的效率提高了约 50%。事实上，它所使用的内存频宽达到约 500MB/s，已经是理论频宽（1,066MB/s）的一半左右了。
[Page 8]
前面介绍的计算都是蛮简单的计算。现在我们来看一个比较复杂的例子。假设现在有一堆数字，都是单精度的浮点数（float），而且范围都在正负 0.5 之间。然后现在要计算它们的自然指数（exp）。当然我们知道 SSE 并没有这个指令。事实上除了简单的四则运算（加、减、乘、除）之外，SSE 和计算有关的指令，只有平方根而已。其它复杂的东西像是各种超越函数都是没有的。
不过，这并不表示要算超越函数就不能用 SSE 了。这只是表示我们得自己算而已。以自然指数来说，最简单的算法当然是用 Taylor series。事实上有更快的算法，不过我们为了举例方便，所以就直接用 Taylor series 了。因为我们的数字范围是在正负 0.5 之间，所以只需要取到第九项，也就是 x 的八次方项，就已经很接近单精度浮点数所能表示的精确度了。下面是利用 Taylor series 算 exp(x) 的程序片断：
- int i;
- float result = coeff[8] * x;
-
- for(i = 7; i >= 2; i--) {
- result += coeff[i];
- result *= x;
- }
-
- return (result + 1) * x + 1;
上面的程序片断中，coeff 数组需要事先设定成 1/n! 的值。这可以事先算好，直接放在程序里面，所以不需要花时间。这个程序片断是利用 Horner 法来计算 Taylor series，以减少计算量。在上例中，每个 X 需要 8 个乘法和 8 个加法，共 16 个计算。
这个程序要改成用 SSE 来写，并不会很困难。不过，coeff 数组需要修改成 __m128 的数组，而且每项存放四个相同的数字（也就是 1/n!）。这样可以改成下面的程序：
- int i;
- __m128 X = _mm_load_ps(data);
- __m128 result = _mm_mul_ps(coeff_sse[8], X);
-
- for(i = 7; i >=2; i--) {
- result = _mm_add_ps(result, coeff_sse[i]);
- result = _mm_mul_ps(result, X);
- }
-
- result = _mm_add_ps(result, sse_one);
- result = _mm_mul_ps(result, X);
- result = _mm_add_ps(result, sse_one);
-
- _mm_store_ps(out, result);
上面的程序片断和 x87 的版本长得几乎一样，当然它是一次对四个浮点数做运算，所以应该会快四倍。下图是测试结果：
这些测试结果是对 1,024 个随机产生的数字进行的。第一个结果是直接用 exp 函式，它里面是用 x87 的指令像 F2XM1 和 FYL2X 等指令来完成计算的。当然这个函式并没有范围的限制，而且它得到的结果也是最精确的。不过，它的速度非常的慢，平均每个数字需要 128 个 cycles 才能算完。第二个结果是用 x87 版的 Taylor series 来计算，也就是第一个程序片断。它很明显要快得多，平均每个数字只需要 66 个 cycles，几乎是快了一倍。不过，其实这里就产生一个问题：为什么 16 个计算会需要 66 个 cycles 才能完成？
第三个结果是用 SSE 计算 Taylor series，也就是第二个程序片断。它比 x87 版本快了近四倍，平均每个数字只需要 17 cycles。可是，因为 SSE 指令一次对四个数字进行运算，所以其实每四个数字共花了 68 cycles。这里也是同样的问题：为什么 16 个计算会需要 68 cycles？这表示平均一个 SSE 乘法指令或加法指令，花了四个 cycles 执行。但是，理论上 Pentium III 不是可以每 cycle 执行一个 SSE 指令吗？
也许问题是出在上面的程序片断里面的循环（for(i = 7; i >=2; i--)）。所以，也许可以把它展开试试。它会得到上面的第四个结果，实际上也是快了不少，平均每个四个数字花了约 43 cycles。也就是说，平均每个运算花了 2.7 cycles。这样已经算是相当不错了。可是，它和「每 cycle 执行一个 SSE 指令」还有一段距离。
事实上，这是因为这些指令都有相依性。在上面的程序片断中，result 这个变量不但从头用到尾，而且每个指令都用到这个变量，也会修改这个变量。因此，第二个指令要等到第一个指令执行完后，才能执行。在实际执行时，因为 Pentium III 的 instruction reorder buffer 够大，所以可以达到一部份程度的平行执行，但是效率还不是非常理想。
要尽可能达到平行执行的效率，就要用另一个方法来展开循环。也就是说，不展开内部的循环，而是一次对多组数据进行计算。因为多组数据的计算之间并没有相依性，要平行执行就容易多了。上图中的第五个结果，就是一次对四组数字进行计算（也就是共 16 个数字），而内部的循环并没有展开。这样平均每四个数字只需要 34 cycles，也就是平均每个运算花了约 2.1 cycles。比起第四个结果，又要再快了一些。
那如果一次对八组数字进行计算，会再变快吗？答案是不会，它反而会变慢。为什么呢？因为 SSE 只有八个浮点数缓存器，所以如果一次对八组数字进行运算，那一定会有些东西无法放到缓存器中，而需要放到内存里面。虽然这些数据几乎是一定会放在快速的 L1 cache 中，但是它还是会需要一些额外的 load/store 指令。所以反而不会增加效率。不过，我们还是可以把内部的循环也展开，因为这个循环并不大，所以展开并不会太夸张。把内部的循环也展开后，得到的就是第六个结果，也就是平均每四个数字只需要 22.6 cycles，平均每个运算只需要 1.4 cycles！这已经离「每 cycle 执行一个 SSE 指令」不远了！
还有办法再加快吗？恐怕不容易了。因为它已经达到约 2,837MFLOPS（在 Pentium III 1.0B Ghz 上执行）。事实上，这个结果比前面计算向量内积的速度还要快。所以，要再加快应该是很难了。
所以，这里要讨论的重点就是，当计算相当冗长，而且相依性很高的时候，光是展开内部的循环是不够的，应该要考虑也展开外部的循环（也就是一次对多组数据进行同样的计算）。这样通常比展开内部的循环要更有效。当然，如果情况许可的话，也可以连内部的循环也展开，通常可以达到更高的效率。
另外，可能会有人问：这样计算自然指数，精确度有多高呢？和用 exp 函数比较又是如何？在这个例子中，1,024 个随机数字中，误差都小于 1.2E-7，也就是单精度浮点数的 epsilon。这表示有误差的情形，都只有最后一个位不同。因此，这应该已经是非常高的精确度了。
| 
中搜索 