使用 MMX 指令集计算 f(x)=2*(x^2)+5 饱和度，用于从二进制文件加载的 128 个大小为 2 字节的数字
Posted
技术标签:
【中文标题】使用 MMX 指令集计算 f(x)=2*(x^2)+5 饱和度，用于从二进制文件加载的 128 个大小为 2 字节的数字【英文标题】：Calculate f(x)=2*(x^2)+5 with saturation using MMX instruction set for 128 numbers with size of 2 bytes loaded from a binary file 【发布时间】：2019-01-17 16:03:59 【问题描述】：我有这个问题，我需要使用 MMX 指令集计算函数 f(x)=2*(x^2)+5。我有两个问题。这是我现在的代码：
section .data
print_fmt db '%d', 10, 0
my_loaded_data times 128 dw 0
fives times 4 dw 5
twos times 4 dw 2
my_loaded_data_file_name db 'test_numbers.bin', 0
mod_f db 'r', 0
section .text
extern printf, fopen, fread
global main
main:
PUSH rbp
mov rbp, rsp
mov rax, 0
mov rdi, my_loaded_data_file_name
mov rsi, mod_f
call fopen
cmp rax, 0
je .end
PUSH rax
PUSH rdi
PUSH rsi
PUSH rdx
PUSH rcx
mov rdi, my_loaded_data
mov rsi, 2
mov rdx, 128
mov rcx, rax
mov rax, 0
call fread
POP rcx
POP rdx
POP rsi
POP rdi
POP rax
mov rsi, my_loaded_data
mov rcx, 32
jmp .square_loop
.square_loop:
movq mm0, [rsi]
movq mm1, [rsi]
pmulhw mm0, mm1
movq [rsi], mm0
add rsi, 8
loop .square_loop
mov rcx, 32
mov rsi, my_loaded_data
movq mm1, [twos]
jmp .mult_with2_loop
.mult_with2_loop:
movq mm0, [rsi]
pmulhw mm0, mm1
movq [rsi], mm0
add rsi, 8
loop .mult_with2_loop
mov rcx, 32
mov rsi, my_loaded_data
movq mm1, [fives]
jmp .add_five_loop
.add_five_loop:
movq mm0, [rsi]
paddusw mm0, mm1
movq [rsi], mm0
add rsi, 8
loop .add_five_loop
jmp .write_it
.write_it:
mov r8, my_loaded_data
mov rcx, 128
.write_loop:
mov rax, 0
mov ax, [r8]
PUSH r8
PUSH rcx
PUSH rdi
PUSH rsi
PUSH rax
mov rdi, print_fmt
mov rsi, rax
mov rax, 0
call printf
POP rax
POP rsi
POP rdi
POP rcx
POP r8
add r8, 2
loop .write_loop
.end:
mov rax, 0
POP rbp
ret
我的第一个问题是乘法指令。我使用哪个指令来进行饱和。首先，我虽然会有像pmulsw 这样的指令，但似乎没有。 pmulhw 保存 32 位结果的高 16 位。我找不到任何可以给出 16 位结果的指令。这是保存 32 位结果的唯一方法吗？
第二个问题是 printf。它不断给出分段错误，我不知道为什么。这是来自我的终端：
Program received signal SIGSEGV, Segmentation fault.
__GI___tcgetattr (fd=1, termios_p=termios_p@entry=0x7ffffffed9a8) at ../sysdeps/unix/sysv/linux/tcgetattr.c:42
42 ../sysdeps/unix/sysv/linux/tcgetattr.c: No such file or directory.
这是生成文件：
zad7.o: zad7.asm
nasm -f elf64 -g -F dwarf zad7.asm -o zad7.o
zad7: zad7.o
gcc -o zad7 zad7.o -no-pie -ggdb
为了您的方便，这里有一个小 C 程序，它可以生成二进制文件以供阅读：
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
void make_data_file(unsigned int number_of_data, unsigned int size_in_bytes, char* file_name)
FILE *write_ptr = fopen(file_name, "wb");
if(write_ptr==NULL)
printf("Error creating file '%s'.\n", file_name);
return;
double n_bits = size_in_bytes * 8;
unsigned int max_num = pow(2, n_bits);
unsigned short random_number;
for(int i=0; i< number_of_data; i++)
random_number = i;
fwrite(&random_number, size_in_bytes, 1, write_ptr);
fclose(write_ptr);
int main()
make_data_file(128, 2, "test_numbers.bin");
return 0;
【问题讨论】：
printf 通常是未对齐堆栈的段错误。请注意，如果您使用的是 x86-64，那么您有 SSE2，因此使用 MMX 有点毫无意义。
@Jester 未对齐的堆栈到底是什么意思？如果我只推送 8 字节的值，它怎么会错位呢？我如何解决它？这是我正在学习的课程的练习任务，所以他们希望我使用 MMX。
【参考方案1】：
如果您关心性能，请避免在现代 CPU 上使用 loop 指令。 Why is the loop instruction slow? Couldn't Intel have implemented it efficiently?。也使用 SSE2 代替 MMX；您的数组大小是 16 和 8 的倍数，并且您使用的是 x86-64，它保证具有 SSE2。 MMX 对此毫无意义，除非您还为 Pentium III / Athlon-XP 及更早版本制作 32 位版本。
（我的答案中的所有代码都将使用 8 字节 MMX 寄存器而不是 16 字节 XMM 寄存器，因为我使用的所有指令都有 MMX 版本。根据appendix B to the NASM manual、pmullw、@ 987654328@、pcmpgtw 和 paddusw 都在原始 P5 Pentium MMX 中可用。英特尔手册中列为“MMX”的一些指令（如 pmulhuw 和 pshufw）只是后来添加的，就像 Pentium II 可能，或者与 Pentium III 中的 SSE 一起使用，但对于此处有用的任何指令，情况并非如此。）
请参阅https://***.com/tags/x86/info 了解性能/优化指南，以及解释调用函数所需的 16 字节堆栈对齐的 ABI/调用约定链接。
mov rax, 0 / mov ax, [r8] 真的很傻。像普通人一样使用movzx eax, word [r8]。您也不需要 jmp 到下一个源代码行，例如 jmp .square_loop / .square_loop: 如果没有分支指令，执行总是会自行落到下一行。
x86 SIMD 没有饱和乘法，只有饱和有符号/无符号加法和饱和打包到更窄的元素。 （MMX/SSE2 paddsw / paddusw）。既然您使用%d 打印，也许您想要签名饱和度？但这只是在解压缩为 32 位之后，您的公式将始终产生积极的结果，因此您可以使用无符号饱和度。我知道这就是您的代码使用的paddusw。
此外，在公式的每个步骤之间使用 3 个单独的循环来存储/重新加载您的数据真的很糟糕。您（几乎）总是希望增加计算强度（每个内存/缓存带宽对数据完成的 ALU 工作量）。您也不需要乘法指令来将数字加倍：只需将其添加到自身。 padd* 比 pmul* 在更多端口上运行，并且具有更好的延迟和（在这种情况下相关的）吞吐量。
default rel
;;; Efficient but *doesn't* saturate the multiply
lea rcx, [rsi + length] ; rcx = end_pointer
movdqa xmm5, [fives]
.loop: ; do
movdqu xmm0, [rsi] ; use movdqa if your data is aligned, faster on very old CPUs.
pmullw xmm0, xmm0 ; x*x ; does NOT saturate. will wrap.
paddusw xmm0, xmm0 ; *= 2 ; with unsigned saturation
paddusw xmm0, xmm5 ; += 5 ; xmm5 = _mm_set1_epi16(5) outside the loop
movdqu [rsi], xmm0
add rsi, 16
cmp rsi, rcx ; while(p<endp);
jb .loop
...
section .rodata
align 16
fives: times 8 dw 5
对于饱和度，您也许可以使用 SSSE3 https://www.felixcloutier.com/x86/pmaddubsw，但这只需要字节输入。它使 i8 x u8 => i16 乘积对的水平总和饱和。
否则，您可能必须解压缩为 dwords 和 packssdw（有符号）或 packusdw（无符号饱和度）返回单词。但是 SSE4.1 pmulld 的双字乘法很慢（Haswell 和更高版本为 2 微秒）。不过，在一些较旧的 CPU 上实际上只有 1 uop。当然，拆包会因拥有更广泛的元素而产生两倍的工作量。
在这种情况下，您的公式与输入的大小是单调的，因此您只需对输入进行比较并手动饱和。
如果我们假设您的输入也是无符号的，那么我们不必在比较之前做一个绝对值。但是（直到 AVX512）我们没有无符号整数比较，只有有符号大于，所以大的无符号输入将比较为负。
2 * 0x00b5^2 + 5 = 0xfff7 适合 16 位
2 * 0x00b6^2 + 5 = 0x102cd 没有，我们希望它饱和到 0xffff
溢出截止点是偶数，因此我们可以通过右移来处理有符号比较问题。那将是无符号除以 2，让 use 安全地将结果视为非负有符号整数。 0xb6 >> 1 = 0x5b。但是pcmpgtw 是> 的比较对象，而不是>=。
如果我们将操作数反转为 pcmpgtw 以比较 (x>>1) < (0xb6>>1)，那么我们必须将常量 movdqa 以避免破坏它，但我们仍然需要使用 movdqa+ 右移 x psrlw。当需要饱和度（否则为 0）时，使用 0xffff 的向量会更有效，因为我们可以使用 POR 或 PADDUSW 直接应用它。
所以我们最好的选择是简单地将 x 和 0xb5 范围转移到有符号，并使用 pcmpgtw SIMD 符号比较来执行 (x-0x8000) > (0xb5 - 0x8000)。
其他更糟糕的选择包括：
检查与pmulhuw 相乘的溢出以计算高半部分（并检查它是否非零）。我们将面临乘数吞吐量遇到瓶颈的危险，而使用 pcmpeqw 检查非零值与我们想要的条件相反。
psubusw x, 0xb5 并检查它是否 == 0。pcmpeqw 会给我们一个反转掩码，但我们不能使用 pcmpgtw 来检查 usat16(x-0xb5) > 0因为大输入（设置高位）将保持“负数” " 减去 0xb5 这样的小数后。
paddusw 并检查 == 0xffff：只有足够小的输入才会使最终结果不饱和。一些 CPU 可以在比padd* 更多的端口上运行pxor，并且它不需要任何更少的非零向量常量，所以这无论如何都不是更好。但在 Skylake 上也同样出色。
default rel
;;; With a check on the input to saturate the output
lea rcx, [rsi + length] ; rcx = end_pointer
movdqa xmm4, [input_saturation_cutoff]
movdqa xmm5, [fives]
pcmpeqd xmm3, xmm3
psllw xmm3, 15 ; xmm3 = set1_epi16(0x8000) for range-shifting to signed
.loop:
movdqu xmm0, [rsi]
movdqa xmm1, xmm0
; if x>0xb5 (unsigned), saturate output to 0xffff
pxor xmm1, xmm3 ; x - 0x8000 = range shift to signed for compare
pcmpgtw xmm1, xmm4 ; xmm1 = (x > 0xb5) ? -1 : 0
pmullw xmm0, xmm0 ; x*x
paddusw xmm0, xmm0 ; *= 2 ; saturating add or not doesn't matter here
por xmm1, xmm5 ; 0xffff (saturation needed) else 5. Off the critical path to give OoO exec an easier time.
paddusw xmm0, xmm1 ; += 5 or += 0xffff with unsigned saturation.
movdqu [rsi], xmm0
add rsi, 16
cmp rsi, rcx
jb .loop
...
section .rodata
align 16
input_saturation_cutoff: times 8 dw (0x00b5 - 0x8000) ; range-shifted to signed for pcmpgtw
fives: times 8 dw 5
; 5 = 0xb6 >> 5 or 0xb6 >> 5 but we'd need to knock off the high bit from input_saturation_cutoff
; Or we could materialize constants like this:
; movdqa xmm4, [set1w_0xb5]
; pcmpeqd xmm3, xmm3
; psllw xmm3, 15 ; rangeshift_mask = set1(0x8000)
; movdqa xmm5, xmm4
; psrlw xmm5, 5 ; fives = set1(5)
; pxor xmm4, xmm3 ; input_sat_cutoff = set1(0x80b5)
;; probably not worth it since we already need to load 1 from memory.
我对此进行了测试，例如 paddusw 确实可以做到 0x2 + 0xffff = 0xffff。
我们可以只用 0 或 0xffff 对最终结果进行 POR 以使其保持不变或将其设置为 0xffff，但将输入修改为最后一个 paddusw 会在一次迭代中创建更多指令级并行性。所以乱序执行不必重叠尽可能多的独立迭代来隐藏循环体的延迟。 （如果我们实际上是为有序的 Atom 或 P5 Pentium-MMX 安排这个，我们会交错更多的两个 dep 链。）
实际上，右移 1 会起作用：我们只需要比较来捕获如此大的输入，以至于乘法换行到一个小的结果。 0xb6 * 0xb6 不会换行，所以它会从paddubsw 自行饱和。
如果我们检查 (x>>1) > (0xb6>>1) 和 pcmpgtw（而不是 >=）来捕获像 0xffff 这样的输入（其中 pmullw 和 0xffff 给我们 0x0001），那就没问题了。所以我们可以保存 1 个向量常量，否则这不是更好。
pxor + pcmpgtw 至少和psrlw xmm, 1 + pcmpgtw 一样便宜，除非我们正在调整英特尔 P6 系列（Core2/Nehalem）并遇到寄存器读取端口停顿。但这不太可能：xmm0、xmm1 和 rsi 应该始终是热的（最近写入并因此从 ROB 中读取，而不是从永久寄存器文件中读取）。我们在循环中的第一组 4 条指令中只读取了 2 个向量常量，之后又读取了 1 个。
正如我在下面所说的，在许多 Intel CPU 上，psrlw 只能在与pmullw 相同的端口上运行，在 vec-int shift+multiply 执行单元上。不过，这里可能不是吞吐量瓶颈。
但是pcmp 和padd 可以在有限的端口上运行（在 Skylake 之前的英特尔上），而pxor 可以在任何矢量 ALU 端口上运行。 纯padd/@ 的混合987654395@/pmul/psrl` 微指令会留下一个向量 ALU 端口未使用。
交替饱和度检查思路
（当我寻找不会溢出的最高输入时，我忘记了公式中的 *2。）
如果公式是(0x00ff)^2 + 5，饱和度检查会更简单。
我们可以只检查位位置。
(0x00ff)^2 + 5 = 0xfe06 适合 16 位
(0x0100)^2 + 5 = 0x10005 没有，我们希望它饱和到 0xffff
所以我们需要检查高16位是否全为零，或者x&0xFF == x，或者x>>8 == 0。
这需要更少的常量，但实际上比使用 PXOR 进行有符号范围移位更糟糕，因为在某些 CPU 上，向量移位和向量乘法执行单元位于同一端口上。 （因此psrlw 和pmullw 相互竞争吞吐量。这是足够的总微指令，我们实际上不会在 Nehalem / Sandybridge / Haswell 的端口 0 上出现瓶颈，但它不会受到伤害。）
lea rcx, [rsi + length] ; rcx = end_pointer
movq xmm5, [fives]
punpcklqdq xmm5, xmm5 ; or with SSE3, movddup xmm5, [fives] to broadcast-load
pxor xmm4, xmm4 ; xmm4 = 0
.loop:
movdqu xmm0, [rsi]
movdqa xmm1, xmm0
; if x>0xffU, i.e. if the high byte >0, saturate output to 0xffff
psrlw xmm1, 8 ; x>>8 (logical right shift)
pcmpgtw xmm1, xmm4 ; xmm1 = ((x>>8) > 0) ? -1 : 0
pmullw xmm0, xmm0 ; x*x ; does NOT saturate. will wrap.
paddusw xmm0, xmm0 ; *= 2 ; with unsigned saturation
por xmm1, xmm5 ; 0xffff (saturation needed) or 5 (normal). Off the critical path to give OoO exec an easier time.
paddusw xmm0, xmm1 ; += 5 or += 0xffff with unsigned saturation.
movdqu [rsi], xmm0
add rsi, 16
cmp rsi, rcx
jb .loop
使用 AVX512BW (Skylake-X) 进行无符号比较，使用屏蔽寄存器
我们终于可以用 AVX512F 做无符号整数比较，用 AVX512BW 做字元大小比较。但是结果是在一个掩码寄存器而不是一个向量中，所以我们不能只用vpor 它和set1(5) 的向量来创建一个用于饱和加法的输入。
相反，我们可以根据比较掩码混合5 和0xffff 的向量。
;; AVX512BW version
;; On a Skylake Xeon you probably only want to use YMM regs unless you spend a lot of time in this
;; to avoid reducing max turbo much.
;; Even with xmm or ymm regs (AVX512VL + BW), this demonstrates
;; that we gain even more efficiency than just widening the vectors
;; Just having non-destructive AVX encodings saves the `movdqa xmm1,xmm0` in the SSE2 version.
;; With YMM or XMM regs, most of these instructions can still use shorter VEX encoding (AVX), not the longer EVEX (AVX512)
;; (Use vmovdqa/u instead of vmovdqu64. The 64 is element-size, irrelevant when not masking.)
;;;;;;;;;;; UNTESTED ;;;;;;;;;;;;;;;;;
mov eax, 0xb5 ;; materialize vector constants from an immediate
vpbroadcastd zmm4, eax ; largest input that won't overflow/saturate
vpsrlw zmm5, zmm4, 5 ; fives = 0xb5 >> 5 = 5
;vpcmpeqd xmm3, xmm3 ; all-ones: result on saturation
vpternlogd zmm3,zmm3,zmm3, 0xff ; alternative for ZMM regs, where there's no compare-into-vector
.loop:
; alignment recommended for 512-bit vectors, but `u` lets it still work (slower) on unaligned.
vmovdqu64 zmm0, [rsi]
;; if x>0xb5 (unsigned), saturate output to 0xffff
vpcmpuw k1, zmm0, zmm4, 2 ; k1 = x <= 0xb5. 2 = LE predicate
; k1 set for elements that WON'T overflow
vpmullw zmm0, zmm0 ; x*x
vpaddusw zmm0, zmm0 ; *= 2 ; saturating add or not doesn't matter here
vmovdqa64 zmm1, zmm3 ; 0xffff
vpaddusw zmm1k1, zmm0, zmm5 ; 0xffff or 2*x^2 + 5 merge masking
vmovdqu64 [rsi], zmm1
add rsi, 64
cmp rsi, rcx
jb .loop
（NASM 允许 vpmullw a, b 作为 vpaddusw a, a, b 的快捷方式，当您不想像 imul eax, 123 那样利用非破坏性目标 3 操作数编码时。）
应用饱和度的早期想法是使用vpblendmw 根据掩码在5 和0xffff 的向量之间进行选择。
vpcmpuw k1, xmm4, xmm0, 1 ; k1 = 0xb5<x = x>0xb5. 1 = LT predicate numerically because NASM doesn't seem to accept vpcmpltuw the way it accepts vpcmpeqb
; k1 = 1 for elements that WILL overflow.
multiply and add as usual
...
vpblendmw xmm1k1, xmm5, xmm3 ; select (k1==0) ? 5 : 0xffff
vpaddusw xmm0, xmm1 ; += 5 or += 0xffff with unsigned saturation.
复制寄存器仍需要前端微指令，但不需要后端 ALU 微指令。所以（特别是对于 512 位寄存器，其中端口 1 为 SKX 上的向量微指令关闭），这个vpblendmb 想法比复制和合并掩码更糟糕。
除此之外，IACA 认为 vpblendmw xmm1k1, xmm3, xmm5 对 XMM1 有输出依赖性，即使它实际上是只写的。 （我通过将其中的 8 个放在一个循环中进行测试，有/没有一个 dep-break vpxor）。掩码混合指令是一种特殊情况：对于未设置的掩码位意味着它从 src1 复制（或零掩码为零），而对于设置掩码位它从 src2 复制。
但是机器编码使用合并掩码，因此硬件可能会将其视为具有合并掩码的任何其他 ALU 操作。 （其中输出向量是第三个输入依赖项vpaddw xmm1k1, xmm2, xmm3：如果掩码有任何零，则 XMM1 中的结果将是该元素的输入值。）
这可能不是问题：根据 IACA，SKX 可以每 2.24 个周期运行一次迭代（在前端有瓶颈），因此只要通过 XMM1 的循环承载的 dep 链就不是问题它只有 1 个周期延迟。 （如果展开以减少循环开销/前端瓶颈，您应该为混合目标使用单独的向量来解耦迭代，即使您无法在每次迭代接近 1 个周期时将其降低。）
（并且使用复制 + 合并掩码到0xffff 的向量的版本也以该吞吐量运行，即使对于 ZMM 向量也是如此。但 IACA 认为使用 ZMM 的 vpblendmb 版本会更慢，即使它说两个瓶颈都在前端...）
【讨论】：
【参考方案2】：好的，所以我找到了解决方案。使用说明是pmullw。指令pmullw mm0, mm1将依次计算寄存器内4个字的乘积并将它们存储在mm0中。对于printf 问题，我只是在调用它之前推送了另一个寄存器rdx，现在它可以工作了。我认为这与提到的堆栈未对齐有关。如果有人可以更详细地向我解释这是如何工作的，那就太好了。
【讨论】：
ABI 要求 16 字节堆栈指针对齐。pmullw 可以乘法，但不会饱和。以上是关于使用 MMX 指令集计算 f(x)=2*(x^2)+5 饱和度，用于从二进制文件加载的 128 个大小为 2 字节的数字的主要内容，如果未能解决你的问题，请参考以下文章
MultiMedia eXtensions - MMX：第一套应用于英特尔 80x86 指令集的 SIMD 扩展