优化汇编例程（13）

12.8. Core2上相同的例子

归功于更好的μop融合与更强大的执行单元，相同的循环在Core2运行得更高效。例子12.6c与12.6d都是Core2上运行的好候选。在例子12.6c中的j1 L1指令应该改为jb L1，以启用32位模式中的宏操作融合。这个改变是可能的，因为在例子12.6c中，循环计数器不能是负数。

现在，我们将分析对一个Core2处理器，例子12.6d中的循环的资源使用。循环中每条指令使用的资源列在下面的表12.2中。

指令		指令长度		μop融合域		每个执行端口的μop
						端口0	端口1	端口5	端口2	端口3	端口4
movapd xmm1, [esi+eax]		5		1					1
mulpd xmm1, xmm2		4		1		1
addpd xmm1, [edi+eax]			5	1			1		1
movapd [edi+eax], xmm1			5	1						1	1
add eax, 16		3		1		x	x	x
js L1		2		1				1
Total	24				6	1.33	1.33	1.33	2	1	1

表12.2. Core2上DAXPY循环的μop总数

这里预解码不是一个问题，因为循环的大小小于64字节。没有必要对齐循环，因为其大小小于50字节。如果循环超过64字节，那么我们将注意到预解码器在一个时钟周期里仅可以处理前3条指令，因为它不能处理超过16字节。对解码器没有限制，因为所有的指令每条仅产生一个μop。解码器每时钟周期可以处理4条指令，因此解码时间将是每迭代6/5 = 1.5时钟周期。

端口0、1与2每个有一个不去别处的μop。另外，指令ADD EAX, 16可去往这三个端口的其中一个。如果这些指令在端口0、1与2间均匀分布，那么这些端口平均每迭代忙碌1.33时钟周期。

端口3每迭代接收2个μop，用于读X[i]及Y[i]。因此，我们可以断定端口3上的内存读是DAXPY循环的瓶颈，执行时间是每迭代2时钟周期。

例子12.6c比例子12.6d多一条指令，我们已经分析过它了。这个额外的指令是CMP EAX, ECX，它去往端口0、1与2的其中一个。这将增加这些端口的压力到1.67，这仍然小于端口3上的2个μop。在32位模式中，额外的μop可以通过宏操作融合消除。浮点加法与乘法单元都有吞吐率1。因此，这些单元不是DAXPY循环中的瓶颈。

我们需要考虑循环中寄存器读暂停的可能性。例子12.6c循环里有4个读，但不修改的寄存器。它们是寄存器ESI，EDI，XMM2与ECX。这会导致一次寄存器读暂停，如果在4个连续μop内读所有这些寄存器。不过读这些寄存器的μop被分得如此开，在4个类型μop中，不可能读超过其中3个寄存器。例子12.6d少读一个寄存器，因为ECX不再使用。因此，在这些循环里，寄存器读暂停不太可能造成麻烦。

结论是，Core2处理器运行DAXPY循环的时间只需PM的一半。在PM上使优化特别困难的μop冲突问题，在Core2设计中被消除了。

12.9. ​​​​​​​Sandy Bridge上相同的例子

向量寄存器中元素数量，在具有AVX指令集的处理器上，使用YMM寄存器加倍了， 如下例所示：

; Example 12.6f. Loop of DAXPY for processors with AVX

.data

SignBit DD 0, 80000000H                                           ; qword with sign bit set

n = 100                                                                           ; Define constant n (divisible by 4)

 

.code

      vbroadcastsd ymm0, [SignBit]                             ; preload sign bit x 4

      mov eax, -n * 8                                                       ; Index = -n * sizeof(double)

      lea rsi, [X + 8 * n]                                                    ; Point to end of array X (aligned)

      lea rdi, [Y + 8 * n]                                                    ; Point to end of array Y (aligned)

      vbroadcastsd ymm2, [DA]                                     ; Load DA x 4

      vxorpd ymm2, ymm2, ymm0                                ; Change sign

 

      align 32

L1: vmulpd ymm1, ymm2, [rsi+rax]                           ; X[i]*(-DA)

      vaddpd ymm1, ymm1, [rdi+rax]                           ; Y[i]-X[i]*DA

      vmovapd [rdi+rax], ymm1                                     ; Store (vmovupd if not aligned by 32)

      add rax, 32                                                                ; Add size of four elements to index

      jl L1                                                                             ; Loop back

      vzeroupper                                                                ; If subsequent code is non-VEX

现在我们将分析，对Sandy Bridge处理器，例子12.6f中循环的资源使用。循环中的每条指令的资源列出在下面表12.3中。

指令		指令长度		μop融合域		每个执行端口的μop
						端口0	端口1	端口5	端口2	端口3	端口4
movapd xmm1, [esi+eax]		5		1		1			1+
mulpd xmm1, xmm2		5		1			1			1+
addpd xmm1, [edi+eax]			5	1					1		1+
movapd [edi+eax], xmm1			4	½				½
add eax, 16		2		½				½
Total	21				4	1	1	1	2	1+	1+

表12.3. Sandy Bridge上DAXPY循环的μop总数（例子12.6f）

指令长度不重要，如果我们假定循环在μop缓存中。ADD与JL指令很可能融合在一起，一起仅使用一个μop。有两个256位读操作，每个连续两个时钟周期使用一个读端口，在表中表示为1+。使用两个读端口（端口2与3），在两个时钟周期里，将有两次256位读的吞吐率。其中一个读端口，对在第二个时钟周期里的写，进行地址计算。写端口（端口4）被256位写占据两个时钟周期。限制因素将是读与写操作，将两个读端口以及写端口用到最大限度。端口0、1与5用了一半的容量。没有循环携带依赖链。因此，预期吞吐率是两个时钟周期，一次迭代的四个计算。

不幸的是，由于大量的内存操作，缓存库冲突的可能性很高。很难预测哪些操作将同时读或写，因为同时进行许多迭代。我在Sandy Bridge上的测试显示，每迭代2时钟周期的理论时间，从来没有得到过。最好的结果是每迭代2.2时钟周期，当RSI与RDI间的距离是4k字节倍数时。在其他情形下，时间上每得到2.6时钟周期或更多。

12.10. ​​​​​​​使用FMA4相同的例子

AMD Bulldozer是第一个提供融合乘加指令的x86处理器。这些指令可以显著提高吞吐率。

; Example 12.6g. Loop of DAXPY for processors with FMA4

.code

      mov eax, -n * 8                                                       ; Index = -n * sizeof(double)

      lea rsi, [X + 8 * n]                                                    ; Point to end of array X (aligned)

      lea rdi, [Y + 8 * n]                                                    ; Point to end of array Y (aligned)

      vmovddup xmm2, [DA]                                          ; Load DA x 2

 

      align 32

L1: vmovapd xmm1, [rsi+rax]                                      ; X[i]

      vfnmaddpd xmm1, xmm1, xmm2, [rdi+rax]       ; Y[i]-X[i]*DA

      vmovapd [rdi+rax], xmm1                                      ; Store

      add rax, 16                                                                 ; Add size of two elements to index

      jl L1                                                                              ; Loop back

这个循环不受执行单元的限制，而是内存访问。不幸的是，在Bulldozer上缓存库冲突非常频繁。这将消耗时间从理论上的每迭代2时钟周期，增加到几乎4。在Bulldozer上，使用更大的YMM寄存器没有优势。另一方面，YMM双倍加长指令的解码有些低效。

12.11. ​​​​​​​使用FMA3相同的例子

Intel Haswell处理器提供3操作数形式、称为FMA3的融合乘加指令。AMD Piledriver同时支持FMA3与FMA4。

; Example 12.6h. Loop of DAXPY with FMA3, using xmm registers

.code

      mov eax, -n * 8                                                             ; Index = -n * sizeof(double)

      lea rsi, [X + 8 * n]                                                          ; Point to end of array X (aligned)

      lea rdi, [Y + 8 * n]                                                          ; Point to end of array Y (aligned)

      vmovddup xmm2, [DA]                                                ; Load DA x 2

      align 32

 

L1: vmovapd xmm1, [rdi+rax]                                           ; Y[i]

      vfnmadd231pd xmm1, xmm2, [rsi+rax]                    ; Y[i]-X[i]*DA

      vmovapd [rdi+rax], xmm1                                            ; Store

      add rax, 16                                                                       ; Add size of two elements to index

      jl L1                                                                                    ; Loop back

在Intel处理器上，这个循环每迭代需要1 ~ 2时钟周期。使用YMM寄存器，吞吐率几乎加倍：

; Example 12.6i. Loop of DAXPY with FMA3, using ymm registers

.code

      mov eax, -n * 8                                                                 ; Index = -n * sizeof(double)

      lea rsi, [X + 8 * n]                                                              ; Point to end of array X

      lea rdi, [Y + 8 * n]                                                              ; Point to end of array Y

      vbroadcastsd xmm2, [DA]                                               ; Load DA x 4

      align 32

 

L1: vmovupd ymm1, [rsi+rax]                                               ; X[i]

      vfnmadd231pd ymm1, ymm2, [rdi+rax]                       ; Y[i]-X[i]*DA

      vmovupd [rdi+rax], ymm1                                               ; Store

      add rax, 32                                                                          ; Add size of two elements to index

      jl L1                                                                                       ; Loop back

例子12.6i无需假定数组对齐到32。因此，VMOVAPD可以改为VMOVUPD。在Haswell与Broadwell上，这个循环的执行时间测得为2时钟周期，在Skylake上是1.5时钟周期。瓶颈不是这里执行的单元，而是缓存效果以及可能分支指令。

12.12. ​​​​​​​循环展开

一个重复n次的循环可以被一个重复n / r次、每次执行r次计算的循环所替代，其中r是展开因子。n最好能被r整除。

循环展开可用于以下目的：

• 减少循环开销。每次计算的循环开销除展开因子r。仅在循环开销对计算时间有显著贡献时，这才有用。如果其他某个瓶颈限制了执行速度，没有理由展开一个循环。例如，上面例子12.6e中的循环就不能从展开获益。
• 向量化。为了使用有r元素的向量寄存器，循环必须展开r或r的倍数。例子12.6e中的循环展开2次，以使用有两个双精度值的向量。如果我们使用了单精度数，那么可以展开循环4次，使用有4个元素的向量。
• 改进分支预测。通过展开循环，使重复计数n / r不超过特定CPU可以预测的最大重复计数，可以改进循环退出分支的预测。
• 改进缓存。如果循环遭遇许多数据缓存不命中或缓存竞争，以对特定处理器最优的方式安排内存读写，可能是有利的。在现代处理器上，这很少需要。细节参考微处理器提供商提供的优化手册。
• 消除整数除法。如果循环包含一个表达式，其中循环计数器i除一个整数r，或者计算i模r，通过展开循环r次，可以避免这个整数除法。
• 消除循环内分支。如果带有周期r重复模式的循环内有一个分支或switch语句，可以通过展开该循环r次消除之。例如，如果一个if-else分支每二次去往同一个方向，这个分支可以通过展开2次消除。
• 打破循环携带依赖链。在某些情形里，可以通过使用多个累加器，打破循环携带依赖链。展开因子r等于累加器数。参考第59页，例子9.3b。
• 减少对归纳变量的依赖。如果从前面的迭代值计算一个归纳变量的时延很高、成为一个瓶颈，那么有可能可以通过展开r次，从序列中后r个位置处的值计算归纳变量的每个值。
• 完全展开。在r = n时，完全展开一个循环。这完全消除了循环开销。每个是循环计数器函数的表达式可被常量替代。可以消除仅依赖循环计数器的每个分支。参考第111页例子。

在重复计数n不能被展开因子r整除时，循环指令存在一个问题。将存在n modulo r的余数个、不能在循环里完成的额外计算。这些额外计算必须在主循环前或后完成。

在重复计数不能被展开因子整除时，正确获取额外计算，可能相当困难；如果使用如例子12.6d与e中的负索引，它变得特别棘手。下面的例子再次展示DAXPY算法，这次使用单精度并展开4次。在这个例子中，n是不一定被4整除的变量。数组x与y必须对齐到16。（为了清晰起见，忽略特定于PM处理器的优化）。

; Example 12.7. Unrolled Loop of DAXPY, single precision.

.data

align 16

SignBitS DD 80000000H                                     ; dword with sign bit set

 

.code

      mov eax, n                                                      ; Number of calculations, n

      sub eax, 4                                                        ; n - 4

      lea esi, [X + eax*4]                                         ; Point to X[n-4]

      lea edi, [Y + eax*4]                                        ; Point to Y[n-4]

      movss xmm2, DA                                           ; Load DA

      xorps xmm2, SignBitS                                   ; Change sign

      shufps xmm2, xmm2, 0                                ; Get -DA into all four dwords of xmm2

      neg eax                                                            ; -(n-4)

      jg L2                                                                  ; Skip main loop if n < 4

 

L1:                                                                           ; Main loop rolled out by 4

      movaps xmm1, [esi+eax*4]                         ; Load 4 values from X

      mulps xmm1, xmm2                                     ; Multiply with -DA

      addps xmm1, [edi+eax*4]                           ; Add 4 values from Y

      movaps [edi+eax*4],xmm1                         ; Store 4 results in Y

      add eax, 4                                                       ; i += 4

      jle L1                                                                ; Loop as long as <= 0

 

L2:                                                                          ; Check for remaining calculations

      sub eax, 4                                                       ; = -remainder

      jns L4                                                               ; Skip extra loop if remainder = 0

 

L3:                                                                          ; Extra loop for up to 3 remaining calculations

      movss xmm1, [esi+eax*4+16]                    ; Load 1 value from X

      mulss xmm1, xmm2                                     ; Multiply with -DA

      addss xmm1, [edi+eax*4+16]                    ; Add 1 value from Y

      movss [edi+eax*4+16],xmm1                    ; Store 1 result in Y

      add eax, 1                                                      ; i += 1

      js L3                                                                ; Loop as long as negative

L4:

展开在数组上进行计算的循环的另一个解决方案是，将数组扩展到具有r-1个不使用的空间，将重复计数n取整到展开因子r最接近的倍数。这消除了计算余数（n mod r）以及余下计算的额外循环的需要。不使用的数组元素必须初始化为0或其他有效的浮点值，以避免次正规数、NAN、溢出、下溢，或其他任何会减慢浮点计算的情形。如果数组是整数类型，只需避免除0。

仅在有原因这样做，且可以获得速度的显著提升时，才进行循环展开。应该避免过分的循环展开。循环展开的坏处有：

• 代码变得更大，在代码缓存中占据更多空间。这会导致代码缓存不命中，使展开得不偿失。注意，在孤立测试循环时，不能检测到代码缓存不命中。
• 许多处理器有一个循环缓冲来提升非常小循环的速度。循环缓冲被限制为20 ~ 40条指令或64字节代码，依赖于处理器。如果超出循环缓冲的大小，循环展开很可能降低性能，（处理器特定细节在手册3《Intel，AMD与VIA CPU微架构》中提供）。
• 某些处理器有一个有限大小的μop缓存。这个μop缓存非常有价值，应该经济地使用之。展开的循环在μop缓存中占据更多的空间。注意，在孤立测试循环时，不能检测到μop缓存不命中。
• 在n不能被r整除时，在展开循环外部需要进行额外计算，使得代码更复杂、笨拙且增加了分支的数量。
• 展开的循环可能需要更多的寄存器，如用于多个累加器。

12.13. 使用掩码寄存器的向量循环（AVX512）

AVX512指令集提供了受掩蔽的操作。在循环计数不能被向量大小整除时，这个特性对掩蔽超出的向量元素是有用的。

; Example 12.8. Vectorized DAXPY loop using mask registers

.data

align 64

countdown dd 15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0

sixteen dd 16

 

.code

      lea rsi, [X]                                                                                 ; point to beginning of X

      lea rdi, [Y]                                                                                 ; point to beginning of Y

      mov edx, n                                                                               ; number of elements, n

      vmovd xmm0, edx

      vpbroadcastd zmm0, xmm0                                                 ; broadcast n

      vpaddd zmm0, zmm0, zmmword [countdown]                ; counts = n+countdown

      vpbroadcastd zmm1, dword [sixteen]                                ; broadcast 16

      vbroadcastss zmm2, dword [DA]                                         ; broadcast DA

      xor eax, eax                                                                              ; loop counter i = 0

L1:                                                                                                    ; Main loop rolled out by 16

      vpcmpnltud k1, zmm0, zmm1                                              ; mask k1 = (counts >= 16)

      vpsubd zmm0, zmm0, zmm1                                                ; counts -= 16

      vmovups zmm3 {k1}{z}, zmmword [rdi+rax*4]                 ; load Y[i]

      vfnmadd231ps zmm3 {k1}, zmm2, zword [rsi+rax*4]      ; Y[i]-DA*X[I]

      vmovups zmmword [rdi+rax*4] {k1}, zmm3                      ; store Y[i]

      add rax, 16                                                                                ; i += 16

      cmp eax, rdx                                                                             ; while i < n

      jb L1

这里，在k1寄存器中的掩码，对所有有效的元素置1比特，对超出n的元素是0。保存指令被k1掩蔽，避免保存超出n个Y元素。使用相同的掩码，掩蔽读与计算指令同样也是好的。这里，必须使用0掩蔽读指令，以避免对向量寄存器先前值的一个假依赖。掩蔽计算节省了能耗，避免可能的次正规值等导致的异常及性能损失。向一条512位向量指令添加掩码，没有代价。

产生掩码的指令不需要与计算相同的向量大小与粒度。例如，如果计算有双精度（64位粒度），我们可以从具有32位粒度的256位寄存器，甚至16位粒度的128位寄存器，产生所需的8位掩码，如果n保证小于216且支持AVX512BW。

要性能最好，数组X与Y最好对齐到64。

如果向量加法的执行单元是一个瓶颈，制作掩码可以使用整数指令而不是向量指令。这在下面的例子中展示：

; Example 12.9. Vectorized DAXPY, mask calculated with integer instr.

; Note: This version works only for n < 256!

 

.code

      lea rsi, [X]                                                                               ; point to beginning of X

      lea rdi, [Y]                                                                               ; point to beginning of Y

      vbroadcastss zmm2, dword[DA]                                        ; broadcast DA

      mov edx, n                                                                             ; number of elements, n

      mov ecx, -1                                                                             ; fill with all 1's

      xor eax, eax                                                                            ; loop counter i = 0

L1:                                                                                                  ; Main loop rolled out by 16

      bzhi ecx, ecx, edx                                                                  ; zero bit positions > edx

      kmovw k1, ecx                                                                       ; copy 16 bits to mask register

      vmovups zmm3 {k1}{z}, zmmword [rdi+rax*4]               ; load Y[i]

      vfnmadd231ps zmm3 {k1}, zmm2, zword [rsi+rax*4]    ; Y[i]-DA*X[I]

      vmovups zmmword [rdi+rax*4] {k1}, zmm3                    ; store Y[i]

      add rax, 16                                                                              ; i += 16

      sub edx, 16                                                                             ; count down n

      ja L1                                                                                         ; repeat if n positive

例子12.9倒数在edx中余下的元素数，如果edx < 16，在最后一次迭代中，使用bzhi清除ecx里余下的比特。指令bzhi属于BMI2指令集，所有已知支持AVX512的处理器都支持它。注意，这仅对n < 256可工作，因为bzhi仅读edx的低8位。

如果迭代计数很高且计算掩码的指令减慢了循环的执行，主循环不带掩码，在主循环后进行要求掩码的操作，会更好：

; Example 12.10. Vectorized DAXPY, masking only after main loop

 

.code

      lea rsi, [X]                                                                                     ; point to beginning of X

      lea rdi, [Y]                                                                                     ; point to beginning of Y

      vbroadcastss zmm2, dword[DA]                                              ; broadcast DA

      mov edx, n                                                                                   ; number of elements, n

      and edx, -16                                                                                 ; round down n to multiple of 16

      lea rsi, [rsi+rdx*4]                                                                       ; point to the end

      lea rdi, [rdi+rdx*4]                                                                     ; point to the end

      neg rdx                                                                                         ; use negative index from the end

L1:                                                                                                       ; Main loop rolled out by 16

      vmovups zmm3, zmmword [rdi+rdx*4]                                 ; load Y[i]

      vfnmadd231ps zmm3, zmm2, zword [rsi+rdx*4]                 ; Y[i]-DA*X[I]

      vmovups zmmword [rdi+rdx*4] {k1}, zmm3                         ; store Y[i]

      add rdx, 16                                                                                   ; next 16 elements

      jnz L1                                                                                             ; count up to zero

      mov edx, n                                                                                   ; calculate remaining elements

      and edx, 15                                                                                  ; n modulo 16

      jz L2                                                                                               ; optionally skip if zero

      mov eax, -1                                                                                  ; all 1's

      bzhi eax, eax, edx                                                                       ; zero bit positions > edx

      kmovw k1, eax                                                                            ; copy 16 bits to mask register

      ; rsi and rdi are pointing to where the main loop ended.

      ; Now do the last 0-15 elements

      vmovups zmm3 {k1}{z}, zmmword [rdi]                                 ; load Y[i]

      vfnmadd231ps zmm3 {k1}, zmm2, zword [rsi]                      ; Y[i]-DA*X[I]

      vmovups zmmword [rdi] {k1}, zmm3                                      ; store Y[i]

L2:

例子12.10包含循环体两次。第一个主体在L1循环重复n/16次。第二个主体在n不能被16整除时，计算余下的0 ~ 15个元素。如果这个分支预测不好，没有必要使用jz L2跳转到最后的主体。

例子12.10在循环体中仅有5条指令。如果循环计数很高，这将是最快的解决方案。如果循环主体很大，我们不希望包含它两次，或者循环计数小，或者瓶颈是内存访问或循环携带依赖链，例子12.8与12.9是有用的。例子12.9仅用在n < 256时。

12.14. ​​​​​​​优化缓存

在访问非缓存内存的循环中，内存访问很可能需要最多时间。如果可能且连续访问，数据应该连续保存，如第11章77页所述。

在循环中数组访问数不应该超出微处理器中读写缓冲数量。减少数据流的一个方式是将多个数组合并为一个结构体数组，使多个数据流被交织为单个流。

现代微处理器有先进的预取机制。这些机制可以检测数据访问模式中的规律，比如以特定步长访问数据。建议尽量减少不同数据流数，尽可能保持固定步长，以利用这样的预取机制。如果数据访问模式足够规律，自动数据预取通常比显式数据预取工作得更好。

在数据访问模式不够规律，自动预取机制无法预测时，使用prefetch指令显式预取数据可能是必须的。对一个以不规律方式访问数据的程序，找出最优的预取策略，通常需要大量的实验。

如果系统有多个CPU核，可以将数据预取放入一个独立的线程。Intel C++有这样的一个特性。

使用一个2指数步长的数据访问很可能导致缓存行竞争。这可以通过改变步长或循环分块来避免。细节参考手册《优化C++软件》中优化内存访问的章节。

对写不太可能很快访问的非缓存内存，非临时写指令是有用的。可以使用向量指令尽量减少非临时写指令，

12.15. ​​​​​​​并行

提升CPU密集代码性能的最重要方式是并行。并行的主要方法有：

• 提升CPU进行乱序执行的可能性。这通过打破长依赖链（参考第59页），以及在不同执行单元或执行端口间均匀分配μop（参考第90页）来完成。
• 使用向量指令。参考第13章，113页。
• 使用多线程。参考第14章，142页。

可以使用多个累加器，打破循环携带依赖链，如第59页所述。如果CPU没有别的事情可做，最优累加器数是依赖链中最关键指令的时延除该指令的吞吐率倒数。例如，在AMD处理器上，浮点加法的时延是4时钟周期，吞吐率倒数是1。这意味着最优累加器数是4。下面的例子12.1展示看使用4个浮点寄存器作为累加器，执行加法的循环。

// Example 12.11a, Loop-carried dependency chain

// (Same as example 9.3a page 65)

double list[100], sum = 0.;

for (int i = 0; i < 100; i++) sum += list[i];

使用4个浮点寄存器作为累加器的实现，看起来像这样：

; Example 12.11b, Four floating point accumulators

      lea esi, list                                                                                       ; Pointer to list

      fld qword ptr [esi]                                                                         ; accum1 = list[0]

      fld qword ptr [esi+8]                                                                     ; accum2 = list[1]

      fld qword ptr [esi+16]                                                                   ; accum3 = list[2]

      fld qword ptr [esi+24]                                                                   ; accum4 = list[3]

      fxch st(3)                                                                                         ; Get accum1 to top

      add esi, 800                                                                                    ; Point to end of list

      mov eax, 32-800                                                                            ; Index to list[4] from end of list

L1:

      fadd qword ptr [esi+eax]                                                              ; Add list[i]

      fxch st(1)                                                                                          ; Swap accumulators

      fadd qword ptr [esi+eax+8]                                                          ; Add list[i+1]

      fxch st(2)                                                                                          ; Swap accumulators

      fadd qword ptr [esi+eax+16]                                                        ; Add list[i+2]

      fxch st(3)                                                                                          ; Swap accumulators

      add eax, 24                                                                                      ; i += 3

      js L1                                                                                                   ; Loop

      faddp st(1), st(0)                                                                             ; Add two accumulators together

      fxch st(1)                                                                                          ; Swap accumulators

      faddp st(2), st(0)                                                                             ; Add the two other accumulators

      faddp st(1), st(0)                                                                             ; Add these sums

      fstp qword ptr [sum]                                                                      ; Store the result

在例子12.11中，从list向4个累加器载入前4个值。

在循环中执行的加法数恰好被展开因子3整除。使用浮点寄存器作为累加器有趣的地方是，累加器数等于展开因子加1。这是交换累加器的fxch指令的使用方式的结果。你必须密切注意浮点寄存器栈上的每个累加器，验证在循环的每次迭代后，4个累加器旋转一个位置，使每个累加器每4个加法得到使用，虽然循环仅展开3次。

例子12.11b中的循环每次加法需要1时钟周期，这是浮点加法器的最大吞吐率。浮点加法4时钟周期的时延，由4个累加器处理。指令fxch有零时延，因为在大多数Intel、AMD与VIA处理器（除了Intel Atom）上，它们被翻译为寄存器重命名。

在具有SSE2指令集的处理器上，通过使用XMM寄存器，而不是浮点栈寄存器，可以避免fxch指令，如例子12.11c所示。在AMD上，浮点向量加法的时延是4，吞吐率倒数是2，因此，最优累加器数是2个向量寄存器。

; Example 12.11c, Two XMM vector accumulators

      lea esi, list                                                                         ; list must be aligned by 16

      movapd xmm0, [esi]                                                       ; list[0], list[1]

      movapd xmm1, [esi+16]                                                ; list[2], list[3]

      add esi, 800                                                                      ; Point to end of list

      mov eax, 32-800                                                              ; Index to list[4] from end of list

L1:

      addpd xmm0, [esi+eax]                                                  ; Add list[i], list[i+1]

      addpd xmm1, [esi+eax+16]                                           ; Add list[i+2], list[i+3]

      add eax, 32                                                                       ; i += 4

      js L1                                                                                    ; Loop

      addpd xmm0, xmm1                                                       ; Add the two accumulators together

      movhlps xmm1, xmm0                                                   ; There is no movhlpd instruction

      addsd xmm0, xmm1                                                       ; Add the two vector elements

      movsd [sum], xmm0                                                       ; Store the result

在AMD处理器上，例子12.11b与12.11c一样快，因为两者都受到浮点加法器吞吐率的限制。

在Intel Core2处理器上，例子12.11c比12.11b更快，因为这个处理器有可以在一次操作中处理整个向量的128位浮点加法器、

12.16. ​​​​​​​分析依赖性

循环可能有几个互锁的依赖链。这样复杂的情形要求仔细的分析。

下一个例子是泰勒展开。你可能知道，许多函数可通过这个形式的泰勒多项式逼近

fx≈ i=0ncixi

每个指数xi通过之前指数xi-1乘以x，方便地计算。系数ci保存在一个表中。

; Example 12.12a. Taylor expansion

.data

      align 16

      one dq 1.0                                                                                ; 1.0

      x dq ?                                                                                        ; x

      coeff dq c0, c1, c2, ...                                                             ; Taylor coefficients

      coeff_end label qword                                                          ; end of coeff. list

 

.code

      movsd xmm2, [x]                                                                   ; xmm2 = x

      movsd xmm1, [one]                                                              ; xmm1 = x^i

      xorps xmm0, xmm0                                                              ; xmm0 = sum. init. to 0

      mov eax, offset coeff                                                            ; point to c[i]

L1: movsd xmm3, [eax]                                                              ; c[i]

      mulsd xmm3, xmm1                                                             ; c[i] * x^i

      mulsd xmm1, xmm2                                                             ; x^(i+1)

      addsd xmm0, xmm3                                                             ; sum += c[i] * x^i

      add eax, 8                                                                               ; point to c[i+1]

      cmp eax, offset coeff_end                                                   ; stop at end of list

      jb L1                                                                                         ; loop

（如果你的汇编器混淆了movsd指令与同名的字符串指令，把它编码为DB 0F2H / movups）。

现在关注分析。系数列表很短，我们可以预期它待在缓存中。

为了检查时延是否重要，我们必须查看代码中的依赖。这些依赖显示在图12.1中。

图12.1. 例子12.12泰勒展开中的依赖

有两个连续的依赖链，一个用于计算xi，一个用于计算和。在Intel处理器上，指令mulsd的时延比addsd长。因此，垂直乘法链比加法链更关键。加法必须等待cixi，这个值比xi晚一个乘法时延，且比前面的加法晚。如果没有别的限制性能，那么可以预期这个代码每迭代需要一个乘法时延，这是4 ~ 7时钟周期，取决于处理器。

吞吐率不是限制因素，因为在具有128位执行单元的处理器上，乘法单元的吞吐率是每时钟周期一次乘法。

测量例子12.12a的时间，我们发现循环比乘法多用了至少一个时钟周期。对此解释如下：两个乘法都要等待前面等待在xmm1中的值xi-1。因此，两个乘法在同一时间就绪。我们希望图12.1中阶梯里的垂直乘法先开始，因为它是最关键依赖链的部分。但是微处理器看不到交换这两个乘法的理由，因此图12.1中的水平乘法先开始。垂直乘法被推迟1时钟周期，这是浮点乘法单元的吞吐率倒数。可以通过把图12.1中的水平乘法放在垂直乘法后，可以解决这个问题：

; Example 12.12b. Taylor expansion

      movsd xmm2, [x]                                                             ; xmm2 = x

      movsd xmm1, [one]                                                        ; xmm1 = x^i

      xorps xmm0, xmm0                                                        ; xmm0 = sum. initialize to 0

      mov eax, offset coeff                                                      ; point to c[i]

L1: movapd xmm3, xmm1                                                   ; copy x^i

      mulsd xmm1, xmm2                                                       ; x^(i+1) (vertical multipl.)

      mulsd xmm3, [eax]                                                         ; c[i]*x^i (horizontal multipl.)

      add eax, 8                                                                         ; point to c[i+1]

      cmp eax, offset coeff_end                                             ; stop at end of list

      addsd xmm0, xmm3                                                       ; sum += c[i] * x^i

      jb L1                                                                                   ; loop

这里，我们需要额外的指令来拷贝xi（movapd xmm3, xmm1），因此我们可以在水平乘法前进行垂直乘法。这使得循环每迭代的执行快1时钟周期。

我们仍然仅使用了xmm寄存器低半部。同时执行两个乘法，一次进行泰勒展开的两步，增益更大。技巧是从xi-2乘上x2计算xi的每个值：

; Example 12.12c. Taylor expansion, double steps

      movsd xmm4, [x]                                                                         ; x

      mulsd xmm4, xmm4                                                                   ; x^2

      movlhps xmm4, xmm4                                                               ; x^2, x^2

      movapd xmm1, xmmword ptr [one]                                        ; xmm1(L)=1.0, xmm1(H)=x

      xorps xmm0, xmm0                                                                     ; xmm0 = sum. init. to 0

      mov eax, offset coeff                                                                   ; point to c[i]

L1: movapd xmm3, xmm1                                                                ; copy x^i, x^(i+1)

      mulpd xmm1, xmm4                                                                    ; x^(i+2), x^(i+3)

      mulpd xmm3, [eax]                                                                      ; c[i]*x^i, c[i+1]*x^(i+1)

      addpd xmm0, xmm3                                                                    ; sum += c[i] * x^i

      add eax, 16                                                                                    ; point to c[i+2]

      cmp eax, offset coeff_end                                                          ; stop at end of list

      jb L1                                                                                                ; loop

      haddpd xmm0, xmm0                                                                 ; join the two parts of sum

现在速度加倍了。循环每迭代仍然需要乘法时延，但迭代次数减半。乘法时延是4 ~ 5时钟周期，取决于处理器。每迭代执行两个向量乘法，我们仅使用了最大乘法器吞吐率的40或50%。这意味着提升并行度，有更大的增益。可以一次进行4步泰勒展开，xi-4乘上x4计算xi每个值来做到：

; Example 12.12d. Taylor expansion, quadruple steps

      movsd xmm4, [x]                                                                   ; x

      mulsd xmm4, xmm4                                                             ; x^2

      movlhps xmm4, xmm4                                                         ; x^2, x^2

      movapd xmm2, xmmword ptr [one]                                  ; xmm2(L)=1.0, xmm2(H)=x

      movapd xmm1, xmm2                                                          ; xmm1 = 1, x

      mulpd xmm2, xmm4                                                             ; xmm2 = x^2, x^3

      mulpd xmm4, xmm4                                                             ; xmm4 = x^4, x^4

      xorps xmm5, xmm5                                                               ; xmm5 = sum. init. to 0

      xorps xmm6, xmm6                                                               ; xmm6 = sum. init. to 0

      mov eax, offset coeff                                                             ; point to c[i]

L1: movapd xmm3, xmm1                                                          ; copy x^i, x^(i+1)

      movapd xmm0, xmm2                                                          ; copy x^(i+2), x^(i+3)

      mulpd xmm1, xmm4                                                             ; x^(i+4), x^(i+5)

      mulpd xmm2, xmm4                                                             ; x^(i+6), x^(i+7)

      mulpd xmm3, xmmword ptr [eax]                                      ; term(i), term(i+1)

      mulpd xmm0, xmmword ptr [eax+16]                               ; term(i+2), term(i+3)

      addpd xmm5, xmm3                                                             ; add to sum

      addpd xmm6, xmm0                                                             ; add to sum

      add eax, 32                                                                             ; point to c[i+2]

      cmp eax, offset coeff_end                                                   ; stop at end of list

      jb L1                                                                                         ; loop

      addpd xmm5, xmm6                                                            ; join two accumulators

      haddpd xmm5, xmm5                                                          ; final sum in xmm5

在例子12.12d中，我们使用两个寄存器xmm1与xmm2来保存x的连续指数，以将图12.1中的垂直依赖链分解为4个并行的链。每个xi指数从xi-4计算。类似的，我们对和使用两个累加器，xmm5与xmm6，以将加法链分解为4条并行的链。在循环后，四个和被加起来。在Core2上的一个测量显示，对例子12.12d，每迭代5.04时钟周期，或每泰勒项1.26时钟周期。这接近于由乘法时延确定的、每迭代5时钟周期的期望。这个解决方案使用了乘法器最大吞吐率的80%，这是一个非常令人满意的结果。进一步展开几乎没有增益。向量加法与其他指令不会增加执行时间，因为它们不需要使用与乘法相同的执行单元。在这个情形里，指令获取与解码不是瓶颈。

在带有256位YMM寄存器、支持AVX指令集的处理器上，吞吐率可以进一步提高：

; Example 12.12e. Taylor expansion using AVX

      vmovddup xmm2, [x]                                                             ; x, x

      vmulpd xmm5, xmm2, xmm2                                               ; x^2, x^2

      vmulpd xmm4, xmm5, xmm5                                               ; x^4, x^4

      vmovsd xmm2, [one]                                                              ; 1

      vmovhpd xmm2, xmm2, [x]                                                   ; 1, x

      vmulpd xmm0, xmm2, xmm5                                               ; x^2, x^3

      vinsertf128 ymm2, ymm2, xmm0, 1                                    ; 1, x, x^2, x^3

      vinsertf128 ymm4, ymm4, xmm4, 1                                    ; x^4, x^4, x^4, x^4

      vmulpd ymm3, ymm2, ymm4                                               ; x^4, x^5, x^6, x^7

      vmulpd ymm4, ymm4, ymm4                                               ; x^8, x^8, x^8, x^8

      vxorps xmm0, xmm0                                                              ; sum0 = 0

      vxorps xmm1, xmm1                                                              ; sum1 = 0

      lea rax, [coeff_end]                                                                 ; point to end of coeff[i]

      mov rcx, -n                                                                                ; counter

      jmp L2                                                                                        ; jump into loop

      align 32                                                                                      ; loop

L1: vmulpd ymm2, ymm2, ymm4                                               ; multiply powers of x by x^8

      vmulpd ymm3, ymm3, ymm4                                               ; multiply powers of x by x^8

L2: vmulpd ymm5, ymm2, [rax+rcx]                                          ; first four terms

      vaddpd ymm0, ymm0, ymm5                                               ; sum0

      vmulpd ymm5, ymm3, [rax+rcx+32]                                    ; next four terms

      vaddpd ymm1, ymm1, ymm5                                               ; sum1

      add rcx, 64

      jl L1

      ; make total sum

      vaddpd ymm0, ymm0, ymm1                                               ; join two accumulators

      vextractf128 xmm5, ymm0, 1                                               ; get high part

      vaddpd xmm0, xmm0, xmm5                                               ; add high and low part

      vhaddpd xmm0, xmm0, xmm0                                             ; final sum in xmm0

      vzeroupper                                                                               ; if subsequent code is non-VEX

例子12.12e每迭代计算8项。在Sandy Bridge上，测得执行时间为每迭代5.5时钟周期，接近由5时钟周期乘法时延确定的理论最小值。在Skylake上，测得执行时间大致相同，尽管在这个处理器上乘法时延仅是4时钟周期。

当然，在例子12.12d与12.12e中，花在循环外的时间比前面的例子多，但性能增益仍然可观，即使对于不多的迭代。最后的一小点别扭是，在循环前将第一个系数放入累加器，从x1开始，而不是x0。可以通过最后加第一项得到更好的精度，但代价是一个额外的加法时延。

可以在循环中使用融合乘加指令（FMA3或FMA4），进一步改进例子12.12e的代码，如果支持这些指令。在Skylake上测得，在使用FMA指令时，执行时间从5.4减为4.5时钟周期。

总之，垂直乘法应该在水平乘法前完成。如果循环迭代次数很高，最佳累加器寄存器数是使用执行单元最大可能吞吐率的值。如果循环迭代次数很低，循环前的设置时间与循环后的加法都要紧，那么最佳累加器寄存器数可能要小一些。

一个泰勒展开通常在项变成负数时结束。不过，总是包括最坏情形的最大项数，以避免检查停止条件所需的浮点比较，可能更好。重复次数常量也将防止循环退出的误预测。计算比所需更多的项的代价，在每项可能使用不到1时钟周期的优化例子中实际上相当小。

记住将mxcsr寄存器置为“Flush to zero模式”，以避免在执行比所需更多次迭代时，可能下溢导致的性能损失。

1. 1. 没有乱序执行处理器上的循环

P1与PMMX处理器没有乱序执行能力。当然，这些处理器今天已经过时了，但这里描述的原则可能适用于更小的嵌入式处理器。

我选择了一个简单的例子——一个从数组读整数，改变每个整数的符号，把结果保存到另一个数组的过程。这个过程的C++语言代码应该是：

// Example 12.13a. Loop to change sign

void ChangeSign (int * A, int * B, int N) {

     for (int i = 0; i < N; i++) B[i] = -A[i];

}

一个对P1优化的汇编实现看起来像这样：

; Example 12.13b

      mov esi, [A]

      mov eax, [N]

      mov edi, [B]

      xor ecx, ecx

      lea esi, [esi+4*eax]                                                 ; point to end of array a

      sub ecx, eax                                                             ; -n

      lea edi, [edi+4*eax]                                                ; point to end of array b

      jz short L3

      xor ebx, ebx                                                             ; start first calculation

      mov eax, [esi+4*ecx]

      inc ecx

      jz short L2

L1: sub ebx, eax                                                             ; u

      mov eax, [esi+4*ecx]                                              ; v (pairs)

      mov [edi+4*ecx-4], ebx                                          ; u

      inc ecx                                                                        ; v (pairs)

      mov ebx, 0                                                                ; u

      jnz L1                                                                          ; v (pairs)

L2: sub ebx, eax                                                              ; end last calculation

      mov [edi+4*ecx-4], ebx

L3:

这里，迭代是重叠的，以提升成对的机会。在保存第一个值之前，开始读第二个值。指令mov ebx, 0放在inc ecx与jnz L1之间，不是改进成对，而是第一个指令对使用ecx作为地址索引，避免在递增ecx后导致的AGI暂停。

带有浮点操作的循环有点不同，因为浮点指令是流水线化重叠的，而不是成对。考虑第90页例子12.6的DAXPY循环。对P1优化的一个解决方案如下：

; Example 12.14. DAXPY optimized for P1 and PMMX

      mov eax, [n]                                                             ; number of elements

      mov esi, [X]                                                              ; pointer to X

      mov edi, [Y]                                                              ; pointer to Y

      xor ecx, ecx

      lea esi, [esi+8*eax]                                                  ; point to end of X

      sub ecx, eax                                                              ; -n

      lea edi, [edi+8*eax]                                                 ; point to end of Y

      jz short L3                                                                  ; test for n = 0

      fld qword ptr [DA]                                                   ; start first calc.

      fmul qword ptr [esi+8*ecx]                                   ; DA * X[0]

      jmp short L2                                                             ; jump into loop

L1: fld qword ptr [DA]

      fmul qword ptr [esi+8*ecx]                                   ; DA * X[i]

      fxch                                                                            ; get old result

      fstp qword ptr [edi+8*ecx-8]                                ; store Y[i]

L2: fsubr qword ptr [edi+8*ecx]                                 ; subtract from Y[i]

      inc ecx                                                                       ; increment index

      jnz L1                                                                         ; loop

      fstp qword ptr [edi+8*ecx-8]                                ; store last result

L3:

这里我们将循环计数器用作数组索引，从负数数到0。

每个操作在之前一个结束前开始，以提升计算重叠。具有乱序执行的处理器自动这样做，但对没有乱序能力的处理器，我们必须显式进行这个重叠。

这里，浮点操作的交错完美工作：之前结果的FSTP填充到FMUL与FSUBR间2时钟的暂停。循环开销与下一个操作的前两条指令填充到FSUBR与FSTP之间3时钟暂停。在递增索引后的第一个时钟周期中读不依赖于索引计数器参数，避免了地址生成互锁（AGI）暂停。

这个解决方案每迭代需要6时钟周期，比别处发布的展开解决方案要好。

12.18. ​​​​​​​宏循环

如果一个循环的重复计数小且为常量，完全展开该循环是可能的。这样的好处是，仅依赖于循环计数器的计算可以在汇编时刻完成，而不是在执行时刻。当然坏处是，如果重复次数很高，它占据更多的代码缓存。

MASM语法提供了一个完成支持元编程的强大宏语言，它相当有用。其他汇编器有类似的能力，但语法不同。例如，如果我们需要一个平方数列表，C++代码看起来像这样：

// Example 12.15a. Loop to make list of squares

int squares[10];

for (int i = 0; i < 10; i++) squares[i] = i*i;

在MASM语言中可以通过一个宏循环产生相同的列表：

; Example 12.15b. Macro loop to produce data

.DATA

squares LABEL DWORD                                       ; label at start of array

I = 0                                                                         ; temporary counter

REPT 10                                                                  ; repeat 10 times

      DD I * I                                                              ; define one array element

      I = I + 1                                                              ; increment counter

ENDM                                                                     ; end of REPT loop

这里，I是一个预处理变量。I循环在汇编时刻运行，但不在执行时刻。变量I与语句I = I + 1不会进入最终的代码，因此无需时间执行。实际上，例子12.15b不产生可执行代码，仅数据。宏预处理器将上面的代码翻译为：

; Example 12.15c. Resuls of macro loop expansion

squares LABEL DWORD ; label at start of array

      DD 0

      DD 1

      DD 4

      DD 9

      DD 16

      DD 25

      DD 36

      DD 49

      DD 64

      DD 81

宏循环对产生代码也是有用的。下一个例子计算xn，其中x是一个浮点值，n是一个正整数。通过反复取x平方，与n里对应二进制数字的因子相乘，这可以最高效地执行。这个算法的C++代码：

// Example 12.16a. Calculate pow(x,n) where n is a positive integer

double x, xp, power;

unsigned int n, i;

xp = x; power = 1.0;

for (i = n; i != 0; i >>= 1) {

      if (i & 1) power *= xp;

      xp *= xp;

}

如果n在汇编时刻已知，那么可以使用性能的宏循环实现这个指数函数：

; Example 12.16b.

; This macro will raise a double-precision float in X

; to the power of N, where N is a positive integer constant.

; The result is returned in Y. X and Y must be two different

; XMM registers. X is not preserved.

; (Only for processors with SSE2)

INTPOWER MACRO X, Y, N

      LOCAL I, YUSED                                                                   ; define local identifiers

      I = N                                                                                       ; I used for shifting N

      YUSED = 0                                                                             ; remember if Y contains valid data

      REPT 32                                                                                 ; maximum repeat count is 32

            IF I AND 1                                                                        ; test bit 0

                  IF YUSED                                                                    ; If Y already contains data

                        mulsd Y, X                                                            ; multiply Y with a power of X

                  ELSE                                                                            ; If this is first time Y is used:

                        movsd Y, X                                                           ; copy data to Y

                        YUSED = 1                                                            ; remember that Y now contains data

                  ENDIF                                                                          ; end of IF YUSED

            ENDIF                                                                                ; end of IF I AND 1

            I = I SHR 1                                                                         ; shift right I one place

            IF I EQ 0                                                                            ; stop when I = 0

                  EXITM                                                                         ; exit REPT 32 loop prematurely

            ENDIF                                                                                ; end of IF I EQ 0

            mulsd X, X                                                                        ; square X

      ENDM                                                                                     ; end of REPT 32 loop

ENDM                                                                                           ; end of INTPOWER macro definition

这个宏产生了完成这项工作所需的最少指令。在最终代码中没有循环开销、prolog与epilog。最重要的，没有分支。所有的分支都被宏预处理器解决了。要计算xmm0的12次方，写：

; Example 12.16c. Macro invocation

INTPOWER xmm0, xmm1, 12

这将被展开为：

; Example 12.16d. Result of macro expansion

mulsd xmm0, xmm0                                                                   ; x^2

mulsd xmm0, xmm0                                                                   ; x^4

movsd xmm1, xmm0                                                                  ; save x^4

mulsd xmm0, xmm0                                                                   ; x^8

mulsd xmm1, xmm0                                                                   ; x^4 * x^8 = x^12

这甚至比没有展开的优化汇编循环的指令还要少。当mulpd替换mulsd，movapd替换movsd时，这个宏还可以在向量上工作。