使用从二进制文件加载的 128 个大小为 2 字节的数字的 MMX 指令集计算 f(x)=2*(x^2)+5 饱和度
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
我遇到这个问题，我需要使用 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;
}
好的，所以我找到了解决方案。使用说明为pmullw。指令pmullw mm0, mm1会依次计算寄存器中4个字的乘积，存入mm0。对于 printf 问题，我只是在调用它之前按下另一个寄存器 rdx，现在它可以工作了。我认为它与所提到的堆栈未对齐有关。如果有人能更详细地向我解释它是如何工作的，那就太好了。
如果您关心性能，请避免在现代 CPU 上使用 loop 指令。 Why is the loop instruction slow? Couldn't Intel have implemented it efficiently?。同样使用 SSE2 而不是 MMX；您的数组大小是 16 和 8 的倍数，并且您使用的是保证具有 SSE2 的 x86-64。 MMX 对此完全没有意义，除非您还为 Pentium III / Athlon-XP 和更早版本制作 32 位版本。
（我回答中的所有代码都将使用 8 字节 MMX 寄存器而不是 16 字节 XMM 寄存器，因为我使用的所有指令都有 MMX 版本。根据 appendix B to the NASM manual，[ =16=, pxor, pcmpgtw, 和 paddusw 都在原始的 P5 Pentium MMX 中可用。英特尔手册中列出的一些指令为 "MMX"（如 pmulhuw 和 pshufw) 只是后来添加的，就像 Pentium II 一样，或者在 Pentium III 中与 SSE 一起添加，但对于此处有用的任何指令都不是这种情况。)
请参阅 https://whosebug.com/tags/x86/info 了解性能/优化指南，以及解释调用函数所需的 16 字节堆栈对齐的 ABI/调用约定链接。
mov rax, 0 / mov ax, [r8] 真的很傻。像正常人一样使用 movzx eax, word [r8]。您也不需要 jmp 到下一个源代码行，例如 jmp .square_loop / .square_loop: 如果没有分支指令，执行总是自行进入下一行。
x86 SIMD 没有饱和乘法，只有饱和 signed/unsigned 加法和饱和填充到较窄的元素 。 (MMX/SSE2 paddsw / paddusw)。既然你用 %d 打印，也许你想要有符号的饱和度？但这只是在解压缩为 32 位之后，您的公式将始终产生正结果，因此您可以使用无符号饱和度。我看到这就是您的代码使用的 paddusw.
此外，使用 3 个单独的循环 store/reload 您在公式的每个步骤之间的数据真的很糟糕 。您（几乎）总是希望 增加 计算强度（每个 memory/cache 带宽对您的数据完成的 ALU 工作量）。您也不需要乘法指令来将数字加倍：只需将其与自身相加即可。 padd* 运行s 在比 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 乘积对的水平和饱和。
否则，您可能必须解压为双字，然后 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+psrlw 移动 x。当需要饱和时（否则为 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 可以 运行 pxor 在比 padd* 更多的端口上，并且它不需要任何更少的非零向量常量，所以这在任何方面都不会更好。但它在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.
自行饱和就好了
如果我们用 pcmpgtw（而不是 >=）检查 (x>>1) > (0xb6>>1) 以捕获像 0xffff 这样的输入（其中带有 0xffff 的 pmullw 给我们 0x0001）。所以我们可以保存 1 个向量常量，但除此之外这不是更好。
pxor + pcmpgtw 至少和 psrlw xmm, 1 + pcmpgtw 一样便宜，除非我们针对 Intel P6 系列进行调整（Core2/Nehalem) 和 运行ning 进入寄存器读取端口停顿。但这不太可能：xmm0、xmm1 和 rsi 应该始终是热的（最近写入并因此从 ROB 读取，而不是永久寄存器文件）。我们只在循环的第一组 4 条指令中读取 2 个向量常量，然后再读取 1 个。
正如我在下面所说的，在许多 Intel CPU 上，psrlw 只能在与 pmullw 相同的端口上，在 vec-int shift+multiply 执行单元上。不过，这里可能不是吞吐量瓶颈。
但是 pcmp 和 padd 运行 在有限的端口上（在 Skylake 之前的英特尔上），而 pxor 可以 运行任何矢量 ALU 端口。 纯 padd/pcmp/pmul/psrl` uops 的混合将使一个矢量 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 相互竞争吞吐量。这是足够的总 uops，我们实际上不会在 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 zmm1{k1}, 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 的快捷方式，当您不想像它那样利用非破坏性目标 3 操作数编码时imul eax, 123.)
应用饱和度的早期想法是根据掩码在 5 和 0xffff 向量之间使用 vpblendmw 到 select。
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 xmm1{k1}, xmm5, xmm3 ; select (k1==0) ? 5 : 0xffff
vpaddusw xmm0, xmm1 ; += 5 or += 0xffff with unsigned saturation.
复制寄存器仍然需要前端微指令，但不需要后端 ALU 微指令。所以（特别是对于 512 位寄存器，端口 1 在 SKX 上为矢量微指令关闭），这个 vpblendmb 想法比复制和合并掩码更糟糕。
除此之外，IACA 认为 vpblendmw xmm1{k1}, xmm3, xmm5 对 XMM1 有输出依赖性，即使它实际上只写 。 （我通过将其中的 8 个放在一个循环中进行测试，with/without 是一个 dep-breaking vpxor）。掩码混合指令是一种特殊情况：对于未设置掩码位意味着它从 src1 复制（或零掩码为零），而对于设置掩码位它从 src2 复制。
但机器编码使用合并屏蔽，因此硬件可能会像对待任何其他具有合并屏蔽的 ALU 操作一样对待它。 （输出向量是第三个输入依赖项，vpaddw xmm1{k1}, xmm2, xmm3：如果掩码有任何零，则 XMM1 中的结果将是该元素的输入值。）
这可能不是问题：根据 IACA，SKX 可以 运行 每 2.24 个周期进行一次迭代（在前端出现瓶颈），因此通过 XMM1 的循环承载 dep 链不是'只要只有 1 个周期的延迟就不是问题。 （如果展开以减少循环开销/前端瓶颈，您应该为混合目标使用一个单独的向量来解耦迭代，即使您无法将其降低到每次迭代 1 个循环附近。）
（并且使用复制+合并屏蔽到 0xffff 向量中的版本也在该吞吐量下 运行s，即使对于 ZMM 向量也是如此。但是 IACA 认为 vpblendmb 版本会更慢ZMM，尽管它说前端都有瓶颈...)
我遇到这个问题，我需要使用 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;
}
好的，所以我找到了解决方案。使用说明为pmullw。指令pmullw mm0, mm1会依次计算寄存器中4个字的乘积，存入mm0。对于 printf 问题，我只是在调用它之前按下另一个寄存器 rdx，现在它可以工作了。我认为它与所提到的堆栈未对齐有关。如果有人能更详细地向我解释它是如何工作的，那就太好了。
如果您关心性能，请避免在现代 CPU 上使用 loop 指令。 Why is the loop instruction slow? Couldn't Intel have implemented it efficiently?。同样使用 SSE2 而不是 MMX；您的数组大小是 16 和 8 的倍数，并且您使用的是保证具有 SSE2 的 x86-64。 MMX 对此完全没有意义，除非您还为 Pentium III / Athlon-XP 和更早版本制作 32 位版本。
（我回答中的所有代码都将使用 8 字节 MMX 寄存器而不是 16 字节 XMM 寄存器，因为我使用的所有指令都有 MMX 版本。根据 appendix B to the NASM manual，[ =16=, pxor, pcmpgtw, 和 paddusw 都在原始的 P5 Pentium MMX 中可用。英特尔手册中列出的一些指令为 "MMX"（如 pmulhuw 和 pshufw) 只是后来添加的，就像 Pentium II 一样，或者在 Pentium III 中与 SSE 一起添加，但对于此处有用的任何指令都不是这种情况。)
请参阅 https://whosebug.com/tags/x86/info 了解性能/优化指南，以及解释调用函数所需的 16 字节堆栈对齐的 ABI/调用约定链接。
mov rax, 0 / mov ax, [r8] 真的很傻。像正常人一样使用 movzx eax, word [r8]。您也不需要 jmp 到下一个源代码行，例如 jmp .square_loop / .square_loop: 如果没有分支指令，执行总是自行进入下一行。
x86 SIMD 没有饱和乘法，只有饱和 signed/unsigned 加法和饱和填充到较窄的元素 。 (MMX/SSE2 paddsw / paddusw)。既然你用 %d 打印，也许你想要有符号的饱和度？但这只是在解压缩为 32 位之后，您的公式将始终产生正结果，因此您可以使用无符号饱和度。我看到这就是您的代码使用的 paddusw.
此外，使用 3 个单独的循环 store/reload 您在公式的每个步骤之间的数据真的很糟糕 。您（几乎）总是希望 增加 计算强度（每个 memory/cache 带宽对您的数据完成的 ALU 工作量）。您也不需要乘法指令来将数字加倍：只需将其与自身相加即可。 padd* 运行s 在比 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 乘积对的水平和饱和。
否则，您可能必须解压为双字，然后 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+psrlw 移动 x。当需要饱和时（否则为 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. 这样的小数字后将保持 "negative"
paddusw并检查== 0xffff：只有足够小的输入才会使最终结果不饱和。有些 CPU 可以 运行pxor在比padd*更多的端口上，并且它不需要任何更少的非零向量常量，所以这在任何方面都不会更好。但它在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.
如果我们用 pcmpgtw（而不是 >=）检查 (x>>1) > (0xb6>>1) 以捕获像 0xffff 这样的输入（其中带有 0xffff 的 pmullw 给我们 0x0001）。所以我们可以保存 1 个向量常量，但除此之外这不是更好。
pxor + pcmpgtw 至少和 psrlw xmm, 1 + pcmpgtw 一样便宜，除非我们针对 Intel P6 系列进行调整（Core2/Nehalem) 和 运行ning 进入寄存器读取端口停顿。但这不太可能：xmm0、xmm1 和 rsi 应该始终是热的（最近写入并因此从 ROB 读取，而不是永久寄存器文件）。我们只在循环的第一组 4 条指令中读取 2 个向量常量，然后再读取 1 个。
正如我在下面所说的，在许多 Intel CPU 上，psrlw 只能在与 pmullw 相同的端口上，在 vec-int shift+multiply 执行单元上。不过，这里可能不是吞吐量瓶颈。
但是 pcmp 和 padd 运行 在有限的端口上（在 Skylake 之前的英特尔上），而 pxor 可以 运行任何矢量 ALU 端口。 纯 padd/pcmp/pmul/psrl` uops 的混合将使一个矢量 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 相互竞争吞吐量。这是足够的总 uops，我们实际上不会在 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 zmm1{k1}, 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 的快捷方式，当您不想像它那样利用非破坏性目标 3 操作数编码时imul eax, 123.)
应用饱和度的早期想法是根据掩码在 5 和 0xffff 向量之间使用 vpblendmw 到 select。
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 xmm1{k1}, xmm5, xmm3 ; select (k1==0) ? 5 : 0xffff
vpaddusw xmm0, xmm1 ; += 5 or += 0xffff with unsigned saturation.
复制寄存器仍然需要前端微指令，但不需要后端 ALU 微指令。所以（特别是对于 512 位寄存器，端口 1 在 SKX 上为矢量微指令关闭），这个 vpblendmb 想法比复制和合并掩码更糟糕。
除此之外，IACA 认为 vpblendmw xmm1{k1}, xmm3, xmm5 对 XMM1 有输出依赖性，即使它实际上只写 。 （我通过将其中的 8 个放在一个循环中进行测试，with/without 是一个 dep-breaking vpxor）。掩码混合指令是一种特殊情况：对于未设置掩码位意味着它从 src1 复制（或零掩码为零），而对于设置掩码位它从 src2 复制。
但机器编码使用合并屏蔽，因此硬件可能会像对待任何其他具有合并屏蔽的 ALU 操作一样对待它。 （输出向量是第三个输入依赖项，vpaddw xmm1{k1}, xmm2, xmm3：如果掩码有任何零，则 XMM1 中的结果将是该元素的输入值。）
这可能不是问题：根据 IACA，SKX 可以 运行 每 2.24 个周期进行一次迭代（在前端出现瓶颈），因此通过 XMM1 的循环承载 dep 链不是'只要只有 1 个周期的延迟就不是问题。 （如果展开以减少循环开销/前端瓶颈，您应该为混合目标使用一个单独的向量来解耦迭代，即使您无法将其降低到每次迭代 1 个循环附近。）
（并且使用复制+合并屏蔽到 0xffff 向量中的版本也在该吞吐量下 运行s，即使对于 ZMM 向量也是如此。但是 IACA 认为 vpblendmb 版本会更慢ZMM，尽管它说前端都有瓶颈...)