; Copyright (C) 2026 Kiyotsugu Arai
; SPDX-License-Identifier: LGPL-3.0-or-later
;
; mpn_x64_mul.asm — mul_basecase small-size specialization (n×n, n=1..7 and n=8)
;
; Functions:
;   mpn_mul_small_asm(rp, ap, bp, n) — n×n fully expanded multiplication (n=1..7)
;   mpn_mul_8x8_asm(rp, ap, bp)      — 8×8 fully expanded, MULX + ADCX/ADOX dual chain
;
; Purpose: eliminate the push/pop + loop overhead of mul_basecase
;   n=2: 4 MULX (no push/pop)
;   n=3: 9 MULX (2 push/pop)
;   n=4: 16 MULX (4 push/pop)
;   n=5: 25 MULX (cyclic 6-reg, ADCX/ADOX dual chain, 5 push/pop)
;   n=6: 36 MULX (cyclic 7-reg, dual chain, 6 push/pop)
;   n=7: 49 MULX (cyclic 8-reg, dual chain, 7 push/pop)
;   n=8: 64 MULX (cyclic 9-reg, dual chain, 8 push/pop)
;
; BMI2 + ADX required (MULX/ADCX/ADOX)
;
; Windows x64 calling convention:
;   rcx = rp, rdx = ap, r8 = bp, r9 = n (1..4)
;   return value: none
;   non-volatile: rbx, rbp, rdi, rsi, r12-r15
;   volatile: rax, rcx, rdx, r8, r9, r10, r11

.code

; ============ Macros for n=5,6,7 dual-chain row body ============
; ROWN_DUAL expands Row j (j ≥ 1) of the cyclic accumulator pattern:
;   1) zero out p_n (new top = reg stored in the previous row) via xor + clear CF/OF
;   2) MULX × n + ADCX/ADOX dual chain
;   3) absorb remaining CF/OF into p_n
;   4) store r[j] = p0
; Precondition: rbx=bp, r9=ap, rcx=rp, p0..p_{n-1} are p[j..j+n-1] (accumulated values up to the previous row),
;        p_n is a discardable reg (stored in the previous row, value no longer needed)
;
; ROW5_DUAL: 5 MULX
ROW5_DUAL MACRO j_off:REQ, p0:REQ, p1:REQ, p2:REQ, p3:REQ, p4:REQ, p5:REQ
    xor     p5, p5
    mov     rdx, [rbx + 8*j_off]
    mulx    r8, rax, [r9]
    adcx    p0, rax
    adox    p1, r8
    mulx    r8, rax, [r9+8]
    adcx    p1, rax
    adox    p2, r8
    mulx    r8, rax, [r9+16]
    adcx    p2, rax
    adox    p3, r8
    mulx    r8, rax, [r9+24]
    adcx    p3, rax
    adox    p4, r8
    mulx    r8, rax, [r9+32]
    adcx    p4, rax
    adox    p5, r8
    mov     rax, 0
    adcx    p5, rax
    adox    p5, rax
    mov     [rcx + 8*j_off], p0
ENDM

; ROW6_DUAL: 6 MULX
ROW6_DUAL MACRO j_off:REQ, p0:REQ, p1:REQ, p2:REQ, p3:REQ, p4:REQ, p5:REQ, p6:REQ
    xor     p6, p6
    mov     rdx, [rbx + 8*j_off]
    mulx    r8, rax, [r9]
    adcx    p0, rax
    adox    p1, r8
    mulx    r8, rax, [r9+8]
    adcx    p1, rax
    adox    p2, r8
    mulx    r8, rax, [r9+16]
    adcx    p2, rax
    adox    p3, r8
    mulx    r8, rax, [r9+24]
    adcx    p3, rax
    adox    p4, r8
    mulx    r8, rax, [r9+32]
    adcx    p4, rax
    adox    p5, r8
    mulx    r8, rax, [r9+40]
    adcx    p5, rax
    adox    p6, r8
    mov     rax, 0
    adcx    p6, rax
    adox    p6, rax
    mov     [rcx + 8*j_off], p0
ENDM

; ROW7_DUAL: 7 MULX
ROW7_DUAL MACRO j_off:REQ, p0:REQ, p1:REQ, p2:REQ, p3:REQ, p4:REQ, p5:REQ, p6:REQ, p7:REQ
    xor     p7, p7
    mov     rdx, [rbx + 8*j_off]
    mulx    r8, rax, [r9]
    adcx    p0, rax
    adox    p1, r8
    mulx    r8, rax, [r9+8]
    adcx    p1, rax
    adox    p2, r8
    mulx    r8, rax, [r9+16]
    adcx    p2, rax
    adox    p3, r8
    mulx    r8, rax, [r9+24]
    adcx    p3, rax
    adox    p4, r8
    mulx    r8, rax, [r9+32]
    adcx    p4, rax
    adox    p5, r8
    mulx    r8, rax, [r9+40]
    adcx    p5, rax
    adox    p6, r8
    mulx    r8, rax, [r9+48]
    adcx    p6, rax
    adox    p7, r8
    mov     rax, 0
    adcx    p7, rax
    adox    p7, rax
    mov     [rcx + 8*j_off], p0
ENDM

; =====================================================================
; void mpn_mul_small_asm(uint64_t* rp, const uint64_t* ap,
;                         const uint64_t* bp, size_t n)
;
; rp[0..2n-1] = ap[0..n-1] × bp[0..n-1]
; precondition: 1 <= n <= 7, CPU supporting BMI2 + ADX
;   n=1..4: simple ADD/ADC unrolling up to 4×4
;   n=5..7: cyclic accumulator (N=n+1) + MULX/ADCX/ADOX dual chain
;
; Dispatch per n to the fully unrolled version via the jump table.
; =====================================================================
mpn_mul_small_asm PROC
    lea     r10, [mul_s_jt]
    jmp     qword ptr [r10 + r9*8 - 8]

mul_s_jt:
    dq      mul_s_1x1
    dq      mul_s_2x2
    dq      mul_s_3x3
    dq      mul_s_4x4
    dq      mul_s_5x5
    dq      mul_s_6x6
    dq      mul_s_7x7

; --- 1×1: 1 MULX ---
    ALIGN   16
mul_s_1x1:
    mov     r9, rdx             ; r9 = ap
    mov     rdx, [r8]           ; rdx = b[0]
    mulx    r10, rax, [r9]      ; r10:rax = a[0]*b[0]
    mov     [rcx], rax
    mov     [rcx+8], r10
    ret

; --- 2×2: 4 MULX, no push/pop ---
;
; r[0..3] = a[0..1] × b[0..1]
;
; Registers:
;   r9=ap, rdx=MULX multiplier, rcx=rp
;   r10,r11=accumulator, rax,r8=temp
    ALIGN   16
mul_s_2x2:
    mov     r9, rdx             ; r9 = ap

    ; Row 0: a × b[0]
    mov     rdx, [r8]           ; rdx = b[0]
    mulx    r11, r10, [r9]      ; r11:r10 = a[0]*b[0]
    mov     [rcx], r10          ; r[0]
    mulx    rax, r10, [r9+8]    ; rax:r10 = a[1]*b[0]
    add     r10, r11            ; partial r[1]
    adc     rax, 0              ; partial r[2]

    ; Row 1: a × b[1]
    mov     rdx, [r8+8]         ; rdx = b[1]
    mulx    r11, r8, [r9]       ; r11:r8 = a[0]*b[1]
    add     r10, r8             ; r[1]
    mov     [rcx+8], r10
    mulx    r10, r8, [r9+8]     ; r10:r8 = a[1]*b[1]
    adc     rax, r11            ; partial r[2] += h(a0*b1)
    adc     r10, 0              ; r10 = h(a1*b1) + carry
    add     rax, r8             ; r[2] += l(a1*b1)
    mov     [rcx+16], rax
    adc     r10, 0              ; r[3]
    mov     [rcx+24], r10
    ret

; --- 3×3: 9 MULX, 2 push/pop ---
;
; r[0..5] = a[0..2] × b[0..2]
;
; Registers:
;   r9=ap, rbx=bp (callee-saved), rdx=MULX multiplier, rcx=rp
;   r10,r11,r8=accumulator, rax,rdi=temp
    ALIGN   16
mul_s_3x3:
    push    rbx
    push    rdi
    mov     r9, rdx             ; r9 = ap
    mov     rbx, r8             ; rbx = bp

    ; Row 0: a × b[0]
    mov     rdx, [rbx]          ; rdx = b[0]
    mulx    r11, r10, [r9]      ; r11:r10 = a[0]*b[0]
    mov     [rcx], r10          ; r[0]
    mulx    rdi, r10, [r9+8]    ; rdi:r10 = a[1]*b[0]
    add     r10, r11            ; partial r[1]
    mulx    r8, r11, [r9+16]    ; r8:r11 = a[2]*b[0]
    adc     r11, rdi            ; partial r[2]
    adc     r8, 0               ; partial r[3]

    ; Row 1: a × b[1]
    mov     rdx, [rbx+8]        ; rdx = b[1]
    mulx    rdi, rax, [r9]      ; rdi:rax = a[0]*b[1]
    add     r10, rax            ; r[1]
    mov     [rcx+8], r10

    mulx    rax, r10, [r9+8]    ; rax:r10 = a[1]*b[1]
    adc     r11, rdi
    adc     r8, 0
    add     r11, r10            ; partial r[2] += l(a1*b1)

    mulx    rdi, r10, [r9+16]   ; rdi:r10 = a[2]*b[1]
    adc     r8, rax             ; partial r[3] += h(a1*b1)
    adc     rdi, 0              ; partial r[4]
    add     r8, r10             ; partial r[3] += l(a2*b1)
    adc     rdi, 0

    ; Row 2: a × b[2]
    mov     rdx, [rbx+16]       ; rdx = b[2]
    mulx    rax, r10, [r9]      ; rax:r10 = a[0]*b[2]
    add     r11, r10            ; r[2]
    mov     [rcx+16], r11

    mulx    r10, r11, [r9+8]    ; r10:r11 = a[1]*b[2]
    adc     r8, rax
    adc     rdi, 0
    add     r8, r11             ; r[3]
    mov     [rcx+24], r8

    mulx    rax, r11, [r9+16]   ; rax:r11 = a[2]*b[2]
    adc     rdi, r10            ; partial r[4] += h(a1*b2)
    adc     rax, 0              ; partial r[5]
    add     rdi, r11            ; r[4] += l(a2*b2)
    mov     [rcx+32], rdi
    adc     rax, 0              ; r[5]
    mov     [rcx+40], rax

    pop     rdi
    pop     rbx
    ret

; --- 4×4: 16 MULX, 4 push/pop ---
;
; r[0..7] = a[0..3] × b[0..3]
;
; Registers:
;   r9=ap, rbx=bp, rdx=MULX, rcx=rp
;   rdi,rsi=callee-saved temp
;   r10,r11,r8,rax=accumulator
    ALIGN   16
mul_s_4x4:
    push    rbx
    push    rdi
    push    rsi
    push    rbp
    mov     r9, rdx             ; r9 = ap
    mov     rbx, r8             ; rbx = bp

    ; Row 0: a × b[0]
    mov     rdx, [rbx]          ; rdx = b[0]
    mulx    r11, r10, [r9]      ; a[0]*b[0]
    mov     [rcx], r10          ; r[0]
    mulx    rdi, r10, [r9+8]    ; a[1]*b[0]
    add     r10, r11            ; partial r[1]
    mulx    rsi, r11, [r9+16]   ; a[2]*b[0]
    adc     r11, rdi            ; partial r[2]
    mulx    r8, rdi, [r9+24]    ; a[3]*b[0]
    adc     rdi, rsi            ; partial r[3]
    adc     r8, 0               ; partial r[4]
    ; r10=r[1]p, r11=r[2]p, rdi=r[3]p, r8=r[4]p

    ; Row 1: a × b[1]
    mov     rdx, [rbx+8]
    mulx    rsi, rax, [r9]      ; a[0]*b[1]
    add     r10, rax            ; r[1]
    mov     [rcx+8], r10

    mulx    rax, r10, [r9+8]    ; a[1]*b[1]
    adc     r11, rsi
    adc     rdi, 0
    adc     r8, 0
    add     r11, r10

    mulx    r10, rsi, [r9+16]   ; a[2]*b[1]
    adc     rdi, rax
    adc     r8, 0
    add     rdi, rsi

    mulx    rsi, rax, [r9+24]   ; a[3]*b[1]
    adc     r8, r10
    adc     rsi, 0
    add     r8, rax
    adc     rsi, 0
    ; r11=r[2]p, rdi=r[3]p, r8=r[4]p, rsi=r[5]p

    ; Row 2: a × b[2]
    mov     rdx, [rbx+16]
    mulx    rax, r10, [r9]      ; a[0]*b[2]
    add     r11, r10            ; r[2]
    mov     [rcx+16], r11

    mulx    r10, r11, [r9+8]    ; a[1]*b[2]
    adc     rdi, rax
    adc     r8, 0
    adc     rsi, 0
    add     rdi, r11

    mulx    r11, rax, [r9+16]   ; a[2]*b[2]
    adc     r8, r10
    adc     rsi, 0
    add     r8, rax

    mulx    rbp, rax, [r9+24]   ; a[3]*b[2]
    adc     rsi, r11
    adc     rbp, 0
    add     rsi, rax
    adc     rbp, 0
    ; rdi=r[3]p, r8=r[4]p, rsi=r[5]p, rbp=r[6]p

    ; Row 3: a × b[3]
    mov     rdx, [rbx+24]
    mulx    rax, r10, [r9]      ; a[0]*b[3]
    add     rdi, r10            ; r[3]
    mov     [rcx+24], rdi

    mulx    r10, r11, [r9+8]    ; a[1]*b[3]
    adc     r8, rax
    adc     rsi, 0
    adc     rbp, 0
    add     r8, r11             ; r[4]
    mov     [rcx+32], r8

    mulx    r11, rax, [r9+16]   ; a[2]*b[3]
    adc     rsi, r10
    adc     rbp, 0
    add     rsi, rax            ; r[5]
    mov     [rcx+40], rsi

    mulx    rax, r10, [r9+24]   ; a[3]*b[3]
    adc     rbp, r11
    adc     rax, 0
    add     rbp, r10            ; r[6]
    mov     [rcx+48], rbp
    adc     rax, 0              ; r[7]
    mov     [rcx+56], rax

    pop     rbp
    pop     rsi
    pop     rdi
    pop     rbx
    ret

; --- 5×5: 25 MULX, cyclic 6-reg accumulator, ADCX/ADOX dual chain ---
;
; r[0..9] = a[0..4] × b[0..4]
; accumulator reg (cyclic, N=6): [r10, r11, rdi, rbp, r12, r13]
; partial[k] = regs[k mod 6]; Row 0 initializes p[1..5] = r11..r13
; Row 1..4: dual chain. At the start of each row, zero out the new top via xor
    ALIGN   16
mul_s_5x5:
    push    rbx
    push    rdi
    push    rbp
    push    r12
    push    r13
    mov     r9, rdx                 ; r9 = ap
    mov     rbx, r8                 ; rbx = bp

    ; ============ Row 0: a × b[0] (MULX-direct initialization) ============
    ; p[1..5] = r11, rdi, rbp, r12, r13
    mov     rdx, [rbx]              ; b[0]
    mulx    r11, rax, [r9]
    mov     [rcx], rax              ; r[0]
    mulx    rdi, rax, [r9+8]
    add     r11, rax
    mulx    rbp, rax, [r9+16]
    adc     rdi, rax
    mulx    r12, rax, [r9+24]
    adc     rbp, rax
    mulx    r13, rax, [r9+32]
    adc     r12, rax
    adc     r13, 0                  ; p[5]

    ; Rows 1..4 (cyclic dual chain)
    ; p_n (= new top) is zeroed out by the xor inside the macro
    ROW5_DUAL 1, r11, rdi, rbp, r12, r13, r10
    ROW5_DUAL 2, rdi, rbp, r12, r13, r10, r11
    ROW5_DUAL 3, rbp, r12, r13, r10, r11, rdi
    ROW5_DUAL 4, r12, r13, r10, r11, rdi, rbp

    ; ============ Final stores: r[5..9] ============
    ; After row 4: p[5..9] = [r13, r10, r11, rdi, rbp]
    mov     [rcx+40], r13           ; r[5]
    mov     [rcx+48], r10           ; r[6]
    mov     [rcx+56], r11           ; r[7]
    mov     [rcx+64], rdi           ; r[8]
    mov     [rcx+72], rbp           ; r[9]

    pop     r13
    pop     r12
    pop     rbp
    pop     rdi
    pop     rbx
    ret

; --- 6×6: 36 MULX, cyclic 7-reg accumulator, ADCX/ADOX dual chain ---
;
; r[0..11] = a[0..5] × b[0..5]
; accumulator reg (cyclic, N=7): [r10, r11, rdi, rbp, r12, r13, r14]
    ALIGN   16
mul_s_6x6:
    push    rbx
    push    rdi
    push    rbp
    push    r12
    push    r13
    push    r14
    mov     r9, rdx
    mov     rbx, r8

    ; ============ Row 0: p[1..6] = r11, rdi, rbp, r12, r13, r14 ============
    mov     rdx, [rbx]
    mulx    r11, rax, [r9]
    mov     [rcx], rax
    mulx    rdi, rax, [r9+8]
    add     r11, rax
    mulx    rbp, rax, [r9+16]
    adc     rdi, rax
    mulx    r12, rax, [r9+24]
    adc     rbp, rax
    mulx    r13, rax, [r9+32]
    adc     r12, rax
    mulx    r14, rax, [r9+40]
    adc     r13, rax
    adc     r14, 0                  ; p[6]

    ; Rows 1..5
    ROW6_DUAL 1, r11, rdi, rbp, r12, r13, r14, r10
    ROW6_DUAL 2, rdi, rbp, r12, r13, r14, r10, r11
    ROW6_DUAL 3, rbp, r12, r13, r14, r10, r11, rdi
    ROW6_DUAL 4, r12, r13, r14, r10, r11, rdi, rbp
    ROW6_DUAL 5, r13, r14, r10, r11, rdi, rbp, r12

    ; ============ Final stores: r[6..11] ============
    ; After row 5: p[6..11] = [r14, r10, r11, rdi, rbp, r12]
    mov     [rcx+48], r14           ; r[6]
    mov     [rcx+56], r10           ; r[7]
    mov     [rcx+64], r11           ; r[8]
    mov     [rcx+72], rdi           ; r[9]
    mov     [rcx+80], rbp           ; r[10]
    mov     [rcx+88], r12           ; r[11]

    pop     r14
    pop     r13
    pop     r12
    pop     rbp
    pop     rdi
    pop     rbx
    ret

; --- 7×7: 49 MULX, cyclic 8-reg accumulator, ADCX/ADOX dual chain ---
;
; r[0..13] = a[0..6] × b[0..6]
; accumulator reg (cyclic, N=8): [r10, r11, rdi, rbp, r12, r13, r14, r15]
    ALIGN   16
mul_s_7x7:
    push    rbx
    push    rdi
    push    rbp
    push    r12
    push    r13
    push    r14
    push    r15
    mov     r9, rdx
    mov     rbx, r8

    ; ============ Row 0: p[1..7] = r11, rdi, rbp, r12, r13, r14, r15 ============
    mov     rdx, [rbx]
    mulx    r11, rax, [r9]
    mov     [rcx], rax
    mulx    rdi, rax, [r9+8]
    add     r11, rax
    mulx    rbp, rax, [r9+16]
    adc     rdi, rax
    mulx    r12, rax, [r9+24]
    adc     rbp, rax
    mulx    r13, rax, [r9+32]
    adc     r12, rax
    mulx    r14, rax, [r9+40]
    adc     r13, rax
    mulx    r15, rax, [r9+48]
    adc     r14, rax
    adc     r15, 0                  ; p[7]

    ; Rows 1..6
    ROW7_DUAL 1, r11, rdi, rbp, r12, r13, r14, r15, r10
    ROW7_DUAL 2, rdi, rbp, r12, r13, r14, r15, r10, r11
    ROW7_DUAL 3, rbp, r12, r13, r14, r15, r10, r11, rdi
    ROW7_DUAL 4, r12, r13, r14, r15, r10, r11, rdi, rbp
    ROW7_DUAL 5, r13, r14, r15, r10, r11, rdi, rbp, r12
    ROW7_DUAL 6, r14, r15, r10, r11, rdi, rbp, r12, r13

    ; ============ Final stores: r[7..13] ============
    ; After row 6: p[7..13] = [r15, r10, r11, rdi, rbp, r12, r13]
    mov     [rcx+56], r15           ; r[7]
    mov     [rcx+64], r10           ; r[8]
    mov     [rcx+72], r11           ; r[9]
    mov     [rcx+80], rdi           ; r[10]
    mov     [rcx+88], rbp           ; r[11]
    mov     [rcx+96], r12           ; r[12]
    mov     [rcx+104], r13          ; r[13]

    pop     r15
    pop     r14
    pop     r13
    pop     r12
    pop     rbp
    pop     rdi
    pop     rbx
    ret

mpn_mul_small_asm ENDP

; =====================================================================
; void mpn_mul_8x8_asm(uint64_t* rp, const uint64_t* ap,
;                       const uint64_t* bp)
;
; rp[0..15] = ap[0..7] × bp[0..7]
; precondition: CPU supporting BMI2 + ADX (MULX/ADCX/ADOX)
;
; Design:
;   register-resident 9-slot cyclic accumulator [r10, r11, rdi, rbp, r12,
;   r13, r14, r15, rsi], partial[k] = regs[k mod 9].
;   Row 0 is initialized MULX-direct; Rows 1..7 accumulate with a MULX +
;   ADCX/ADOX dual carry chain. At the start of each row, zero out the new top
;   (= reg stored in the previous row) via XOR, simultaneously clearing CF/OF.
;
; Windows x64: rcx=rp, rdx=ap, r8=bp
; Registers:
;   rcx=rp, rdx=MULX multiplier, r9=ap, rbx=bp
;   rax/r8=MULX output (lo/hi)
;   accumulator (cyclic, 9 slots): r10, r11, rdi, rbp, r12, r13, r14, r15, rsi
; Push/pop: rbx, rsi, rdi, rbp, r12-r15 = 8 (maintains 16-byte stack alignment)
; =====================================================================
mpn_mul_8x8_asm PROC
    push    rbx
    push    rsi
    push    rdi
    push    rbp
    push    r12
    push    r13
    push    r14
    push    r15
    mov     r9, rdx                 ; r9 = ap
    mov     rbx, r8                 ; rbx = bp

    ; ============ Row 0: a × b[0] (initialized MULX-direct) ============
    ; partial[k] = regs9[k mod 9], regs9 = [r10,r11,rdi,rbp,r12,r13,r14,r15,rsi]
    ; Row 0 initializes p[1..8] = r11, rdi, rbp, r12, r13, r14, r15, rsi
    mov     rdx, [rbx]              ; b[0]
    mulx    r11, rax, [r9]          ; rax=lo, r11=hi → p[1]
    mov     [rcx], rax              ; r[0]
    mulx    rdi, rax, [r9+8]        ; rdi = hi → p[2] init
    add     r11, rax                ; p[1] += lo
    mulx    rbp, rax, [r9+16]
    adc     rdi, rax
    mulx    r12, rax, [r9+24]
    adc     rbp, rax
    mulx    r13, rax, [r9+32]
    adc     r12, rax
    mulx    r14, rax, [r9+40]
    adc     r13, rax
    mulx    r15, rax, [r9+48]
    adc     r14, rax
    mulx    rsi, rax, [r9+56]
    adc     r15, rax
    adc     rsi, 0                  ; p[8]

    ; r10 = regs9[0] is unused in Row 0 (becomes Row 1's new top p[9])
    xor     r10, r10                ; clear (= 0); also clears CF/OF

    ; ============ Row 1: a × b[1] ============
    ; p[1..9] = r11, rdi, rbp, r12, r13, r14, r15, rsi, r10
    ; r10 (new top) is already cleared (the xor above)
    mov     rdx, [rbx+8]            ; b[1]
    mulx    r8, rax, [r9]
    adcx    r11, rax                ; p[1] += lo
    adox    rdi, r8                 ; p[2] += hi
    mulx    r8, rax, [r9+8]
    adcx    rdi, rax
    adox    rbp, r8
    mulx    r8, rax, [r9+16]
    adcx    rbp, rax
    adox    r12, r8
    mulx    r8, rax, [r9+24]
    adcx    r12, rax
    adox    r13, r8
    mulx    r8, rax, [r9+32]
    adcx    r13, rax
    adox    r14, r8
    mulx    r8, rax, [r9+40]
    adcx    r14, rax
    adox    r15, r8
    mulx    r8, rax, [r9+48]
    adcx    r15, rax
    adox    rsi, r8
    mulx    r8, rax, [r9+56]
    adcx    rsi, rax
    adox    r10, r8                 ; p[9] (new top)
    mov     rax, 0
    adcx    r10, rax                ; absorb remaining CF
    adox    r10, rax                ; absorb remaining OF
    mov     [rcx+8], r11            ; r[1]

    ; ============ Row 2: a × b[2] ============
    ; p[2..10] = rdi, rbp, r12, r13, r14, r15, rsi, r10, r11
    xor     r11, r11
    mov     rdx, [rbx+16]
    mulx    r8, rax, [r9]
    adcx    rdi, rax
    adox    rbp, r8
    mulx    r8, rax, [r9+8]
    adcx    rbp, rax
    adox    r12, r8
    mulx    r8, rax, [r9+16]
    adcx    r12, rax
    adox    r13, r8
    mulx    r8, rax, [r9+24]
    adcx    r13, rax
    adox    r14, r8
    mulx    r8, rax, [r9+32]
    adcx    r14, rax
    adox    r15, r8
    mulx    r8, rax, [r9+40]
    adcx    r15, rax
    adox    rsi, r8
    mulx    r8, rax, [r9+48]
    adcx    rsi, rax
    adox    r10, r8
    mulx    r8, rax, [r9+56]
    adcx    r10, rax
    adox    r11, r8
    mov     rax, 0
    adcx    r11, rax
    adox    r11, rax
    mov     [rcx+16], rdi           ; r[2]

    ; ============ Row 3: a × b[3] ============
    ; p[3..11] = rbp, r12, r13, r14, r15, rsi, r10, r11, rdi
    xor     rdi, rdi
    mov     rdx, [rbx+24]
    mulx    r8, rax, [r9]
    adcx    rbp, rax
    adox    r12, r8
    mulx    r8, rax, [r9+8]
    adcx    r12, rax
    adox    r13, r8
    mulx    r8, rax, [r9+16]
    adcx    r13, rax
    adox    r14, r8
    mulx    r8, rax, [r9+24]
    adcx    r14, rax
    adox    r15, r8
    mulx    r8, rax, [r9+32]
    adcx    r15, rax
    adox    rsi, r8
    mulx    r8, rax, [r9+40]
    adcx    rsi, rax
    adox    r10, r8
    mulx    r8, rax, [r9+48]
    adcx    r10, rax
    adox    r11, r8
    mulx    r8, rax, [r9+56]
    adcx    r11, rax
    adox    rdi, r8
    mov     rax, 0
    adcx    rdi, rax
    adox    rdi, rax
    mov     [rcx+24], rbp           ; r[3]

    ; ============ Row 4: a × b[4] ============
    ; p[4..12] = r12, r13, r14, r15, rsi, r10, r11, rdi, rbp
    xor     rbp, rbp
    mov     rdx, [rbx+32]
    mulx    r8, rax, [r9]
    adcx    r12, rax
    adox    r13, r8
    mulx    r8, rax, [r9+8]
    adcx    r13, rax
    adox    r14, r8
    mulx    r8, rax, [r9+16]
    adcx    r14, rax
    adox    r15, r8
    mulx    r8, rax, [r9+24]
    adcx    r15, rax
    adox    rsi, r8
    mulx    r8, rax, [r9+32]
    adcx    rsi, rax
    adox    r10, r8
    mulx    r8, rax, [r9+40]
    adcx    r10, rax
    adox    r11, r8
    mulx    r8, rax, [r9+48]
    adcx    r11, rax
    adox    rdi, r8
    mulx    r8, rax, [r9+56]
    adcx    rdi, rax
    adox    rbp, r8
    mov     rax, 0
    adcx    rbp, rax
    adox    rbp, rax
    mov     [rcx+32], r12           ; r[4]

    ; ============ Row 5: a × b[5] ============
    ; p[5..13] = r13, r14, r15, rsi, r10, r11, rdi, rbp, r12
    xor     r12, r12
    mov     rdx, [rbx+40]
    mulx    r8, rax, [r9]
    adcx    r13, rax
    adox    r14, r8
    mulx    r8, rax, [r9+8]
    adcx    r14, rax
    adox    r15, r8
    mulx    r8, rax, [r9+16]
    adcx    r15, rax
    adox    rsi, r8
    mulx    r8, rax, [r9+24]
    adcx    rsi, rax
    adox    r10, r8
    mulx    r8, rax, [r9+32]
    adcx    r10, rax
    adox    r11, r8
    mulx    r8, rax, [r9+40]
    adcx    r11, rax
    adox    rdi, r8
    mulx    r8, rax, [r9+48]
    adcx    rdi, rax
    adox    rbp, r8
    mulx    r8, rax, [r9+56]
    adcx    rbp, rax
    adox    r12, r8
    mov     rax, 0
    adcx    r12, rax
    adox    r12, rax
    mov     [rcx+40], r13           ; r[5]

    ; ============ Row 6: a × b[6] ============
    ; p[6..14] = r14, r15, rsi, r10, r11, rdi, rbp, r12, r13
    xor     r13, r13
    mov     rdx, [rbx+48]
    mulx    r8, rax, [r9]
    adcx    r14, rax
    adox    r15, r8
    mulx    r8, rax, [r9+8]
    adcx    r15, rax
    adox    rsi, r8
    mulx    r8, rax, [r9+16]
    adcx    rsi, rax
    adox    r10, r8
    mulx    r8, rax, [r9+24]
    adcx    r10, rax
    adox    r11, r8
    mulx    r8, rax, [r9+32]
    adcx    r11, rax
    adox    rdi, r8
    mulx    r8, rax, [r9+40]
    adcx    rdi, rax
    adox    rbp, r8
    mulx    r8, rax, [r9+48]
    adcx    rbp, rax
    adox    r12, r8
    mulx    r8, rax, [r9+56]
    adcx    r12, rax
    adox    r13, r8
    mov     rax, 0
    adcx    r13, rax
    adox    r13, rax
    mov     [rcx+48], r14           ; r[6]

    ; ============ Row 7: a × b[7] ============
    ; p[7..15] = r15, rsi, r10, r11, rdi, rbp, r12, r13, r14
    xor     r14, r14
    mov     rdx, [rbx+56]
    mulx    r8, rax, [r9]
    adcx    r15, rax
    adox    rsi, r8
    mulx    r8, rax, [r9+8]
    adcx    rsi, rax
    adox    r10, r8
    mulx    r8, rax, [r9+16]
    adcx    r10, rax
    adox    r11, r8
    mulx    r8, rax, [r9+24]
    adcx    r11, rax
    adox    rdi, r8
    mulx    r8, rax, [r9+32]
    adcx    rdi, rax
    adox    rbp, r8
    mulx    r8, rax, [r9+40]
    adcx    rbp, rax
    adox    r12, r8
    mulx    r8, rax, [r9+48]
    adcx    r12, rax
    adox    r13, r8
    mulx    r8, rax, [r9+56]
    adcx    r13, rax
    adox    r14, r8
    mov     rax, 0
    adcx    r14, rax
    adox    r14, rax
    mov     [rcx+56], r15           ; r[7]

    ; ============ Final stores: r[8..15] ============
    ; p[8..15] = rsi, r10, r11, rdi, rbp, r12, r13, r14
    mov     [rcx+64], rsi           ; r[8]
    mov     [rcx+72], r10           ; r[9]
    mov     [rcx+80], r11           ; r[10]
    mov     [rcx+88], rdi           ; r[11]
    mov     [rcx+96], rbp           ; r[12]
    mov     [rcx+104], r12          ; r[13]
    mov     [rcx+112], r13          ; r[14]
    mov     [rcx+120], r14          ; r[15]

    pop     r15
    pop     r14
    pop     r13
    pop     r12
    pop     rbp
    pop     rdi
    pop     rsi
    pop     rbx
    ret
mpn_mul_8x8_asm ENDP

; =====================================================================
; uint64_t mpn_addmul_1_small_asm(uint64_t* rp, const uint64_t* ap,
;                                   size_t n, uint64_t b)
;
; rp[0..n-1] += ap[0..n-1] * b
; return value: carry (rax)
; precondition: 1 <= n <= 4, CPU supporting BMI2
;
; Eliminates the 7 push/pop + loop control of the general mpn_addmul_1_mulx.
; Uses volatile registers only (no push/pop).
;
; Windows x64: rcx=rp, rdx=ap, r8=n, r9=b
; =====================================================================
mpn_addmul_1_small_asm PROC
    lea     r10, [am1_s_jt]
    jmp     qword ptr [r10 + r8*8 - 8]

am1_s_jt:
    dq      am1_s_1
    dq      am1_s_2
    dq      am1_s_3
    dq      am1_s_4

    ALIGN   16
am1_s_1:
    ; rp[0] += a[0]*b, return carry
    mov     r10, rdx            ; r10 = ap
    mov     rdx, r9             ; rdx = b (MULX)
    mulx    r11, rax, [r10]     ; r11:rax = a[0]*b
    add     rax, [rcx]
    adc     r11, 0
    mov     [rcx], rax
    mov     rax, r11            ; return carry
    ret

    ALIGN   16
am1_s_2:
    mov     r10, rdx            ; r10 = ap
    mov     rdx, r9             ; rdx = b
    mulx    r11, rax, [r10]     ; r11:rax = a[0]*b
    add     rax, [rcx]
    adc     r11, 0
    mov     [rcx], rax

    mulx    r8, rax, [r10+8]    ; r8:rax = a[1]*b
    add     rax, r11
    adc     r8, 0
    add     rax, [rcx+8]
    adc     r8, 0
    mov     [rcx+8], rax
    mov     rax, r8
    ret

    ALIGN   16
am1_s_3:
    mov     r10, rdx            ; r10 = ap
    mov     rdx, r9             ; rdx = b
    mulx    r11, rax, [r10]     ; a[0]*b
    add     rax, [rcx]
    adc     r11, 0
    mov     [rcx], rax

    mulx    r8, rax, [r10+8]    ; a[1]*b
    add     rax, r11
    adc     r8, 0
    add     rax, [rcx+8]
    adc     r8, 0
    mov     [rcx+8], rax

    mulx    r9, rax, [r10+16]   ; a[2]*b
    add     rax, r8
    adc     r9, 0
    add     rax, [rcx+16]
    adc     r9, 0
    mov     [rcx+16], rax
    mov     rax, r9
    ret

    ALIGN   16
am1_s_4:
    mov     r10, rdx            ; r10 = ap
    mov     rdx, r9             ; rdx = b
    mulx    r11, rax, [r10]     ; a[0]*b
    add     rax, [rcx]
    adc     r11, 0
    mov     [rcx], rax

    mulx    r8, rax, [r10+8]    ; a[1]*b
    add     rax, r11
    adc     r8, 0
    add     rax, [rcx+8]
    adc     r8, 0
    mov     [rcx+8], rax

    mulx    r9, rax, [r10+16]   ; a[2]*b
    add     rax, r8
    adc     r9, 0
    add     rax, [rcx+16]
    adc     r9, 0
    mov     [rcx+16], rax

    mulx    r11, rax, [r10+24]  ; a[3]*b
    add     rax, r9
    adc     r11, 0
    add     rax, [rcx+24]
    adc     r11, 0
    mov     [rcx+24], rax
    mov     rax, r11
    ret
mpn_addmul_1_small_asm ENDP

END
