; 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