; Copyright (C) 2026 Kiyotsugu Arai
; SPDX-License-Identifier: LGPL-3.0-or-later
;
; mpn_x64_mod1.asm — dedicated 1-limb mod (preinv) function
;
; Function:
;   mpn_mod_1_preinv_asm(ap, an, d_norm, dinv, shift) -> uint64_t remainder
;
; Use: inside IntPrime::isDivisibleBySmallPrime trial division,
;       fast remainder of multi-limb / single-limb (equivalent to GMP's mpn_mod_1_1).
;
; Algorithm:
;   Returns (Σ ap[i] * 2^{64i}) mod d_actual (d_actual = d_norm >> shift)
;   For each limb, using the Möller-Granlund (2011) preinv quotient estimate,
;   iterate r ← ((r << 64) | n1) mod d_norm.
;   When shift > 0, combine adjacent limbs to form n1 (normalized n1).
;
; Windows x64 calling convention:
;   rcx = ap, rdx = an, r8 = d_norm, r9 = dinv
;   [rsp+40] = shift (uint64_t, low byte is 0..63)
;   Return value: rax = remainder
;   Preserved: rbx, rbp, rdi, rsi, r12-r15
;   Clobbered: rax, rcx, rdx, r8-r11

.code

mpn_mod_1_preinv_asm PROC
    ; Prologue: push 5 callee-save regs (40 bytes), + 8 for 16-byte alignment
    push    rbx
    push    rsi
    push    rdi
    push    r12
    push    r13
    sub     rsp, 8

    ; Stack at this point:
    ;   [rsp+0..7]    - alignment pad
    ;   [rsp+8..47]   - 5 pushed regs (rbx, rsi, rdi, r12, r13)
    ;   [rsp+48..55]  - return addr
    ;   [rsp+56..87]  - Win64 shadow space (caller's home for rcx, rdx, r8, r9)
    ;   [rsp+88..95]  - arg5 (shift)
    SHIFT_OFFSET = 88

    ; Argument setup
    ;   rcx = ap            → rsi
    ;   rdx = an            → r10
    ;   r8  = d_norm        (stay)
    ;   r9  = dinv          → r11 (RAX/RDX are used by MUL, so r9 is not freed)
    ;   shift on stack      → rcx (shift in CL, or rshift in rdi)
    mov     rsi, rcx                   ; rsi = ap
    mov     r10, rdx                   ; r10 = an
    mov     r11, r9                    ; r11 = dinv
    movzx   ecx, byte ptr [rsp + SHIFT_OFFSET] ; cl = shift (0..63)

    ; an == 0 → return 0
    test    r10, r10
    jz      ret_zero

    ; r12 = r (running remainder)
    ; r13 = n1 (current digit)
    ; r8  = d_norm (unchanged)
    ; r11 = dinv (unchanged)
    ; rcx = shift
    ; rdi = rshift = 64 - shift (only used when shift > 0)

    ; shift == 0 → simpler loop
    test    cl, cl
    jnz     shift_nonzero

    ; ─────────────── shift == 0 ───────────────
    xor     r12, r12                   ; r = 0
    ; Used as loop counter rcx, so shift is unneeded (cl=0 so rcx=0)
    mov     rcx, r10                   ; rcx = an (loop counter, downto 0)

shift0_loop:
    ; n1 = ap[rcx-1]
    mov     r13, [rsi + rcx*8 - 8]
    ; mod_step(r, n1, d_norm, dinv) → r
    ; rdx:rax = r * dinv
    mov     rax, r12
    mul     r11                        ; rdx:rax = r12 * r11 (high:low)
    add     rax, r13                   ; q0 = p_lo + n1, CF = carry
    adc     rdx, r12                   ; q1 = p_hi + r + CF
    inc     rdx                        ; q1++ (tentative)
    ; rem = n1 - q1 * d_norm
    mov     rdi, rdx                   ; rdi = q1
    imul    rdi, r8                    ; rdi = q1 * d_norm (low 64 bits)
    mov     r9, r13                    ; r9 = n1
    sub     r9, rdi                    ; r9 = rem = n1 - q1*d_norm
    ; if rem > q0: rem += d_norm
    cmp     r9, rax                    ; rem vs q0
    jbe     shift0_no_corr1
    add     r9, r8                     ; rem += d_norm
shift0_no_corr1:
    ; if rem >= d_norm: rem -= d_norm
    cmp     r9, r8
    jb      shift0_no_corr2
    sub     r9, r8                     ; rem -= d_norm
shift0_no_corr2:
    mov     r12, r9                    ; r = rem

    dec     rcx
    jnz     shift0_loop

    mov     rax, r12
    jmp     epilogue

    ; ─────────────── shift != 0 ───────────────
shift_nonzero:
    ; rdi = rshift = 64 - shift
    mov     rdi, 64
    sub     rdi, rcx
    ; rcx (shift) is preserved; we'll use it as CL for shifts
    ; rdi as CL for rshift (we'll switch via mov)

    ; r = ap[an-1] >> rshift
    mov     r12, [rsi + r10*8 - 8]     ; r12 = ap[an-1]
    push    rcx                         ; save shift
    mov     rcx, rdi                    ; cl = rshift
    shr     r12, cl                     ; r = ap[an-1] >> rshift
    pop     rcx                         ; restore shift in cl

    ; if an == 1: skip loop, just process last limb
    cmp     r10, 1
    je      shift_last

    ; Loop: i = an-1 downto 1
    ; rax (or some reg) = i; we use r10 as i (and decrement)
    ; n1 = (ap[i] << shift) | (ap[i-1] >> rshift)

shift_loop:
    ; i = r10
    mov     r13, [rsi + r10*8 - 8]     ; high = ap[i-1+1] = ap[i] (we treat r10 = i+1 for indexing convenience)
    ; Wait: we want ap[i] where i goes from an-1 down to 1.
    ; r10 currently = an. ap[i] = ap[an-1] when r10 = an.
    ; Use offset: ap[r10-1] = ap[i] where i = r10-1.
    ; So loop: while r10 > 1: i = r10-1; process ap[i] and ap[i-1]; r10--
    ; Currently r10 = an (>= 2 here), so i = an-1.
    ; ap[i] = [rsi + (r10-1)*8] = [rsi + r10*8 - 8]
    ; ap[i-1] = [rsi + (r10-2)*8] = [rsi + r10*8 - 16]

    mov     r13, [rsi + r10*8 - 8]     ; ap[i]
    mov     r9, [rsi + r10*8 - 16]     ; ap[i-1]
    ; n1 = (ap[i] << shift) | (ap[i-1] >> rshift)
    ;   SHLD dst, src, cl: dst = (dst << cl) | (src >> (64-cl))
    ;   So SHLD r13, r9, cl gives r13 = (r13 << cl) | (r9 >> (64-cl))
    shld    r13, r9, cl                ; r13 = n1
    ; mod_step(r, n1, d_norm, dinv)
    mov     rax, r12
    mul     r11                         ; rdx:rax = r * dinv
    add     rax, r13                    ; q0
    adc     rdx, r12                    ; q1 (with carry)
    inc     rdx
    mov     r9, rdx
    imul    r9, r8                      ; q1 * d_norm
    mov     rdx, r13
    sub     rdx, r9                     ; rem = n1 - q1*d_norm
    cmp     rdx, rax
    jbe     shiftN_no_corr1
    add     rdx, r8
shiftN_no_corr1:
    cmp     rdx, r8
    jb      shiftN_no_corr2
    sub     rdx, r8
shiftN_no_corr2:
    mov     r12, rdx                    ; r = rem

    dec     r10
    cmp     r10, 1
    jne     shift_loop

shift_last:
    ; n1 = ap[0] << shift
    mov     r13, [rsi]
    shl     r13, cl
    ; mod_step
    mov     rax, r12
    mul     r11
    add     rax, r13
    adc     rdx, r12
    inc     rdx
    mov     r9, rdx
    imul    r9, r8
    mov     rdx, r13
    sub     rdx, r9
    cmp     rdx, rax
    jbe     last_no_corr1
    add     rdx, r8
last_no_corr1:
    cmp     rdx, r8
    jb      last_no_corr2
    sub     rdx, r8
last_no_corr2:
    mov     r12, rdx

    ; rax = r12 >> shift  (denormalize)
    mov     rax, r12
    shr     rax, cl
    jmp     epilogue

ret_zero:
    xor     rax, rax

epilogue:
    add     rsp, 8
    pop     r13
    pop     r12
    pop     rdi
    pop     rsi
    pop     rbx
    ret
mpn_mod_1_preinv_asm ENDP

END