; Copyright (C) 2026 Kiyotsugu Arai
; SPDX-License-Identifier: LGPL-3.0-or-later
;
; mpn_x64_div.asm — div_basecase 3-by-2 ループの ASM 実装
;
; 関数:
;   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
;
; スタックフレーム (ループ内の push/pop を排除):
;   [rsp+0..7]   = q0 一時保存
;   [rsp+8..15]  = n0 一時保存
;   [rsp+16..23] = n1 一時保存

; スタックフレーム:
;   [rsp+0..31]   = CALL 用シャドウスペース (32 bytes, callee が使用可)
;   [rsp+32..39]  = q0 一時保存
;   [rsp+40..47]  = n0 一時保存
;   [rsp+48..55]  = n1 一時保存
;   [rsp+56..63]  = dn 一時保存 (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         ; スタックフレーム確保

    ; パラメータ取得
    ; 呼び出し側の rsp から見た stack args:
    ;   5th arg (dn)    = [original_rsp + 32]
    ;   6th arg (dinv3) = [original_rsp + 40]
    ; 被呼出し側: 8 pushes (64) + sub rsp,FRAME_SIZE (64) = 128 bytes
    ; よって: 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)

    ; 定数ロード
    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 は j = qn downto 0 (qn+1 回) で処理。
    ; j=qn の反復は C++ の qh check と重複するが、
    ; np[nn]=0 sentinel により q_hat ≈ 0 で実質 no-op。
    ; ただし np[nn-2], np[nn-1] を udiv_qrnnd_3by2 結果に上書きするため
    ; dcpi1 等の内部呼び出しでは qp[qn] の保護が必要。

outer_loop:
    ; r15 = current j (qn-1, qn-2, ..., 0)
    ; j < 0 なら終了
    test    r15, r15
    js      done

    ; 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

    ; 特殊ケース: 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

    ; dn = rbp + 2 をスタックに保存 (レジスタ不足のため)
    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 (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 (極めてまれ)
    dec     r11
    xor     ecx, ecx
    clc
    ; addback: dn 回ループ + キャリー伝播 (dn+1 limb 目)
    mov     rcx, [rsp + 56]         ; rcx = dn (ループカウンタ)
    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]            ; CF 非破壊
    dec     rcx                     ; CF を保存 (dec は CF 非破壊)
    jnz     qmax_addback
    ; キャリーを a[j+dn] に伝播 (C++ の 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              ; 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++

    ; 条件付き補正
    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
    ; submul_1(a+j, dp, dn-2, q_hat) を C++ ASM 関数で実行
    mov     [rsp+40], rax            ; n0 保存
    mov     [rsp+48], rcx           ; n1 保存
    ; r11 (q_hat) は volatile なので CALL 前に保存
    mov     [rsp+32], r11           ; q_hat 保存 (q0 は不要になったので再利用)

    ; 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 = 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) — add-back 必要か判定
    ; CF を保存してから 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 (ループカウンタ)
    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 は 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