; Copyright (C) 2026 Kiyotsugu Arai
; SPDX-License-Identifier: LGPL-3.0-or-later
;
; mpn_x64.asm — x86-64 optimized routines
;
; Functions (BMI2+ADX required):
;   mpn_addmul_1_mulx(rp, ap, n, b) -> carry
;   mpn_mul_1_mulx(rp, ap, n, b)    -> carry
;   mpn_submul_1_mulx(rp, ap, n, b) -> borrow
;   mpn_mul_basecase_mulx(rp, ap, an, bp, bn)
;   mpn_sqr_basecase_mulx(rp, ap, n)
;
; Functions (base x86-64 only):
;   mpn_add_n_asm(rp, ap, bp, n) -> carry
;   mpn_sub_n_asm(rp, ap, bp, n) -> borrow
;   mpn_lshift_asm(rp, ap, n, shift) -> overflow
;   mpn_rshift_asm(rp, ap, n, shift) -> underflow (shifted-out bits)
;
; Precondition: the CPU must support BMI2 + ADX (caller has done the CPUID check)
;
; Windows x64 calling convention:
;   rcx = 1st, rdx = 2nd, r8 = 3rd, r9 = 4th
;   Return value: rax
;   Non-volatile: rbx, rbp, rdi, rsi, r12-r15
;   Volatile: rax, rcx, rdx, r8, r9, r10, r11
;
; Advantages of MULX+ADCX/ADOX:
;   MULX does not touch flags. ADCX uses only CF, ADOX uses only OF.
;   Two independent carry chains decouple the data dependency and
;   shorten the critical path per iteration (~2 cycles vs ~4 cycles).
;
; Loop control:
;   Uses LEA (flag-preserving) and JRCXZ (flag-preserving), so
;   CF/OF are never clobbered inside the loop.

.code

; =====================================================================
; uint64_t mpn_addmul_1_mulx(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
;
; 4x unrolled version:
;   Issue 4 independent MULX ahead of time and fold the results with the
;   dual ADCX/ADOX carry chains. Improves ~2 cycles/elem of the 1x version → ~1.75 cycles/elem.
;   - ADCX chain (CF): lo += prev_hi (carry propagation of the product)
;   - ADOX chain (OF): lo += rp[i]   (accumulation of the memory value)
;   Collect the carry every 4 elements, reset CF/OF, then continue the loop.
;
; Register allocation (4x loop):
;   rdx = b (MULX implicit), rsi = ap, rbx = rp, rcx = remaining element count
;   r8  = 0 (constant), r9 = carry (between groups)
;   MULX outputs: r10:r11 (elem0), r12:r13 (elem1), r14:r15 (elem2), rdi:rax (elem3)
; =====================================================================
mpn_addmul_1_mulx PROC
    push    rbx
    push    rsi
    push    rdi
    push    r12
    push    r13
    push    r14
    push    r15

    mov     rsi, rdx            ; rsi = ap
    mov     rbx, rcx            ; rbx = rp
    mov     rcx, r8             ; rcx = n
    mov     rdx, r9             ; rdx = b (MULX multiplier)

    xor     r9d, r9d            ; carry = 0
    xor     r8d, r8d            ; r8 = 0 (constant)

    test    rcx, rcx
    jz      am1_zero

    cmp     rcx, 8
    jb      am1_check_4x

    ; ---- 8x main loop ----
    ; Let CF/OF flow naturally from Group 1 → Group 2 (no mid-group merge needed)
    ; Collect the carry only at the end of Group 2
    ALIGN   16
am1_8x_entry:
    ; --- Group 1: elements [0..3] ---
    xor     eax, eax            ; CF = 0, OF = 0

    mulx    r10, r11, [rsi]
    mulx    r12, r13, [rsi + 8]
    mulx    r14, r15, [rsi + 16]
    mulx    rdi, rax, [rsi + 24]

    adcx    r11, r9
    adox    r11, [rbx]
    mov     QWORD PTR [rbx], r11

    adcx    r13, r10
    adox    r13, [rbx + 8]
    mov     QWORD PTR [rbx + 8], r13

    adcx    r15, r12
    adox    r15, [rbx + 16]
    mov     QWORD PTR [rbx + 16], r15

    adcx    rax, r14
    adox    rax, [rbx + 24]
    mov     QWORD PTR [rbx + 24], rax

    ; CF, OF flow as-is into Group 2 (rdi = hi3)

    ; --- Group 2: elements [4..7] ---
    mulx    r10, r11, [rsi + 32]
    mulx    r12, r13, [rsi + 40]
    mulx    r14, r15, [rsi + 48]
    mulx    r9, rax, [rsi + 56]

    adcx    r11, rdi
    adox    r11, [rbx + 32]
    mov     QWORD PTR [rbx + 32], r11

    adcx    r13, r10
    adox    r13, [rbx + 40]
    mov     QWORD PTR [rbx + 40], r13

    adcx    r15, r12
    adox    r15, [rbx + 48]
    mov     QWORD PTR [rbx + 48], r15

    adcx    rax, r14
    adox    rax, [rbx + 56]
    mov     QWORD PTR [rbx + 56], rax

    ; collect carry (r9 = hi7, gather CF/OF)
    adcx    r9, r8
    adox    r9, r8

    lea     rsi, [rsi + 64]
    lea     rbx, [rbx + 64]
    sub     rcx, 8
    cmp     rcx, 8
    jae     am1_8x_entry

    ; ---- remaining 0-7 elements ----
am1_check_4x:
    cmp     rcx, 4
    jb      am1_tail_check

    ; ---- single 4x ----
    xor     eax, eax
    mulx    r10, r11, [rsi]
    mulx    r12, r13, [rsi + 8]
    mulx    r14, r15, [rsi + 16]
    mulx    rdi, rax, [rsi + 24]

    adcx    r11, r9
    adox    r11, [rbx]
    mov     QWORD PTR [rbx], r11
    adcx    r13, r10
    adox    r13, [rbx + 8]
    mov     QWORD PTR [rbx + 8], r13
    adcx    r15, r12
    adox    r15, [rbx + 16]
    mov     QWORD PTR [rbx + 16], r15
    adcx    rax, r14
    adox    rax, [rbx + 24]
    mov     QWORD PTR [rbx + 24], rax

    mov     r9, rdi
    adcx    r9, r8
    adox    r9, r8
    lea     rsi, [rsi + 32]
    lea     rbx, [rbx + 32]
    sub     rcx, 4

am1_tail_check:
    test    rcx, rcx
    jz      am1_done

    ; ---- 1x tail loop ----
am1_tail_loop:
    mulx    r10, r11, [rsi]
    add     r11, r9
    adc     r10, 0
    add     QWORD PTR [rbx], r11
    adc     r10, 0
    mov     r9, r10
    lea     rsi, [rsi + 8]
    lea     rbx, [rbx + 8]
    dec     rcx
    jnz     am1_tail_loop

am1_done:
    mov     rax, r9

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

am1_zero:
    xor     eax, eax
    pop     r15
    pop     r14
    pop     r13
    pop     r12
    pop     rdi
    pop     rsi
    pop     rbx
    ret
mpn_addmul_1_mulx ENDP


; =====================================================================
; uint64_t mpn_mul_1_mulx(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
;
; 4x unrolled version:
;   Unlike addmul_1, no addition into rp[i] is needed, so ADOX is unused.
;   Only a single MULX + ADCX carry chain (rdi/r12-r15 are not needed).
;   Issue 4 independent MULX ahead of time and propagate the product carry with ADCX.
;   Collect the carry every 4 elements and continue the loop.
; =====================================================================
mpn_mul_1_mulx PROC
    push    rbx
    push    rsi
    push    rdi

    mov     rsi, rdx            ; rsi = ap
    mov     rbx, rcx            ; rbx = rp
    mov     rcx, r8             ; rcx = n
    mov     rdx, r9             ; rdx = b

    xor     r9d, r9d            ; carry = 0

    test    rcx, rcx
    jz      m1_zero

    cmp     rcx, 8
    jb      m1_check_4x

    ; ---- 8x main loop ----
m1_8x_entry:
    xor     eax, eax            ; CF = 0

    ; Group 1: elements [0..3]
    mulx    r10, r11, [rsi]
    mulx    rax, r8, [rsi + 8]
    adcx    r11, r9
    mov     [rbx], r11
    adcx    r8, r10
    mov     [rbx + 8], r8
    mulx    r10, r11, [rsi + 16]
    mulx    r9, r8, [rsi + 24]
    adcx    r11, rax
    mov     [rbx + 16], r11
    adcx    r8, r10
    mov     [rbx + 24], r8

    ; Group 2: elements [4..7]
    mulx    r10, r11, [rsi + 32]
    mulx    rax, rdi, [rsi + 40]
    adcx    r11, r9
    mov     [rbx + 32], r11
    adcx    rdi, r10
    mov     [rbx + 40], rdi
    mulx    r10, r11, [rsi + 48]
    mulx    r9, rdi, [rsi + 56]
    adcx    r11, rax
    mov     [rbx + 48], r11
    adcx    rdi, r10
    mov     [rbx + 56], rdi

    ; collect carry
    mov     r8d, 0              ; flag-preserving
    adcx    r9, r8

    lea     rsi, [rsi + 64]
    lea     rbx, [rbx + 64]
    sub     rcx, 8
    cmp     rcx, 8
    jae     m1_8x_entry

    ; ---- remaining 4-7 elements: 4x block ----
m1_check_4x:
    cmp     rcx, 4
    jb      m1_tail_check

    xor     eax, eax            ; CF = 0
    mulx    r10, r11, [rsi]
    mulx    rax, r8, [rsi + 8]
    adcx    r11, r9
    mov     [rbx], r11
    adcx    r8, r10
    mov     [rbx + 8], r8
    mulx    r10, r11, [rsi + 16]
    mulx    r9, r8, [rsi + 24]
    adcx    r11, rax
    mov     [rbx + 16], r11
    adcx    r8, r10
    mov     [rbx + 24], r8
    mov     r8d, 0
    adcx    r9, r8
    lea     rsi, [rsi + 32]
    lea     rbx, [rbx + 32]
    sub     rcx, 4

    ; ---- Tail: remaining 0-3 elements ----
m1_tail_check:
    test    rcx, rcx
    jz      m1_done

m1_tail_loop:
    mulx    r10, r11, [rsi]
    add     r11, r9
    adc     r10, 0
    mov     [rbx], r11
    mov     r9, r10
    lea     rsi, [rsi + 8]
    lea     rbx, [rbx + 8]
    dec     rcx
    jnz     m1_tail_loop

m1_done:
    mov     rax, r9
    pop     rdi
    pop     rsi
    pop     rbx
    ret

m1_zero:
    xor     eax, eax
    pop     rdi
    pop     rsi
    pop     rbx
    ret
mpn_mul_1_mulx ENDP


; =====================================================================
; uint64_t mpn_submul_1_mulx(uint64_t* rp, const uint64_t* ap,
;                              size_t n, uint64_t b)
;
; rp[0..n-1] -= ap[0..n-1] * b
; Return value: borrow (carry/borrow)
;
; The ADCX/ADOX dual carry chains cannot be used (because SUB clobbers OF).
; Instead, a single MULX + ADD/ADC/SUB/ADC chain is 4x/8x unrolled.
; Since MULX is flag-preserving, issue it ahead to exploit the OOO pipeline.
;
; Operation per element (4 dependent instructions):
;   (hi, lo) = a[i] * b       ; MULX (flags unchanged, can pipeline)
;   lo += carry; c1 = CF       ; ADD
;   hi += c1                   ; ADC
;   rp[i] -= lo; c2 = CF       ; SUB
;   carry = hi + c2             ; ADC
;
; 8x unrolling reduces loop overhead to 0.5 cycles/elem.
; Issuing MULX pairs ahead lets OOO run the multiplies in parallel with the carry chain.
; =====================================================================
mpn_submul_1_mulx PROC
    push    rbx
    push    rsi

    mov     rsi, rdx            ; rsi = ap
    mov     rbx, rcx            ; rbx = rp
    mov     rcx, r8             ; rcx = n
    mov     rdx, r9             ; rdx = b (MULX multiplier)

    xor     r9d, r9d            ; carry = 0

    test    rcx, rcx
    jz      submul1_zero

    cmp     rcx, 8
    jb      submul1_check_4x

    ; ---- 8x main loop ----
submul1_8x_entry:
    ; --- Group 1: elements 0-3 ---
    mulx    r10, r11, [rsi]           ; hi0:lo0
    mulx    rax, r8, [rsi + 8]        ; hi1:lo1

    add     r11, r9                   ; lo0 += carry
    adc     r10, 0
    sub     QWORD PTR [rbx], r11
    adc     r10, 0

    add     r8, r10                   ; lo1 += carry1
    adc     rax, 0
    sub     QWORD PTR [rbx + 8], r8
    adc     rax, 0

    mulx    r10, r11, [rsi + 16]      ; hi2:lo2
    mulx    r9, r8, [rsi + 24]        ; hi3:lo3

    add     r11, rax                  ; lo2 += carry2
    adc     r10, 0
    sub     QWORD PTR [rbx + 16], r11
    adc     r10, 0

    add     r8, r10                   ; lo3 += carry3
    adc     r9, 0
    sub     QWORD PTR [rbx + 24], r8
    adc     r9, 0

    ; --- Group 2: elements 4-7 ---
    mulx    r10, r11, [rsi + 32]
    mulx    rax, r8, [rsi + 40]

    add     r11, r9
    adc     r10, 0
    sub     QWORD PTR [rbx + 32], r11
    adc     r10, 0

    add     r8, r10
    adc     rax, 0
    sub     QWORD PTR [rbx + 40], r8
    adc     rax, 0

    mulx    r10, r11, [rsi + 48]
    mulx    r9, r8, [rsi + 56]

    add     r11, rax
    adc     r10, 0
    sub     QWORD PTR [rbx + 48], r11
    adc     r10, 0

    add     r8, r10
    adc     r9, 0
    sub     QWORD PTR [rbx + 56], r8
    adc     r9, 0

    lea     rsi, [rsi + 64]
    lea     rbx, [rbx + 64]
    sub     rcx, 8
    cmp     rcx, 8
    jae     submul1_8x_entry

    ; ---- remaining 0-7 elements ----
submul1_check_4x:
    cmp     rcx, 4
    jb      submul1_tail_check

    ; ---- single 4x ----
    mulx    r10, r11, [rsi]
    mulx    rax, r8, [rsi + 8]

    add     r11, r9
    adc     r10, 0
    sub     QWORD PTR [rbx], r11
    adc     r10, 0

    add     r8, r10
    adc     rax, 0
    sub     QWORD PTR [rbx + 8], r8
    adc     rax, 0

    mulx    r10, r11, [rsi + 16]
    mulx    r9, r8, [rsi + 24]

    add     r11, rax
    adc     r10, 0
    sub     QWORD PTR [rbx + 16], r11
    adc     r10, 0

    add     r8, r10
    adc     r9, 0
    sub     QWORD PTR [rbx + 24], r8
    adc     r9, 0

    lea     rsi, [rsi + 32]
    lea     rbx, [rbx + 32]
    sub     rcx, 4

submul1_tail_check:
    test    rcx, rcx
    jz      submul1_done

    ; ---- 1x tail loop ----
submul1_tail_loop:
    mulx    r10, r11, [rsi]
    add     r11, r9
    adc     r10, 0
    sub     QWORD PTR [rbx], r11
    adc     r10, 0
    mov     r9, r10
    lea     rsi, [rsi + 8]
    lea     rbx, [rbx + 8]
    dec     rcx
    jnz     submul1_tail_loop

submul1_done:
    mov     rax, r9
    pop     rsi
    pop     rbx
    ret

submul1_zero:
    xor     eax, eax
    pop     rsi
    pop     rbx
    ret
mpn_submul_1_mulx ENDP


; =====================================================================
; void mpn_mul_basecase_mulx(uint64_t* rp, const uint64_t* ap,
;                             size_t an, const uint64_t* bp, size_t bn)
;
; rp[0..an+bn-1] = ap[0..an-1] * bp[0..bn-1]
; Precondition: an >= bn >= 1, rp must not overlap ap,bp
;
; The C++ version of mul_basecase calls mul_1 + addmul_1×(bn-1) individually, so
; the prologue/epilogue of each call (addmul_1 has 7 push/pop) dominates
; at small sizes. This function handles all rows with a single push/pop.
;
; Windows x64 calling convention:
;   rcx = rp, rdx = ap, r8 = an, r9 = bp, [rsp+40] = bn
;
; Register allocation:
;   Inner loop (mul_1/addmul_1 body):
;     rbx = rp working pointer, rsi = ap working pointer
;     rdx = multiplier b[j], rcx = inner counter
;     r9 = carry, r8 = 0 constant
;     r10:r11, r12:r13, r14:r15, rdi:rax = 4x MULX hi:lo pairs
;   Outer loop state (saved to stack during inner loop):
;     [rsp+0]  = rp_base
;     [rsp+8]  = ap_base
;     [rsp+16] = an
;     [rsp+24] = bp_base
;     [rsp+32] = row counter (remaining addmul rows)
; =====================================================================
mpn_mul_basecase_mulx PROC
    push    rbx
    push    rsi
    push    rdi
    push    r12
    push    r13
    push    r14
    push    r15
    push    rbp
    sub     rsp, 40                 ; local variables (5 qwords)
    ; Stack layout:
    ;   8 push × 8 = 64 bytes + sub 40 = 104 bytes
    ;   5th arg (bn): [rsp + 40 + 64 + 40] = [rsp + 144]

    ; save parameters
    mov     [rsp+0], rcx            ; rp_base
    mov     [rsp+8], rdx            ; ap_base
    mov     [rsp+16], r8            ; an
    mov     [rsp+24], r9            ; bp_base
    mov     rax, [rsp+144]          ; bn
    mov     [rsp+32], rax           ; row_counter

    ; ========== Row 0: mul_1 (r[0..an] = a * b[0]) ==========
    mov     rbx, rcx                ; rbx = rp
    mov     rsi, rdx                ; rsi = ap
    mov     rdx, [r9]               ; rdx = b[0] (MULX multiplier)
    mov     rcx, r8                 ; rcx = an (inner counter)
    xor     r9d, r9d                ; carry = 0
    xor     r8d, r8d                ; r8 = 0

    cmp     rcx, 8
    jb      bc_r0_check_4x

    ALIGN   16
bc_r0_8x:
    ; --- Group 1: elements [0..3] ---
    xor     eax, eax                ; CF = 0
    mulx    r10, r11, [rsi]
    mulx    rax, r8, [rsi + 8]
    adcx    r11, r9
    mov     [rbx], r11
    adcx    r8, r10
    mov     [rbx + 8], r8
    mulx    r10, r11, [rsi + 16]
    mulx    r9, r8, [rsi + 24]
    adcx    r11, rax
    mov     [rbx + 16], r11
    adcx    r8, r10
    mov     [rbx + 24], r8
    ; --- Group 2: elements [4..7] ---
    mulx    r10, r11, [rsi + 32]
    mulx    rax, rdi, [rsi + 40]
    adcx    r11, r9
    mov     [rbx + 32], r11
    adcx    rdi, r10
    mov     [rbx + 40], rdi
    mulx    r10, r11, [rsi + 48]
    mulx    r9, rdi, [rsi + 56]
    adcx    r11, rax
    mov     [rbx + 48], r11
    adcx    rdi, r10
    mov     [rbx + 56], rdi
    mov     r8d, 0                  ; flag-preserving (xor cannot be used since it clobbers CF)
    adcx    r9, r8
    lea     rsi, [rsi + 64]
    lea     rbx, [rbx + 64]
    sub     rcx, 8
    cmp     rcx, 8
    jae     bc_r0_8x

bc_r0_check_4x:
    cmp     rcx, 4
    jb      bc_r0_tail_check

    xor     eax, eax                ; CF = 0
    mulx    r10, r11, [rsi]
    mulx    rax, r8, [rsi + 8]
    adcx    r11, r9
    mov     [rbx], r11
    adcx    r8, r10
    mov     [rbx + 8], r8
    mulx    r10, r11, [rsi + 16]
    mulx    r9, r8, [rsi + 24]
    adcx    r11, rax
    mov     [rbx + 16], r11
    adcx    r8, r10
    mov     [rbx + 24], r8
    mov     r8d, 0
    adcx    r9, r8
    lea     rsi, [rsi + 32]
    lea     rbx, [rbx + 32]
    sub     rcx, 4

bc_r0_tail_check:
    test    rcx, rcx
    jz      bc_r0_done

bc_r0_tail:
    mulx    r10, r11, [rsi]
    add     r11, r9
    adc     r10, 0
    mov     [rbx], r11
    mov     r9, r10
    lea     rsi, [rsi + 8]
    lea     rbx, [rbx + 8]
    dec     rcx
    jnz     bc_r0_tail

bc_r0_done:
    mov     [rbx], r9               ; rp[an] = carry

    ; done if bn == 1
    dec     QWORD PTR [rsp+32]
    jz      bc_finish

    ; ========== Rows 1..bn-1: addmul_1 (outer loop) ==========
    mov     rbp, 1                  ; j = 1 (row index)

bc_outer:
    ; restore ap, an
    mov     rsi, [rsp+8]            ; ap_base
    mov     rcx, [rsp+16]           ; an

    ; rbx = rp + j*8 (inner base address of this row)
    mov     rbx, [rsp+0]
    lea     rbx, [rbx + rbp*8]

    ; load b[j]
    mov     rax, [rsp+24]           ; bp_base
    mov     rdx, [rax + rbp*8]      ; b[j] → MULX multiplier

    ; skip if b[j] == 0
    test    rdx, rdx
    jz      bc_zero_row

    ; ---- addmul_1 inner loop (8x unrolled) ----
    xor     r9d, r9d                ; carry = 0
    xor     r8d, r8d                ; 0 constant

    cmp     rcx, 8
    jb      bc_am_check_4x

    ALIGN   16
bc_am_8x:
    ; --- Group 1: elements [0..3] ---
    xor     eax, eax
    mulx    r10, r11, [rsi]
    mulx    r12, r13, [rsi + 8]
    mulx    r14, r15, [rsi + 16]
    mulx    rdi, rax, [rsi + 24]
    adcx    r11, r9
    adox    r11, [rbx]
    mov     [rbx], r11
    adcx    r13, r10
    adox    r13, [rbx + 8]
    mov     [rbx + 8], r13
    adcx    r15, r12
    adox    r15, [rbx + 16]
    mov     [rbx + 16], r15
    adcx    rax, r14
    adox    rax, [rbx + 24]
    mov     [rbx + 24], rax
    ; CF, OF flow as-is into Group 2 (rdi = hi3)

    ; --- Group 2: elements [4..7] ---
    mulx    r10, r11, [rsi + 32]
    mulx    r12, r13, [rsi + 40]
    mulx    r14, r15, [rsi + 48]
    mulx    r9, rax, [rsi + 56]
    adcx    r11, rdi
    adox    r11, [rbx + 32]
    mov     [rbx + 32], r11
    adcx    r13, r10
    adox    r13, [rbx + 40]
    mov     [rbx + 40], r13
    adcx    r15, r12
    adox    r15, [rbx + 48]
    mov     [rbx + 48], r15
    adcx    rax, r14
    adox    rax, [rbx + 56]
    mov     [rbx + 56], rax
    ; collect carry (r9 = hi7)
    adcx    r9, r8
    adox    r9, r8

    lea     rsi, [rsi + 64]
    lea     rbx, [rbx + 64]
    sub     rcx, 8
    cmp     rcx, 8
    jae     bc_am_8x

bc_am_check_4x:
    cmp     rcx, 4
    jb      bc_am_tail_check

    xor     eax, eax
    mulx    r10, r11, [rsi]
    mulx    r12, r13, [rsi + 8]
    mulx    r14, r15, [rsi + 16]
    mulx    rdi, rax, [rsi + 24]
    adcx    r11, r9
    adox    r11, [rbx]
    mov     [rbx], r11
    adcx    r13, r10
    adox    r13, [rbx + 8]
    mov     [rbx + 8], r13
    adcx    r15, r12
    adox    r15, [rbx + 16]
    mov     [rbx + 16], r15
    adcx    rax, r14
    adox    rax, [rbx + 24]
    mov     [rbx + 24], rax
    mov     r9, rdi
    adcx    r9, r8
    adox    r9, r8
    lea     rsi, [rsi + 32]
    lea     rbx, [rbx + 32]
    sub     rcx, 4

bc_am_tail_check:
    test    rcx, rcx
    jz      bc_am_done

bc_am_tail:
    mulx    r10, r11, [rsi]
    add     r11, r9
    adc     r10, 0
    add     QWORD PTR [rbx], r11
    adc     r10, 0
    mov     r9, r10
    lea     rsi, [rsi + 8]
    lea     rbx, [rbx + 8]
    dec     rcx
    jnz     bc_am_tail

bc_am_done:
    ; write the carry into rp[an+j]
    ; (rbx has reached rp + j*8 + an*8 = rp + (an+j)*8)
    mov     [rbx], r9
    jmp     bc_next_row

bc_zero_row:
    ; b[j] == 0: rp[an+j] = 0
    mov     rax, [rsp+0]            ; rp_base
    mov     rcx, [rsp+16]           ; an
    add     rcx, rbp                ; an + j
    mov     QWORD PTR [rax + rcx*8], 0

bc_next_row:
    inc     rbp
    dec     QWORD PTR [rsp+32]
    jnz     bc_outer

bc_finish:
    add     rsp, 40
    pop     rbp
    pop     r15
    pop     r14
    pop     r13
    pop     r12
    pop     rdi
    pop     rsi
    pop     rbx
    ret
mpn_mul_basecase_mulx ENDP

; =====================================================================
; uint64_t mpn_add_n_asm(uint64_t* rp, const uint64_t* ap,
;                         const uint64_t* bp, size_t n)
;
; rp[0..n-1] = ap[0..n-1] + bp[0..n-1]
; Return value: carry (0 or 1)
; Precondition: n >= 1
;
; BMI2/ADX not required. Uses only the base x86-64 ADC instruction.
; Whereas MSVC's _addcarry_u64 generates a CF↔general-register conversion every iteration,
; a pure ADC chain keeps the carry in the hardware flag throughout.
;
; 8x unrolled. Loop control:
;   - DEC (CF-preserving) + JNZ to manage the counter
;   - LEA (flag-preserving) to advance pointers
;   - JRCXZ (flag-preserving) for the main loop test
;
; Register allocation:
;   rbx = rp (callee-saved), rdx = ap, r8 = bp
;   r10 = n%8 (remainder counter), rcx = n/8 (main loop iterations)
; =====================================================================
mpn_add_n_asm PROC
    push    rbx
    mov     rbx, rcx                ; rbx = rp (free up rcx)

    mov     r10, r9
    and     r10, 7                  ; r10 = n % 8 (remainder element count)
    mov     rcx, r9
    shr     rcx, 3                  ; rcx = n / 8 (main loop iterations)

    clc                             ; CF = 0

    test    r10, r10                ; CF is 0 from clc, so clobbering it is harmless
    jz      addn_rem_done

addn_rem_loop:
    mov     rax, [rdx]
    adc     rax, [r8]
    mov     [rbx], rax
    lea     rbx, [rbx + 8]         ; lea is flag-preserving
    lea     rdx, [rdx + 8]
    lea     r8,  [r8  + 8]
    dec     r10                     ; dec preserves CF
    jnz     addn_rem_loop

addn_rem_done:
    ; CF is carried over from the remainder loop (0 if there is no remainder)
    jrcxz   addn_done               ; finish if rcx == 0 (flag-preserving!)

addn_main_loop:
    mov     rax, [rdx]
    adc     rax, [r8]
    mov     [rbx], rax

    mov     rax, [rdx + 8]
    adc     rax, [r8  + 8]
    mov     [rbx + 8], rax

    mov     rax, [rdx + 16]
    adc     rax, [r8  + 16]
    mov     [rbx + 16], rax

    mov     rax, [rdx + 24]
    adc     rax, [r8  + 24]
    mov     [rbx + 24], rax

    mov     rax, [rdx + 32]
    adc     rax, [r8  + 32]
    mov     [rbx + 32], rax

    mov     rax, [rdx + 40]
    adc     rax, [r8  + 40]
    mov     [rbx + 40], rax

    mov     rax, [rdx + 48]
    adc     rax, [r8  + 48]
    mov     [rbx + 48], rax

    mov     rax, [rdx + 56]
    adc     rax, [r8  + 56]
    mov     [rbx + 56], rax

    lea     rbx, [rbx + 64]
    lea     rdx, [rdx + 64]
    lea     r8,  [r8  + 64]
    dec     rcx                     ; dec preserves CF
    jnz     addn_main_loop

addn_done:
    setc    al
    movzx   eax, al
    pop     rbx
    ret
mpn_add_n_asm ENDP

; =====================================================================
; uint64_t mpn_sub_n_asm(uint64_t* rp, const uint64_t* ap,
;                         const uint64_t* bp, size_t n)
;
; rp[0..n-1] = ap[0..n-1] - bp[0..n-1]
; Return value: borrow (0 or 1)
; Precondition: n >= 1
;
; Same structure as mpn_add_n_asm. ADC → SBB. 8x unrolled.
; =====================================================================
mpn_sub_n_asm PROC
    push    rbx
    mov     rbx, rcx

    mov     r10, r9
    and     r10, 7
    mov     rcx, r9
    shr     rcx, 3

    clc

    test    r10, r10
    jz      subn_rem_done

subn_rem_loop:
    mov     rax, [rdx]
    sbb     rax, [r8]
    mov     [rbx], rax
    lea     rbx, [rbx + 8]
    lea     rdx, [rdx + 8]
    lea     r8,  [r8  + 8]
    dec     r10
    jnz     subn_rem_loop

subn_rem_done:
    jrcxz   subn_done

subn_main_loop:
    mov     rax, [rdx]
    sbb     rax, [r8]
    mov     [rbx], rax

    mov     rax, [rdx + 8]
    sbb     rax, [r8  + 8]
    mov     [rbx + 8], rax

    mov     rax, [rdx + 16]
    sbb     rax, [r8  + 16]
    mov     [rbx + 16], rax

    mov     rax, [rdx + 24]
    sbb     rax, [r8  + 24]
    mov     [rbx + 24], rax

    mov     rax, [rdx + 32]
    sbb     rax, [r8  + 32]
    mov     [rbx + 32], rax

    mov     rax, [rdx + 40]
    sbb     rax, [r8  + 40]
    mov     [rbx + 40], rax

    mov     rax, [rdx + 48]
    sbb     rax, [r8  + 48]
    mov     [rbx + 48], rax

    mov     rax, [rdx + 56]
    sbb     rax, [r8  + 56]
    mov     [rbx + 56], rax

    lea     rbx, [rbx + 64]
    lea     rdx, [rdx + 64]
    lea     r8,  [r8  + 64]
    dec     rcx
    jnz     subn_main_loop

subn_done:
    setc    al
    movzx   eax, al
    pop     rbx
    ret
mpn_sub_n_asm ENDP

; =====================================================================
; uint64_t mpn_lshift_asm(uint64_t* rp, const uint64_t* ap,
;                          size_t n, unsigned shift)
;
; rp[0..n-1] = ap[0..n-1] << shift
; Return value: the overflowed most-significant bits
; Precondition: n >= 1, 1 <= shift <= 63
;
; Algorithm:
;   Process in the HIGH→LOW direction with the SHLD instruction.
;   SHLD dst, src, cl: dst = (dst << cl) | (src >> (64-cl))
;   For each element, concatenate-shift ap[i] and ap[i-1] and store the result in rp[i].
;   4x unrolled + remainder loop. in-place safe (HIGH→LOW).
;
; Register allocation:
;   rbx = rp (callee-saved), rdx = ap
;   cl = shift, rdi = previous ap[i] (callee-saved)
;   r10 = overflow (return value), r9 = remainder counter, r8 = main loop iterations
; =====================================================================
mpn_lshift_asm PROC
    push    rbx
    push    rdi
    mov     rbx, rcx                ; rbx = rp
    mov     rcx, r9                 ; cl = shift

    ; set pointers to the tail (process HIGH→LOW)
    lea     rbx, [rbx + r8*8 - 8]   ; rbx = &rp[n-1]
    lea     rdx, [rdx + r8*8 - 8]   ; rdx = &ap[n-1]

    ; compute overflow = ap[n-1] >> (64-shift) with SHLD
    mov     rdi, [rdx]               ; rdi = ap[n-1]
    xor     r10d, r10d
    shld    r10, rdi, cl             ; r10 = rdi >> (64-shift) = overflow

    ; process n-1 pairs
    dec     r8                       ; r8 = n-1
    jz      lsh_last                 ; n==1 → only the last element

    ; remainder = (n-1) % 4
    mov     r9, r8
    and     r9, 3                    ; r9 = remainder counter
    shr     r8, 2                    ; r8 = main loop iterations

    test    r9, r9
    jz      lsh_rem_done

lsh_rem_loop:
    mov     rax, [rdx - 8]          ; rax = ap[i-1]
    shld    rdi, rax, cl            ; rdi = (ap[i]<<shift)|(ap[i-1]>>(64-shift))
    mov     [rbx], rdi              ; rp[i] = result
    mov     rdi, rax                ; rdi = ap[i-1] for next iteration
    lea     rbx, [rbx - 8]
    lea     rdx, [rdx - 8]
    dec     r9
    jnz     lsh_rem_loop

lsh_rem_done:
    test    r8, r8
    jz      lsh_last

lsh_main_loop:
    mov     rax, [rdx - 8]
    shld    rdi, rax, cl
    mov     [rbx], rdi
    mov     rdi, rax

    mov     rax, [rdx - 16]
    shld    rdi, rax, cl
    mov     [rbx - 8], rdi
    mov     rdi, rax

    mov     rax, [rdx - 24]
    shld    rdi, rax, cl
    mov     [rbx - 16], rdi
    mov     rdi, rax

    mov     rax, [rdx - 32]
    shld    rdi, rax, cl
    mov     [rbx - 24], rdi
    mov     rdi, rax

    lea     rbx, [rbx - 32]
    lea     rdx, [rdx - 32]
    dec     r8
    jnz     lsh_main_loop

lsh_last:
    ; rp[0] = rdi << shift (rdi is the last-read ap[0])
    shl     rdi, cl
    mov     [rbx], rdi

    mov     rax, r10                ; return overflow
    pop     rdi
    pop     rbx
    ret
mpn_lshift_asm ENDP

; =====================================================================
; uint64_t mpn_rshift_asm(uint64_t* rp, const uint64_t* ap,
;                          size_t n, unsigned shift)
;
; rp[0..n-1] = ap[0..n-1] >> shift
; Return value: the bits overflowed off the bottom (ap[0] << (64-shift))
; Precondition: n >= 1, 1 <= shift <= 63
;
; Algorithm:
;   Process in the LOW→HIGH direction with the SHRD instruction.
;   SHRD dst, src, cl: dst = (src << (64-cl)) | (dst >> cl)
;   4x unrolled + remainder loop. in-place safe (LOW→HIGH).
; =====================================================================
mpn_rshift_asm PROC
    push    rbx
    push    rdi
    mov     rbx, rcx                ; rbx = rp
    mov     rcx, r9                 ; cl = shift

    ; compute underflow = ap[0] << (64-shift) with SHRD
    mov     rdi, [rdx]               ; rdi = ap[0]
    xor     r10d, r10d
    shrd    r10, rdi, cl             ; r10 = rdi << (64-shift) = underflow

    ; process n-1 pairs
    dec     r8                       ; r8 = n-1
    jz      rsh_last

    mov     r9, r8
    and     r9, 3
    shr     r8, 2

    test    r9, r9
    jz      rsh_rem_done

rsh_rem_loop:
    mov     rax, [rdx + 8]          ; rax = ap[i+1]
    shrd    rdi, rax, cl            ; rdi = (ap[i+1]<<(64-shift))|(ap[i]>>shift)
    mov     [rbx], rdi
    mov     rdi, rax
    lea     rbx, [rbx + 8]
    lea     rdx, [rdx + 8]
    dec     r9
    jnz     rsh_rem_loop

rsh_rem_done:
    test    r8, r8
    jz      rsh_last

rsh_main_loop:
    mov     rax, [rdx + 8]
    shrd    rdi, rax, cl
    mov     [rbx], rdi
    mov     rdi, rax

    mov     rax, [rdx + 16]
    shrd    rdi, rax, cl
    mov     [rbx + 8], rdi
    mov     rdi, rax

    mov     rax, [rdx + 24]
    shrd    rdi, rax, cl
    mov     [rbx + 16], rdi
    mov     rdi, rax

    mov     rax, [rdx + 32]
    shrd    rdi, rax, cl
    mov     [rbx + 24], rdi
    mov     rdi, rax

    lea     rbx, [rbx + 32]
    lea     rdx, [rdx + 32]
    dec     r8
    jnz     rsh_main_loop

rsh_last:
    ; rp[n-1] = rdi >> shift
    shr     rdi, cl
    mov     [rbx], rdi

    mov     rax, r10                ; return underflow
    pop     rdi
    pop     rbx
    ret
mpn_rshift_asm ENDP

; =====================================================================
; void mpn_sqr_basecase_mulx(uint64_t* rp, const uint64_t* ap, size_t n)
;
; rp[0..2n-1] = ap[0..n-1]²
; Precondition: n >= 1, rp and ap do not overlap, rp need not be pre-zeroed
;
; Algorithm (3 steps):
;   Step 1: Off-diagonal — upper triangle Σ_{i<j} a[i]*a[j]*B^(i+j)
;           Row 0: mul_1, Rows 1..n-2: addmul_1 (all rows inlined)
;   Step 2: double — inline lshift 1 via an ADC chain
;   Step 3: Diagonal — add a[i]² into r[2i..2i+1] with MULX+ADCX
;
; By symmetry, the multiply count = n(n-1)/2 + n (vs n² for ordinary multiplication)
; Inlining all rows removes addmul_1's push/pop (14 instructions/row)
;
; rcx = rp, rdx = ap, r8 = n
; =====================================================================
mpn_sqr_basecase_mulx PROC
    push    rbx
    push    rsi
    push    rdi
    push    r12
    push    r13
    push    r14
    push    r15
    push    rbp
    sub     rsp, 40                 ; 5 local qwords
    ; Stack: 8 push * 8 = 64 + sub 40 = 104 bytes

    mov     [rsp+0], rcx            ; rp_base
    mov     [rsp+8], rdx            ; ap_base
    mov     [rsp+16], r8            ; n

    ; ===== Step 0: Zero r[0..2n-1] (rep stosq) =====
    mov     rdi, rcx                ; rdi = rp
    mov     rcx, r8
    shl     rcx, 1                  ; rcx = 2n
    xor     eax, eax
    rep     stosq

    ; ===== n == 1: r[0..1] = a[0]² =====
    mov     rcx, [rsp+16]
    cmp     rcx, 1
    ja      sqb_step1

    mov     rsi, [rsp+8]
    mov     rdx, [rsi]
    mulx    r10, r11, rdx
    mov     rbx, [rsp+0]
    mov     [rbx], r11
    mov     [rbx+8], r10
    jmp     sqb_done

    ; ========== Step 1: Off-diagonal ==========
sqb_step1:
    ; --- Row 0: mul_1 (r[1..n-1] = a[1..n-1]*a[0], carry → r[n]) ---
    mov     rbx, [rsp+0]
    lea     rbx, [rbx + 8]         ; rbx = rp + 1
    mov     rsi, [rsp+8]
    mov     rdx, [rsi]              ; a[0]
    lea     rsi, [rsi + 8]          ; ap + 1
    mov     rcx, [rsp+16]
    dec     rcx                     ; n - 1
    xor     r9d, r9d
    xor     r8d, r8d

    cmp     rcx, 4
    jb      sqb_r0_tail

sqb_r0_4x:
    xor     eax, eax
    mulx    r10, r11, [rsi]
    mulx    rax, r8, [rsi + 8]
    adcx    r11, r9
    mov     [rbx], r11
    adcx    r8, r10
    mov     [rbx + 8], r8
    mulx    r10, r11, [rsi + 16]
    mulx    r9, r8, [rsi + 24]
    adcx    r11, rax
    mov     [rbx + 16], r11
    adcx    r8, r10
    mov     [rbx + 24], r8
    mov     r8d, 0
    adcx    r9, r8
    lea     rsi, [rsi + 32]
    lea     rbx, [rbx + 32]
    sub     rcx, 4
    cmp     rcx, 4
    jae     sqb_r0_4x

    test    rcx, rcx
    jz      sqb_r0_done

sqb_r0_tail:
    mulx    r10, r11, [rsi]
    add     r11, r9
    adc     r10, 0
    mov     [rbx], r11
    mov     r9, r10
    lea     rsi, [rsi + 8]
    lea     rbx, [rbx + 8]
    dec     rcx
    jnz     sqb_r0_tail

sqb_r0_done:
    mov     [rbx], r9               ; r[n] = carry

    ; n == 2: no more rows
    mov     rax, [rsp+16]
    cmp     rax, 2
    jbe     sqb_step2

    ; --- Rows 1..n-2: addmul_1 (inline, no push/pop needed) ---
    mov     rbp, 1                  ; i = 1
    mov     rax, [rsp+16]
    sub     rax, 2                  ; row_counter = n - 2
    mov     [rsp+32], rax

sqb_outer:
    ; dst = rp + (2*i + 1)*8
    mov     rbx, [rsp+0]
    lea     rax, [rbp*2 + 1]
    lea     rbx, [rbx + rax*8]

    ; src = ap + (i+1)*8
    mov     rsi, [rsp+8]
    lea     rsi, [rsi + rbp*8 + 8]

    ; multiplier = a[i]
    mov     rax, [rsp+8]
    mov     rdx, [rax + rbp*8]

    ; inner count = n - i - 1
    mov     rcx, [rsp+16]
    sub     rcx, rbp
    dec     rcx

    ; skip if a[i] == 0
    test    rdx, rdx
    jz      sqb_zero_row

    xor     r9d, r9d
    xor     r8d, r8d

    cmp     rcx, 8
    jb      sqb_am_4x_check

    ; ---- 8x unrolled addmul_1 (MULX+ADCX/ADOX) ----
sqb_am_8x:
    xor     eax, eax
    mulx    r10, r11, [rsi]
    mulx    r12, r13, [rsi + 8]
    mulx    r14, r15, [rsi + 16]
    mulx    rdi, rax, [rsi + 24]
    adcx    r11, r9
    adox    r11, [rbx]
    mov     [rbx], r11
    adcx    r13, r10
    adox    r13, [rbx + 8]
    mov     [rbx + 8], r13
    adcx    r15, r12
    adox    r15, [rbx + 16]
    mov     [rbx + 16], r15
    adcx    rax, r14
    adox    rax, [rbx + 24]
    mov     [rbx + 24], rax
    mov     r9, rdi
    adcx    r9, r8
    adox    r9, r8

    mulx    r10, r11, [rsi + 32]
    mulx    r12, r13, [rsi + 40]
    mulx    r14, r15, [rsi + 48]
    mulx    rdi, rax, [rsi + 56]
    adcx    r11, r9
    adox    r11, [rbx + 32]
    mov     [rbx + 32], r11
    adcx    r13, r10
    adox    r13, [rbx + 40]
    mov     [rbx + 40], r13
    adcx    r15, r12
    adox    r15, [rbx + 48]
    mov     [rbx + 48], r15
    adcx    rax, r14
    adox    rax, [rbx + 56]
    mov     [rbx + 56], rax
    mov     r9, rdi
    adcx    r9, r8
    adox    r9, r8

    lea     rsi, [rsi + 64]
    lea     rbx, [rbx + 64]
    sub     rcx, 8
    cmp     rcx, 8
    jae     sqb_am_8x

sqb_am_4x_check:
    cmp     rcx, 4
    jb      sqb_am_tail_check

    xor     eax, eax
    mulx    r10, r11, [rsi]
    mulx    r12, r13, [rsi + 8]
    mulx    r14, r15, [rsi + 16]
    mulx    rdi, rax, [rsi + 24]
    adcx    r11, r9
    adox    r11, [rbx]
    mov     [rbx], r11
    adcx    r13, r10
    adox    r13, [rbx + 8]
    mov     [rbx + 8], r13
    adcx    r15, r12
    adox    r15, [rbx + 16]
    mov     [rbx + 16], r15
    adcx    rax, r14
    adox    rax, [rbx + 24]
    mov     [rbx + 24], rax
    mov     r9, rdi
    adcx    r9, r8
    adox    r9, r8
    lea     rsi, [rsi + 32]
    lea     rbx, [rbx + 32]
    sub     rcx, 4

sqb_am_tail_check:
    test    rcx, rcx
    jz      sqb_am_done

sqb_am_tail:
    mulx    r10, r11, [rsi]
    add     r11, r9
    adc     r10, 0
    add     QWORD PTR [rbx], r11
    adc     r10, 0
    mov     r9, r10
    lea     rsi, [rsi + 8]
    lea     rbx, [rbx + 8]
    dec     rcx
    jnz     sqb_am_tail

sqb_am_done:
    mov     [rbx], r9               ; carry → r[n+i]
    jmp     sqb_next_row

sqb_zero_row:
    ; a[i] == 0: r[n+i] already zeroed in Step 0

sqb_next_row:
    inc     rbp
    dec     QWORD PTR [rsp+32]
    jnz     sqb_outer

    ; ========== Step 2: lshift by 1 (inline ADC chain) ==========
sqb_step2:
    mov     rbx, [rsp+0]           ; rp
    mov     rcx, [rsp+16]          ; n (2 words/iter × n = 2n words)
    clc                             ; CF = 0

sqb_lsh:
    mov     rax, [rbx]
    adc     rax, rax
    mov     [rbx], rax
    mov     rax, [rbx + 8]
    adc     rax, rax
    mov     [rbx + 8], rax
    lea     rbx, [rbx + 16]        ; flag-preserving
    dec     rcx                     ; DEC is CF-preserving
    jnz     sqb_lsh

    ; ========== Step 3: Diagonal — add a[i]² into r[2i..2i+1] ==========
    mov     rbx, [rsp+0]           ; rp
    mov     rsi, [rsp+8]           ; ap
    mov     rcx, [rsp+16]          ; n

    ; n/2 groups + handle 1 leading element if odd
    mov     r8, rcx
    shr     rcx, 1                  ; rcx = n / 2
    and     r8, 1                   ; r8 = n % 2

    xor     eax, eax                ; CF = 0 (start of ADCX chain)

    test    r8, r8                  ; CF = 0 (test clears CF)
    jz      sqb_diag_even

    ; n is odd: leading 1 element
    mov     rdx, [rsi]
    mulx    r10, r11, rdx
    adcx    r11, [rbx]
    mov     [rbx], r11
    adcx    r10, [rbx + 8]
    mov     [rbx + 8], r10
    lea     rsi, [rsi + 8]
    lea     rbx, [rbx + 16]

sqb_diag_even:
    jrcxz   sqb_done                ; JRCXZ is flag-preserving

sqb_diag_2x:
    mov     rdx, [rsi]
    mulx    r10, r11, rdx           ; a[i]²
    adcx    r11, [rbx]
    mov     [rbx], r11
    adcx    r10, [rbx + 8]
    mov     [rbx + 8], r10

    mov     rdx, [rsi + 8]
    mulx    r10, r11, rdx           ; a[i+1]²
    adcx    r11, [rbx + 16]
    mov     [rbx + 16], r11
    adcx    r10, [rbx + 24]
    mov     [rbx + 24], r10

    lea     rsi, [rsi + 16]
    lea     rbx, [rbx + 32]
    dec     rcx                     ; DEC is CF-preserving
    jnz     sqb_diag_2x

sqb_done:
    add     rsp, 40
    pop     rbp
    pop     r15
    pop     r14
    pop     r13
    pop     r12
    pop     rdi
    pop     rsi
    pop     rbx
    ret
mpn_sqr_basecase_mulx ENDP

END
