如何用x86-64 SIMD向量指令与%ymm寄存器实现多项式求值

我需要编写一个x86-64汇编函数，利用x86向量/SIMD指令和%ymm寄存器计算表达式：x*(x*(a*x + b) + c) + d。函数签名为：

void map_poly_double_vec(double* input, double* output, uint64_t length, double a, double b, double c, double d)

其中a、b、c、d为常量，input是存储double类型x值的数组指针，计算结果需存入output数组对应位置，且数组长度可被4整除。我已经用%xmm寄存器实现了单元素版本的函数及点积函数，但不知道如何用向量指令完成批量计算，希望得到帮助。

现有%xmm寄存器实现代码 #

.section .note.GNU-stack,""
    .global map_poly_double, dot_double_vec
    .text

dot_double_vec:
    vxorpd %xmm2, %xmm2, %xmm2
    vxorpd %ymm0, %ymm0, %ymm0
    vxorpd %ymm1, %ymm1, %ymm1
    mov $0, %rcx

ddv_loop:
    cmp %rdx, %rcx
    jae ddv_end

    vmovupd (%rdi, %rcx, 8), %ymm1
    vmulpd (%rsi, %rcx, 8), %ymm1, %ymm0

    vextractf128 $0x1, %ymm0, %xmm1
    vaddpd %xmm1, %xmm0, %xmm0
    vshufpd $0b01, %xmm0, %xmm0, %xmm1
    vaddsd %xmm1, %xmm0, %xmm0

    vaddsd %xmm0, %xmm2, %xmm2

    add $4, %rcx
    jmp ddv_loop

ddv_end:
    movq %xmm2, %xmm0
    ret

#######

map_poly_double:
    mov $0, %rcx

map_poly_double_loop:
    cmp %rcx, %rdx
    je map_poly_double_end

    pxor %xmm5, %xmm5

    movq (%rdi, %rcx, 8), %xmm4
    movq %xmm0, %xmm5
    mulsd %xmm4, %xmm5
    addsd %xmm1, %xmm5

    mulsd %xmm4, %xmm5
    addsd %xmm2, %xmm5

    mulsd %xmm4, %xmm5
    addsd %xmm3, %xmm5

    movq %xmm5, (%rsi, %rcx, 8)

    add $1, %rcx
    jmp map_poly_double_loop

map_poly_double_end:
    ret

用%ymm寄存器实现的向量版本函数 #

核心思路是利用AVX指令一次处理4个double元素（%ymm寄存器可容纳4个64位浮点数），通过批量乘法、加法指令完成多项式计算，循环步长设为4，充分利用SIMD的并行性。

.section .note.GNU-stack,""
    .global map_poly_double_vec
    .text

map_poly_double_vec:
    # 函数参数：
    # rdi = input指针, rsi = output指针, rdx = length
    # xmm0 = a, xmm1 = b, xmm2 = c, xmm3 = d

    # 将标量常量广播到整个ymm寄存器，得到4个相同的a/b/c/d值
    vbroadcastsd %xmm0, %ymm0   # ymm0 = [a, a, a, a]
    vbroadcastsd %xmm1, %ymm1   # ymm1 = [b, b, b, b]
    vbroadcastsd %xmm2, %ymm2   # ymm2 = [c, c, c, c]
    vbroadcastsd %xmm3, %ymm3   # ymm3 = [d, d, d, d]

    mov $0, %rcx                # 初始化循环计数器

mpdv_loop:
    cmp %rdx, %rcx
    jae mpdv_end                # 计数器达到length时退出循环

    # 加载4个x值到ymm4
    vmovupd (%rdi, %rcx, 8), %ymm4  # ymm4 = [x3, x2, x1, x0] (内存顺序)

    # 计算 a*x + b：ymm5 = ymm0 * ymm4 + ymm1
    vmulpd %ymm4, %ymm0, %ymm5
    vaddpd %ymm1, %ymm5, %ymm5

    # 计算 x*(a*x + b) + c：ymm5 = ymm4 * ymm5 + ymm2
    vmulpd %ymm4, %ymm5, %ymm5
    vaddpd %ymm2, %ymm5, %ymm5

    # 计算 x*(x*(a*x + b) + c) + d：ymm5 = ymm4 * ymm5 + ymm3
    vmulpd %ymm4, %ymm5, %ymm5
    vaddpd %ymm3, %ymm5, %ymm5

    # 将结果写入output数组
    vmovupd %ymm5, (%rsi, %rcx, 8)

    add $4, %rcx                # 步长为4，处理下一组元素
    jmp mpdv_loop

mpdv_end:
    vzeroupper                  # 清理ymm寄存器，避免AVX/SSE状态冲突
    ret

关键说明 #

• 广播常量：用vbroadcastsd把单个double常量复制到整个ymm寄存器，这样就能和批量加载的x值做向量运算。
• 向量运算指令：vmulpd和vaddpd分别是批量浮点数乘法和加法，一次处理4个元素，比单元素的mulsd/addsd效率高4倍。
• 内存访问：vmovupd用于非对齐内存访问（如果数组是对齐的，也可以用vmovapd，性能略好）。
• vzeroupper：函数结尾执行该指令，把ymm寄存器的高位清零，避免后续SSE指令出现性能问题或异常。

内容的提问来源于stack exchange，提问作者Gyro