; Copyright (C) 2026 Kiyotsugu Arai
; SPDX-License-Identifier: LGPL-3.0-or-later
;
; mpn_x64_div.asm — ASM implementation of the div_basecase 3-by-2 loop
;
; Functions:
;   mpn_sbpi1_div_qr_asm(qp, ap, qn, dp, dn, dinv3)
;
; Windows x64 calling convention:
;   rcx = qp, rdx = ap, r8 = qn, r9 = dp
;   [rsp+40] = dn (≥3), [rsp+48] = dinv3
;
; Stack frame (eliminates push/pop inside the loop):
;   [rsp+0..7]   = q0 temporary save
;   [rsp+8..15]  = n0 temporary save
;   [rsp+16..23] = n1 temporary save

; Stack frame:
;   [rsp+0..31]   = shadow space for CALL (32 bytes, usable by callee)
;   [rsp+32..39]  = q0 temporary save
;   [rsp+40..47]  = n0 temporary save
;   [rsp+48..55]  = n1 temporary save
;   [rsp+56..63]  = dn temporary save (for qmax)
FRAME_SIZE = 64

EXTERN mpn_submul_1_mulx:PROC

.code

mpn_sbpi1_div_qr_asm PROC
    push    rbx
    push    rsi
    push    rdi
    push    rbp
    push    r12
    push    r13
    push    r14
    push    r15
    sub     rsp, FRAME_SIZE         ; allocate stack frame

    ; Retrieve parameters
    ; stack args as seen from the caller's rsp:
    ;   5th arg (dn)    = [original_rsp + 32]
    ;   6th arg (dinv3) = [original_rsp + 40]
    ; callee side: 8 pushes (64) + sub rsp,FRAME_SIZE (64) = 128 bytes
    ; therefore: 5th arg = [rsp + 128 + 32] = [rsp + 160]
    ;            6th arg = [rsp + 128 + 40] = [rsp + 168]
    mov     rbx, rcx                ; rbx = qp
    mov     rsi, rdx                ; rsi = ap
    mov     r15, r8                 ; r15 = qn
    mov     rdi, r9                 ; rdi = dp
    mov     rbp, [rsp + FRAME_SIZE + 8*8 + 40]  ; rbp = dn (5th arg)
    mov     r14, [rsp + FRAME_SIZE + 8*8 + 48]  ; r14 = dinv3 (6th arg)

    ; Load constants
    mov     r12, [rdi + rbp*8 - 8]  ; r12 = d1 = dp[dn-1]
    mov     r13, [rdi + rbp*8 - 16] ; r13 = d0 = dp[dn-2]
    sub     rbp, 2                  ; rbp = dn-2
    ; r15 = qn = an-bn
    ; ASM processes j = qn downto 0 (qn+1 iterations).
    ; The j=qn iteration overlaps with C++'s qh check, but
    ; thanks to the np[nn]=0 sentinel, q_hat ≈ 0 makes it effectively a no-op.
    ; However, since np[nn-2], np[nn-1] are overwritten with the udiv_qrnnd_3by2 result,
    ; internal calls such as dcpi1 require protecting qp[qn].

outer_loop:
    ; r15 = current j (qn-1, qn-2, ..., 0)
    ; if j < 0, finish
    test    r15, r15
    js      done

    ; Load u2, u1, u0
    lea     rax, [r15 + rbp + 2]
    mov     r8, [rsi + rax*8]       ; r8 = u2
    mov     r9, [rsi + rax*8 - 8]   ; r9 = u1
    mov     r10, [rsi + rax*8 - 16] ; r10 = u0

    ; Special case: u2==d1 && u1==d0 → q=MAX
    cmp     r8, r12
    jne     normal_div
    cmp     r9, r13
    jne     normal_div

    ; === q_hat = UINT64_MAX, submul_1(a+j, dp, dn, UINT64_MAX) ===
    mov     r11, -1
    mov     rdx, r11                ; rdx = UINT64_MAX (MULX)
    xor     r9d, r9d                ; r9 = carry = 0
    xor     ecx, ecx                ; ecx = i = 0

    ; Save dn = rbp + 2 on the stack (due to register shortage)
    lea     rax, [rbp + 2]
    mov     [rsp + 56], rax         ; [rsp+24] = dn

qmax_submul:
    cmp     rcx, [rsp + 56]         ; i < dn ?
    jge     qmax_submul_done
    mulx    r8, r10, [rdi + rcx*8]  ; r8:r10 = dp[i] * UINT64_MAX
    add     r10, r9                 ; lo += carry
    adc     r8, 0                   ; hi += carry_out
    lea     rax, [r15 + rcx]        ; rax = j + i (preserve r9/r8)
    mov     r9, [rsi + rax*8]       ; r9 = a[j+i]
    sub     r9, r10                 ; a[j+i] -= lo
    mov     [rsi + rax*8], r9       ; store
    adc     r8, 0                   ; borrow → hi
    mov     r9, r8                  ; r9 = new carry
    inc     rcx
    jmp     qmax_submul

qmax_submul_done:
    ; carry in r9 → a[j+dn] -= carry
    lea     rax, [r15 + rbp + 2]    ; rax = j + dn
    mov     r8, [rsi + rax*8]
    sub     r8, r9
    mov     [rsi + rax*8], r8
    jnc     store_q
    ; add-back (extremely rare)
    dec     r11
    xor     ecx, ecx
    clc
    ; addback: loop dn times + carry propagation (the dn+1-th limb)
    mov     rcx, [rsp + 56]         ; rcx = dn (loop counter)
    xor     r8d, r8d                ; r8 = index = 0
qmax_addback:
    lea     rax, [r15 + r8]
    mov     r9, [rdi + r8*8]
    adc     r9, [rsi + rax*8]
    mov     [rsi + rax*8], r9
    lea     r8, [r8 + 1]            ; does not destroy CF
    dec     rcx                     ; preserves CF (dec does not destroy CF)
    jnz     qmax_addback
    ; propagate carry to a[j+dn] (equivalent to C++'s add(np+i, np+i, dn+1, dp, dn))
    lea     rax, [r15 + r8]         ; rax = j + dn
    adc     QWORD PTR [rsi + rax*8], 0

normal_div:
    ; === udiv_qrnnd_3by2 ===
    ; {q1, q0} = u2 * dinv + {u2, u1}
    mov     rdx, r14
    mulx    rcx, rax, r8            ; rcx:rax = u2 * dinv3
    add     rax, r9                 ; q0 = lo + u1
    adc     rcx, r8                 ; q1 = hi + u2 + carry
    mov     r11, rcx                ; r11 = q1
    mov     [rsp+32], rax              ; save q0

    ; r1 = u1 - d1*q1
    mov     rax, r12
    imul    rax, r11
    mov     rcx, r9
    sub     rcx, rax                ; rcx = r1

    ; {r1, r0} = {r1, u0} - {d1, d0}
    mov     rax, r10
    sub     rax, r13
    sbb     rcx, r12                ; rcx=r1, rax=r0

    ; {r1, r0} -= d0 * q1
    mov     rdx, r13
    mulx    r8, r9, r11             ; r8:r9 = d0 * q1
    sub     rax, r9
    sbb     rcx, r8                 ; rcx=r1, rax=r0

    inc     r11                     ; q1++

    ; Conditional correction
    mov     r8, [rsp+32]               ; r8 = q0
    cmp     rcx, r8
    jb      no_adjust
    dec     r11
    add     rax, r13
    adc     rcx, r12
no_adjust:
    cmp     rcx, r12
    jb      div_done
    ja      final_adjust
    cmp     rax, r13
    jb      div_done
final_adjust:
    inc     r11
    sub     rax, r13
    sbb     rcx, r12

div_done:
    ; rcx=n1, rax=n0, r11=q_hat
    ; execute submul_1(a+j, dp, dn-2, q_hat) via the C++ ASM function
    mov     [rsp+40], rax            ; save n0
    mov     [rsp+48], rcx           ; save n1
    ; r11 (q_hat) is volatile, so save it before CALL
    mov     [rsp+32], r11           ; save q_hat (reuse, since q0 is no longer needed)

    ; mpn_submul_1_mulx(rp=a+j, ap=dp, n=dn-2, b=q_hat)
    ; Windows x64: rcx=rp, rdx=ap, r8=n, r9=b
    lea     rcx, [rsi + r15*8]      ; rcx = a + j
    mov     rdx, rdi                ; rdx = dp
    mov     r8, rbp                 ; r8 = dn-2
    mov     r9, r11                 ; r9 = q_hat
    call    mpn_submul_1_mulx
    ; rax = borrow (submul_1 return value)
    mov     r9, rax                 ; r9 = submul borrow
    mov     r11, [rsp+32]           ; r11 = restore q_hat

    ; n0 -= borrow, n1 -= carry
    mov     rax, [rsp+40]            ; n0
    mov     rcx, [rsp+48]           ; n1
    sub     rax, r9                  ; n0 -= submul_borrow
    sbb     rcx, 0                   ; n1 -= CF (= cy0 = n0 < borrow)
    ; CF = cy1 = (n1 was < cy0) — determine whether add-back is needed
    ; save CF before writing to a[]
    setc    r9b                      ; r9b = cy1

    ; a[j+dn-2] = n0, a[j+dn-1] = n1, a[j+dn] = 0
    lea     r10, [r15 + rbp]
    mov     [rsi + r10*8], rax
    mov     [rsi + r10*8 + 8], rcx
    mov     QWORD PTR [rsi + r10*8 + 16], 0

    ; add-back: cy1 != 0 ?
    test    r9b, r9b
    jz      store_q
    dec     r11
    lea     rcx, [rbp + 2]          ; rcx = dn (loop counter)
    xor     r8d, r8d                ; r8 = index = 0
    clc
addback_loop:
    lea     r10, [r15 + r8]
    mov     rax, [rdi + r8*8]
    adc     rax, [rsi + r10*8]
    mov     [rsi + r10*8], rax
    lea     r8, [r8 + 1]
    dec     rcx                     ; dec does not destroy CF
    jnz     addback_loop

store_q:
    mov     [rbx + r15*8], r11      ; q[j] = q_hat
    dec     r15                     ; j--
    jmp     outer_loop

done:
    add     rsp, FRAME_SIZE
    pop     r15
    pop     r14
    pop     r13
    pop     r12
    pop     rbp
    pop     rdi
    pop     rsi
    pop     rbx
    ret
mpn_sbpi1_div_qr_asm ENDP

END
