// Copyright (C) 2026 Kiyotsugu Arai
// SPDX-License-Identifier: LGPL-3.0-or-later

// IntModular.cpp
// Modular arithmetic implementation

#include "math/core/mp/Int/IntModular.hpp"
#include "math/core/mp/Int/IntOps.hpp"
#include "math/core/mp/Int/IntGCD.hpp"
#include "math/core/mp/Int/IntSpecialStates.hpp"
#include "math/core/mp/Int/MpnOps.hpp"
#include <vector>
#include <cstring>

#if defined(_MSC_VER) && defined(_M_X64)
#include <intrin.h>
#include <immintrin.h>
#endif

// ASM Montgomery function prototypes (declared at global scope)
#ifdef SANGI_INT_HAS_ASM
extern "C" void mpn_mont_mul_mulx(uint64_t* rp, const uint64_t* ap,
                                   const uint64_t* bp, size_t n,
                                   const uint64_t* mp, uint64_t m_inv,
                                   uint64_t* scratch);
extern "C" void mpn_mont_redc_mulx(uint64_t* rp, uint64_t* tp, size_t tn,
                                    const uint64_t* mp, size_t n,
                                    uint64_t m_inv);
extern "C" void mpn_mont_mul_8(uint64_t* rp, const uint64_t* ap,
                                const uint64_t* bp, const uint64_t* mp,
                                uint64_t m_inv);
extern "C" void mpn_mont_sqr_8(uint64_t* rp, const uint64_t* ap,
                                const uint64_t* mp, uint64_t m_inv);
extern "C" void mpn_mont_redc_16(uint64_t* rp, uint64_t* tp,
                                  const uint64_t* mp, uint64_t m_inv);
extern "C" void mpn_mont_redc_32(uint64_t* rp, uint64_t* tp,
                                  const uint64_t* mp, uint64_t m_inv);
extern "C" void mpn_mont_sqr_16(uint64_t* rp, const uint64_t* ap,
                                  const uint64_t* mp, uint64_t m_inv);
extern "C" void mpn_mont_mul_4(uint64_t* rp, const uint64_t* ap,
                                const uint64_t* bp, const uint64_t* mp,
                                uint64_t m_inv);
extern "C" void mpn_mont_sqr_4(uint64_t* rp, const uint64_t* ap,
                                const uint64_t* mp, uint64_t m_inv);
#if defined(_MSC_VER)
extern "C" void mpn_mont_mul_3(uint64_t* rp, const uint64_t* ap,
                                const uint64_t* bp, const uint64_t* mp,
                                uint64_t m_inv);
#endif
extern "C" void mpn_mont_mul_2(uint64_t* rp, const uint64_t* ap,
                                const uint64_t* bp, const uint64_t* mp,
                                uint64_t m_inv);
extern "C" void mpn_mont_sqr_2(uint64_t* rp, const uint64_t* ap,
                                const uint64_t* mp, uint64_t m_inv);
#endif

namespace sangi {

// ===========================================================================
// Montgomery multiplication (internal implementation)
// ===========================================================================
namespace {

// Compute -m^{-1} mod 2^64 (m is odd)
// Newton iteration: x ← x·(2 - m·x) mod 2^64
inline uint64_t mont_neg_inv(uint64_t m0) {
    uint64_t x = 1;
    for (int i = 0; i < 6; ++i)          // converges to 64-bit precision in 6 iterations
        x *= 2 - m0 * x;                 // mod 2^64 is automatic (uint64_t)
    return static_cast<uint64_t>(0) - x;  // -x mod 2^64
}

// Montgomery reduction: T * R^{-1} mod m
// T: 2n words, m: n words, m_inv = -m^{-1} mod 2^64
// Result: r (n words), return: 0 or 1 (r >= m)
inline void mont_redc(uint64_t* r,
                      uint64_t* t, size_t tn,   // T (destroyed, requires 2n+1 words)
                      const uint64_t* m, size_t n,
                      uint64_t m_inv) {
    // CIOS (Coarsely Integrated Operand Scanning) method
    for (size_t i = 0; i < n; ++i) {
        uint64_t q = t[i] * m_inv;   // q = T[i] · (-m^{-1}) mod 2^64
        uint64_t carry = mpn::addmul_1(t + i, m, n, q);  // T += q · m · B^i
        // propagate carry to higher words
        for (size_t j = i + n; carry && j < tn; ++j) {
            uint64_t sum = t[j] + carry;
            carry = (sum < t[j]) ? 1 : 0;
            t[j] = sum;
        }
    }
    // Result is T[n..2n-1]
    std::memcpy(r, t + n, n * sizeof(uint64_t));
    // if r >= m then r -= m
    if (mpn::cmp(r, n, m, n) >= 0) {
        mpn::sub(r, r, n, m, n);
    }
}

// Compute the 2-word inverse: mip * m ≡ -1 mod B^2
// mip[0] = m_inv (existing), compute mip[1]
// Algorithm equivalent to GMP sec_powm.c
inline uint64_t mont_neg_inv2(uint64_t m0, uint64_t m1, uint64_t m_inv) {
    // since m_inv * m[0] ≡ -1 mod B, lo(m_inv * m[0]) = 0xFFFF...FFFF
    // h = hi(m_inv * m[0])
    uint64_t h;
#if defined(_MSC_VER) && defined(_M_X64)
    _umul128(m_inv, m0, &h);
#else
    unsigned __int128 p = static_cast<unsigned __int128>(m_inv) * m0;
    h = static_cast<uint64_t>(p >> 64);
#endif
    // t = h + m_inv * m[1]  (word 1 of m_inv * m, mod B)
    uint64_t t = h + m_inv * m1;
    // mip[1] = (t + 1) * m_inv  (mod B)
    return (t + 1) * m_inv;
}

// Montgomery reduction k=2: GMP mpn_redc_2 algorithm
// Processes REDC iterations two words at a time, halving the loop count
// tp[0..2n] is destroyed (requires 2n+1 words)
// n must be even (for odd n, the caller handles the leading 1 iteration)
inline void mont_redc_2(uint64_t* rp, uint64_t* tp,
                        const uint64_t* mp, size_t n,
                        uint64_t mip0, uint64_t mip1) {
    // GMP redc_2 algorithm:
    // for each pair (i, i+1):
    //   q = umul2low(mip, tp[i..i+1])  — 2-word quotient
    //   save tp[n], addmul_2(tp, mp, n, q), bookkeeping
    // final: add_n(rp, tp+n, tp, n)

    uint64_t* up = tp;

    // if n is odd: handle the leading 1 iteration normally
    if (n & 1) {
        uint64_t q = up[0] * mip0;
        uint64_t cy = mpn::addmul_1(up, mp, n, q);
        up[n] += cy;
        up++;
    }

    for (size_t j = (n & 1) ? (n - 1) : n; j >= 2; j -= 2) {
        // 2-word quotient: (q1, q0) = umul2low(mip, up[0..1])
        uint64_t q0, q1;
        q0 = up[0] * mip0;
        // q1 = hi(mip0 * up[0]) + mip0 * up[1] + mip1 * up[0]
        uint64_t h;
#if defined(_MSC_VER) && defined(_M_X64)
        _umul128(mip0, up[0], &h);
#else
        unsigned __int128 p = static_cast<unsigned __int128>(mip0) * up[0];
        h = static_cast<uint64_t>(p >> 64);
#endif
        q1 = h + mip0 * up[1] + mip1 * up[0];

        // addmul_2: up[0..n] += mp[0..n-1] * (q0 + q1*B)
        uint64_t upn = up[n];   // save (will be overwritten)
        up[n] = mpn::addmul_1(up, mp, n, q0);   // carry → up[n]
        uint64_t cy = mpn::addmul_1(up + 1, mp, n, q1);  // carry returned
        // GMP bookkeeping:
        up[1] = cy;       // addmul_2 return value
        up[0] = up[n];    // accumulated carry → position 0
        up[n] = upn;      // restore original up[n]
        up += 2;
    }

    // add the upper half + lower half (carry information)
    uint64_t cy = mpn::add(rp, up, n, up - n, n);

    // if cy > 0 or rp >= mp then rp -= mp
    if (cy != 0 || mpn::cmp(rp, n, mp, n) >= 0) {
        mpn::sub(rp, rp, n, mp, n);
    }
}

// Threshold for redc_n: use redc_n (O(M(n))) at or above this size
// since mul_low is now Karatsuba-based, the threshold can be lowered
static constexpr size_t REDC_N_THRESHOLD = 48;

// Compute the Hensel inverse: result * m ≡ -1 (mod B^n)
// Newton iteration: x ← x * (2 + m * x) mod B^{2k}
// (m_inv = -m^{-1} mod B is precomputed)
// scratch: requires 3*n words
inline void hensel_inverse(uint64_t* result, const uint64_t* m, size_t n,
                           uint64_t m_inv, uint64_t* scratch) {
    // since m_inv is -m[0]^{-1} mod B, result[0] = m_inv
    // this satisfies m_inv * m[0] ≡ -1 (mod B)
    std::memset(result, 0, n * sizeof(uint64_t));
    result[0] = m_inv;

    // scratch layout: shares t[0..2k-1], p[0..2k-1], x_new[0..k-1]
    // each stage needs at most 3*n words

    size_t cur = 1;  // current number of significant words
    while (cur < n) {
        size_t next = std::min(2 * cur, n);
        // t = lower next words of m[0..next-1] * result[0..cur-1]
        // = m * x mod B^next
        uint64_t* t = scratch;          // uses next words
        uint64_t* p = scratch + next;   // uses next words
        // compute the lower next words with basecase multiplication (upper not needed)
        // t[0..next-1] = lower next words of m[0..next-1] * result[0..cur-1]
        std::memset(t, 0, next * sizeof(uint64_t));
        for (size_t i = 0; i < cur; ++i) {
            if (result[i] == 0) continue;
            size_t jmax = std::min(next - i, next);
            uint64_t carry = mpn::addmul_1(t + i, m, jmax, result[i]);
            // carry propagates to word next and beyond, but ignored since mod B^next
            (void)carry;
        }
        // t = m * x mod B^next
        // compute the lower next words of 2 + t: t[0] is -1 mod B (m_inv*m[0] = -1 mod B)
        // so 2 + t: t[0] += 2 → t[0] = 1 (mod B) (no carry)
        // in fact t = m*x ≡ -1 (mod B^cur), so all words are 0xFFF...FFF
        // 2 + t mod B^next = 1 (lower cur words) + (2+t)[cur..next-1]*B^cur
        // Newton: x_new = x * (2 + m*x) mod B^next
        // here we compute 2 + t
        t[0] += 2;  // 2 + m*x (mod B^next); the lower cur words become 0...01
        // t[0..cur-1] = 1, 0, 0, ..., 0 (carry propagated)
        // in fact t[0] = FFFFFFFFFFFFFFFFh + 2 = 1 (carry=1)
        // t[1] = FFFFFFFFFFFFFFFFh + 1 = 0 (carry=1) ... t[cur-1]+1=0, carry=1
        // t[cur] += 1 (carry from below)
        // this is tedious, so a different approach: compute t = 2 + m*x mod B^next directly
        // the lower cur words are always [1, 0, 0, ..., 0] (a property of Newton)
        // → only the upper (next-cur) words of x_new = x * t mod B^next are new
        // → result[cur..next-1] = lower (next-cur) words of (x * t[cur..next-1])
        //   + the upper part of x * [1,0,...,0] (= result[0..cur-1] unchanged)
        // concisely: result[0..next-1] = result[0..cur-1] * t[0..next-1] mod B^next

        // p = lower next words of result[0..cur-1] * t[0..next-1]
        std::memset(p, 0, next * sizeof(uint64_t));
        for (size_t i = 0; i < cur; ++i) {
            if (result[i] == 0) continue;
            size_t jmax = std::min(next - i, next);
            mpn::addmul_1(p + i, t, jmax, result[i]);
        }
        // result[0..next-1] = p[0..next-1]
        std::memcpy(result, p, next * sizeof(uint64_t));
        cur = next;
    }
}

// mul_low: (a * b) mod B^n — computes only the lower n words
// Karatsuba version: omits a1*b1, and recurses on a0*b1 + a1*b0 via mul_low(h)
// below the basecase threshold, uses an addmul_1 loop
static constexpr size_t MUL_LOW_THRESHOLD = 24;  // Karatsuba switch threshold

static void mul_low_n(uint64_t* rp, const uint64_t* ap, const uint64_t* bp,
                      size_t n, uint64_t* scratch) {
    if (n < MUL_LOW_THRESHOLD) {
        // basecase: addmul_1 loop (discards the upper words)
        std::memset(rp, 0, n * sizeof(uint64_t));
        for (size_t i = 0; i < n; ++i) {
            if (ap[i] == 0) continue;
            size_t jmax = n - i;
            mpn::addmul_1(rp + i, bp, jmax, ap[i]);
        }
        return;
    }

    // Karatsuba mul_low (combines the cross terms into a single multiplication)
    size_t h = n / 2;
    size_t l = n - h;  // l >= h

    // a = a1·B^h + a0, b = b1·B^h + b0
    // (a*b) mod B^n = a0*b0 + (a0*b1 + a1*b0)·B^h mod B^n
    //
    // Karatsuba trick:
    //   a0*b1 + a1*b0 = (a0+a1)*(b0+b1) - a0*b0 - a1*b1
    //   a1*b1 contributes only at B^{2h} or above, so it is not needed mod B^n
    //   → only the lower l words of cross = (a0+a1)*(b0+b1) - a0*b0 are needed
    //   → reduces 2 mul_low calls → 1 full multiplication

    // scratch layout: s[h+1] + t[h+1] + p[2(h+1)] + mul_scratch
    uint64_t* s = scratch;                          // h+1 words
    uint64_t* t = scratch + h + 1;                  // h+1 words
    uint64_t* p = scratch + 2 * (h + 1);            // 2*(h+1) words
    uint64_t* rec_scratch = scratch + 4 * (h + 1);  // remainder

    // Step 1: rp[0..2h-1] = a0 * b0 (full multiplication)
    size_t mul_s = mpn::multiply_scratch_size(h, h);
    if (mul_s > 0)
        mpn::multiply(rp, ap, h, bp, h, rec_scratch);
    else
        mpn::multiply(rp, ap, h, bp, h, nullptr);

    // Step 2: s = a0 + a1, t = b0 + b1 (h+1 words each)
    s[h] = mpn::add(s, ap, h, ap + h, h);
    t[h] = mpn::add(t, bp, h, bp + h, h);
    size_t sn = h + (s[h] ? 1 : 0);
    size_t tn = h + (t[h] ? 1 : 0);

    // Step 3: p = s * t (full multiplication)
    size_t mul_st = mpn::multiply_scratch_size(sn, tn);
    if (sn >= tn) {
        if (mul_st > 0) mpn::multiply(p, s, sn, t, tn, rec_scratch);
        else mpn::multiply(p, s, sn, t, tn, nullptr);
    } else {
        if (mul_st > 0) mpn::multiply(p, t, tn, s, sn, rec_scratch);
        else mpn::multiply(p, t, tn, s, sn, nullptr);
    }
    size_t pn = mpn::normalized_size(p, sn + tn);

    // Step 4: lower l words of cross = p - a0*b0
    // p -= rp[0..2h-1] (a0*b0)
    size_t sub_n = std::min(pn, 2 * h);
    mpn::sub(p, p, pn, rp, sub_n);
    // add the lower l words of p to rp[h..n-1]
    size_t add_n = std::min(pn, l);
    mpn::add(rp + h, rp + h, l, p, add_n);
}

// redc_n: O(M(n)) Montgomery reduction
// Algorithm equivalent to GMP mpn_redc_n:
//   Q = T[0..n-1] * m_inv_n[0..n-1] mod B^n  (lower n words only)
//   P = Q * m (full multiplication)
//   result = T[n..2n-1] + P[n..2n-1]
//   conditional subtraction
//
// scratch: multiply_scratch_size(n,n) words (for multiply)
// tp is destroyed
// scratch size required by mul_low (accounting for recursion)
static size_t mul_low_scratch_size(size_t n) {
    if (n < MUL_LOW_THRESHOLD) return 0;
    size_t h = n / 2;
    size_t l = n - h;
    // t[l] + recursive (multiply scratch or mul_low scratch)
    size_t mul_s = mpn::multiply_scratch_size(h, h);
    size_t rec = mul_low_scratch_size(l);
    return l + std::max(mul_s, rec);
}

inline void mont_redc_n(uint64_t* rp, uint64_t* tp,
                        const uint64_t* mp, size_t n,
                        const uint64_t* m_inv_n,
                        uint64_t* scratch) {
    // Step 1: Q = T_low * m_inv_n mod B^n (mul_low yields only the lower n words)
    uint64_t* Q = static_cast<uint64_t*>(_alloca(n * sizeof(uint64_t)));
    // use scratch as mul_low's working area
    mul_low_n(Q, tp, m_inv_n, n, scratch);

    // Step 2: P = Q * m (full multiplication, 2n words)
    uint64_t* P = static_cast<uint64_t*>(_alloca(2 * n * sizeof(uint64_t)));
    size_t ms = mpn::multiply_scratch_size(n, n);
    if (ms > 0) {
        mpn::multiply(P, Q, n, mp, n, scratch);
    } else {
        mpn::multiply(P, Q, n, mp, n, nullptr);
    }

    // Step 3: rp = tp[n..2n-1] + P[n..2n-1]
    uint64_t cy = mpn::add(rp, tp + n, n, P + n, n);

    // Step 4: conditional subtraction
    if (cy != 0 || mpn::cmp(rp, n, mp, n) >= 0) {
        mpn::sub(rp, rp, n, mp, n);
    }
}

// FIOS (Finely Integrated Operand Scanning) Montgomery multiplication
// Fuses multiplication and reduction into a single pass. Halves memory accesses.
// May be faster than CIOS at large sizes n≥16.
// No scratch needed (only an internal n+2 word buffer).
// FIOS Montgomery multiplication: works on all platforms
// GCC: optimized with __uint128_t, MSVC: _umul128 intrinsic
inline void mont_mul_fios(uint64_t* r,
                          const uint64_t* a, const uint64_t* b,
                          const uint64_t* m, size_t n,
                          uint64_t m_inv) {
    // t[0..n+1]: shift accumulator
    uint64_t* t = static_cast<uint64_t*>(_alloca((n + 2) * sizeof(uint64_t)));
    std::memset(t, 0, (n + 2) * sizeof(uint64_t));

#if defined(__SIZEOF_INT128__)
    // GCC/Clang: concise and fast with __uint128_t
    for (size_t i = 0; i < n; ++i) {
        uint64_t bi = b[i];

        // Word 0
        __uint128_t uv = (__uint128_t)a[0] * bi + t[0];
        uint64_t S = (uint64_t)uv;
        __uint128_t Ca = uv >> 64;

        uint64_t q = S * m_inv;

        uv = (__uint128_t)m[0] * q + S;
        __uint128_t Cm = uv >> 64;

        // Words 1..n-1
        for (size_t j = 1; j < n; ++j) {
            uv = (__uint128_t)a[j] * bi + t[j] + Ca;
            S = (uint64_t)uv;
            Ca = uv >> 64;

            uv = (__uint128_t)m[j] * q + S + Cm;
            t[j - 1] = (uint64_t)uv;
            Cm = uv >> 64;
        }

        uv = (__uint128_t)t[n] + Ca + Cm;
        t[n - 1] = (uint64_t)uv;
        t[n] = (uint64_t)(uv >> 64) + t[n + 1];
        t[n + 1] = 0;
    }
#else
    // MSVC: _umul128 intrinsic
    for (size_t i = 0; i < n; ++i) {
        uint64_t bi = b[i];
        uint64_t hi_a, lo_a, hi_m, lo_m;

        lo_a = _umul128(a[0], bi, &hi_a);
        uint64_t S = t[0] + lo_a;
        uint64_t Ca = hi_a + (S < t[0] ? 1ULL : 0);
        unsigned Ca_hi = 0;

        uint64_t q = S * m_inv;

        lo_m = _umul128(m[0], q, &hi_m);
        uint64_t tmp = S + lo_m;
        uint64_t Cm = hi_m + (tmp < S ? 1ULL : 0);
        unsigned Cm_hi = 0;

        for (size_t j = 1; j < n; ++j) {
            lo_a = _umul128(a[j], bi, &hi_a);
            uint64_t s1 = t[j] + lo_a;
            unsigned c1 = (s1 < t[j]) ? 1u : 0u;
            S = s1 + Ca;
            c1 += (S < s1) ? 1u : 0u;
            uint64_t new_Ca = hi_a + c1;
            unsigned new_Ca_hi = (new_Ca < c1) ? 1u : 0u;
            new_Ca += Ca_hi;
            new_Ca_hi += (new_Ca < Ca_hi) ? 1u : 0u;

            lo_m = _umul128(m[j], q, &hi_m);
            uint64_t s2 = S + lo_m;
            unsigned c2 = (s2 < S) ? 1u : 0u;
            uint64_t res = s2 + Cm;
            c2 += (res < s2) ? 1u : 0u;
            uint64_t new_Cm = hi_m + c2;
            unsigned new_Cm_hi = (new_Cm < c2) ? 1u : 0u;
            new_Cm += Cm_hi;
            new_Cm_hi += (new_Cm < Cm_hi) ? 1u : 0u;

            t[j - 1] = res;
            Ca = new_Ca;
            Ca_hi = new_Ca_hi;
            Cm = new_Cm;
            Cm_hi = new_Cm_hi;
        }

        uint64_t f = t[n] + Ca;
        unsigned fc = (f < t[n]) ? 1u : 0u;
        f += Cm;
        fc += (f < Cm) ? 1u : 0u;
        t[n - 1] = f;
        t[n] = static_cast<uint64_t>(fc + Ca_hi + Cm_hi) + t[n + 1];
        t[n + 1] = 0;
    }
#endif

    // conditional subtraction
    if (t[n] != 0 || mpn::cmp(t, n, m, n) >= 0) {
        mpn::sub(r, t, n, m, n);
    } else {
        std::memcpy(r, t, n * sizeof(uint64_t));
    }
}

// Montgomery multiplication: (a * b * R^{-1}) mod m
// a, b: n words (Montgomery form), m: n words
// Result: r (n words)
inline void mont_mul(uint64_t* r,
                     const uint64_t* a, const uint64_t* b,
                     const uint64_t* m, size_t n,
                     uint64_t m_inv, uint64_t mip1,
                     uint64_t* scratch,                    // 2n+1 words (product buffer)
                     uint64_t* mul_scratch = nullptr,      // multiply_scratch_size(n,n) words
                     const uint64_t* m_inv_n = nullptr) {  // n-word inverse for redc_n (nullable)
#ifdef SANGI_INT_HAS_ASM
    if (mpn::detail::has_bmi2_adx()) {
        if (n == 2) {
            mpn_mont_mul_2(r, a, b, m, m_inv);
            return;
        }
#if defined(_MSC_VER)
        if (n == 3) {
            // 192-bit specialization (avoids the MULX/ADCX/ADOX performance cliff of MSVC FIOS)
            // Windows MASM only, since the SysV ABI version is not implemented
            mpn_mont_mul_3(r, a, b, m, m_inv);
            return;
        }
#endif
        if (n == 4) {
            mpn_mont_mul_4(r, a, b, m, m_inv);
            return;
        }
        if (n == 8) {
            mpn_mont_mul_8(r, a, b, m, m_inv);
            return;
        }
        // n=5..12: FIOS is faster than SOS (basecase+REDC)
        // (fusing multiplication and reduction halves memory accesses)
        if (n <= 12 && n >= 5) {
            mont_mul_fios(r, a, b, m, n, m_inv);
            return;
        }
        if (n == 16) {
            // SOS separated method: basecase multiplication + specialized REDC
            mpn::mul_basecase(scratch, a, n, b, n);
            scratch[2 * n] = 0;
            mpn_mont_redc_16(r, scratch, m, m_inv);
            return;
        }
        if (n == 32) {
            // SOS separated method: basecase multiplication + specialized REDC
            mpn::mul_basecase(scratch, a, n, b, n);
            scratch[2 * n] = 0;
            mpn_mont_redc_32(r, scratch, m, m_inv);
            return;
        }
        // generic path: SOS separated method + REDC
        if (mul_scratch && n >= 64) {
            mpn::multiply(scratch, a, n, b, n, mul_scratch);
        } else {
            mpn::mul_basecase(scratch, a, n, b, n);
        }
        scratch[2 * n] = 0;
        if (n >= REDC_N_THRESHOLD && m_inv_n && mul_scratch) {
            mont_redc_n(r, scratch, m, n, m_inv_n, mul_scratch);
            return;
        }
        mont_redc_2(r, scratch, m, n, m_inv, mip1);
        return;
    }
#endif
    // C++ fallback: FIOS (n≥3) + n=2 specialization + CIOS (n=1)
#if defined(__SIZEOF_INT128__)
    // n=2 specialization: __uint128_t-based CIOS (2-limb Montgomery multiplication)
    // Since n=2 (128-bit) appears frequently in Miller-Rabin small-prime tests,
    // this uses an integrated loop without function calls for speed.
    if (n == 2) {
        // SOS: 2x2 multiplication (4-limb result)
        __uint128_t ab00 = (__uint128_t)a[0] * b[0];
        __uint128_t ab01 = (__uint128_t)a[0] * b[1];
        __uint128_t ab10 = (__uint128_t)a[1] * b[0];
        __uint128_t ab11 = (__uint128_t)a[1] * b[1];

        // t[0..4]: product accumulation
        uint64_t t0 = (uint64_t)ab00;
        __uint128_t carry = (ab00 >> 64) + (uint64_t)ab01 + (uint64_t)ab10;
        uint64_t t1 = (uint64_t)carry;
        carry = (carry >> 64) + (ab01 >> 64) + (ab10 >> 64) + (uint64_t)ab11;
        uint64_t t2 = (uint64_t)carry;
        uint64_t t3 = (uint64_t)(carry >> 64) + (uint64_t)(ab11 >> 64);

        // REDC iteration 0
        uint64_t q0 = t0 * m_inv;
        __uint128_t qm0 = (__uint128_t)q0 * m[0];
        __uint128_t qm1 = (__uint128_t)q0 * m[1];
        carry = (__uint128_t)t0 + (uint64_t)qm0;
        // carry>>64 == 1 or 0 (carry-out of t0 + qm0_lo)
        carry = (__uint128_t)t1 + (qm0 >> 64) + (uint64_t)qm1 + (carry >> 64);
        t0 = (uint64_t)carry;
        carry = (__uint128_t)t2 + (qm1 >> 64) + (carry >> 64);
        t1 = (uint64_t)carry;
        t2 = t3 + (uint64_t)(carry >> 64);

        // REDC iteration 1
        uint64_t q1 = t0 * m_inv;
        qm0 = (__uint128_t)q1 * m[0];
        qm1 = (__uint128_t)q1 * m[1];
        carry = (__uint128_t)t0 + (uint64_t)qm0;
        carry = (__uint128_t)t1 + (qm0 >> 64) + (uint64_t)qm1 + (carry >> 64);
        r[0] = (uint64_t)carry;
        carry = (__uint128_t)t2 + (qm1 >> 64) + (carry >> 64);
        r[1] = (uint64_t)carry;
        uint64_t top = (uint64_t)(carry >> 64);

        // conditional subtraction
        if (top != 0 || (r[1] > m[1] || (r[1] == m[1] && r[0] >= m[0]))) {
            __uint128_t sub = (__uint128_t)r[0] - m[0];
            r[0] = (uint64_t)sub;
            r[1] -= m[1] + (uint64_t)(sub >> 127);
        }
        return;
    }
#endif
    if (n >= 3) {
        mont_mul_fios(r, a, b, m, n, m_inv);
        return;
    }
    // n=1: simple CIOS
    std::memset(scratch, 0, (2 * n + 1) * sizeof(uint64_t));
    for (size_t i = 0; i < n; ++i) {
        uint64_t carry = mpn::addmul_1(scratch + i, a, n, b[i]);
        scratch[i + n] += carry;
        if (scratch[i + n] < carry) {
            for (size_t j = i + n + 1; j <= 2 * n; ++j) {
                ++scratch[j];
                if (scratch[j] != 0) break;
            }
        }
        uint64_t q = scratch[i] * m_inv;
        carry = mpn::addmul_1(scratch + i, m, n, q);
        scratch[i + n] += carry;
        if (scratch[i + n] < carry) {
            for (size_t j = i + n + 1; j <= 2 * n; ++j) {
                ++scratch[j];
                if (scratch[j] != 0) break;
            }
        }
    }
    if (scratch[2 * n] != 0 || mpn::cmp(scratch + n, n, m, n) >= 0) {
        mpn::sub(r, scratch + n, n, m, n);
    } else {
        std::memcpy(r, scratch + n, n * sizeof(uint64_t));
    }
}

// Montgomery squaring: (a^2 * R^{-1}) mod m
// Accelerated with mpn::square (symmetry: 1.5n² MULX) + mont_redc
//
// Verified: tried an n=8 CIOS unrolled version / CIOS mont_mul(a,a), but on
// Zen 3 the MULX/ADCX/ADOX port contention (all on port 0/3) is the bottleneck,
// so neither loop removal nor method change helped. On Intel Skylake+, ADCX/ADOX
// are spread across separate ports, so different results are expected.
inline void mont_sqr(uint64_t* r,
                     const uint64_t* a,
                     const uint64_t* m, size_t n,
                     uint64_t m_inv, uint64_t mip1,
                     uint64_t* scratch,          // 2n+1 words
                     uint64_t* sq_scratch,       // square_scratch_size(n) words
                     const uint64_t* m_inv_n = nullptr,   // n-word inverse for redc_n (nullable)
                     uint64_t* mul_scratch = nullptr) {    // multiply scratch for redc_n
#ifdef SANGI_INT_HAS_ASM
    if (mpn::detail::has_bmi2_adx()) {
        if (n == 2) {
            mpn_mont_sqr_2(r, a, m, m_inv);
            return;
        }
        if (n == 4) {
            mpn_mont_sqr_4(r, a, m, m_inv);
            return;
        }
        if (n == 8) {
            mpn_mont_sqr_8(r, a, m, m_inv);
            return;
        }
        if (n == 16) {
            mpn_mont_sqr_16(r, a, m, m_inv);
            return;
        }
        // n ≤ 12: FIOS is faster than SOS (square+REDC)
        if (n <= 12 && n >= 3) {
            mont_mul_fios(r, a, a, m, n, m_inv);
            return;
        }
        mpn::square(scratch, a, n, sq_scratch);
        scratch[2 * n] = 0;
        if (n >= REDC_N_THRESHOLD && m_inv_n && mul_scratch) {
            mont_redc_n(r, scratch, m, n, m_inv_n, mul_scratch);
            return;
        }
        mont_redc_2(r, scratch, m, n, m_inv, mip1);
        return;
    }
#endif
    // C++ fallback
#if defined(__SIZEOF_INT128__)
    if (n == 2) {
        // reuse the n=2 specialization of mont_mul (a*a)
        mont_mul(r, a, a, m, n, m_inv, mip1, scratch);
        return;
    }
#endif
    if (n >= 3) {
        mont_mul_fios(r, a, a, m, n, m_inv);
        return;
    }
    mpn::square(scratch, a, n, sq_scratch);
    scratch[2 * n] = 0;
    mont_redc(r, scratch, 2 * n + 1, m, n, m_inv);
}

// Conversion to Montgomery form: aR mod m
// R = 2^(64*n), compute aR mod m
// Method: shift a left by n words and take the remainder modulo m
inline void mont_encode(uint64_t* r,
                        const uint64_t* a, size_t an,
                        const uint64_t* m, size_t n) {
    // divide aR = a << (64*n) by m
    // build a 2n-word dividend and divide
    Int shifted;
    // convert a to Int
    std::vector<uint64_t> awords(a, a + an);
    while (!awords.empty() && awords.back() == 0) awords.pop_back();
    if (awords.empty()) {
        std::memset(r, 0, n * sizeof(uint64_t));
        return;
    }
    Int aInt = Int::fromRawWords(awords, 1);
    // aR = a << (64*n)
    IntOps::leftShift(aInt, static_cast<int>(64 * n));
    // aR mod m
    Int mInt = Int::fromRawWords(std::vector<uint64_t>(m, m + n), 1);
    Int result = aInt % mInt;
    // pad the result to n words
    std::memset(r, 0, n * sizeof(uint64_t));
    size_t rn = result.size();
    if (rn > 0) {
        std::memcpy(r, result.words().data(), std::min(rn, n) * sizeof(uint64_t));
    }
}

// ===========================================================================
// AVX2 radix-2^29 Montgomery (SOS: multiplication + REDC separated)
// ===========================================================================
// Uses VPMULUDQ (4×32→64 bit parallel multiplication) with carry-free accumulation.
// A 29-bit limb product is 58 bit → can accumulate ~64 times in a 64-bit accumulator.
// At 512-bit (n29=18) no carry propagation is needed; even at 4096-bit (n29=142),
// only carry propagation every ~30 iterations.

#if defined(_MSC_VER) && defined(_M_X64)
namespace {
namespace r29 {

constexpr int BITS = 29;
constexpr uint64_t MASK = (1ULL << BITS) - 1;

// Number of safe accumulations (before 64-bit overflow):
// each accumulation is at most (2^29-1)^2 ≈ 2^58, floor(2^64 / 2^58) = 64
// CIOS does 2 addmul per iteration → safe iteration count = 30
constexpr size_t SAFE_ITERS = 30;

inline bool has_avx2() {
    static const bool result = []() {
        int info[4];
        __cpuidex(info, 7, 0);
        return (info[1] >> 5) & 1;  // EBX bit 5 = AVX2
    }();
    return result;
}

inline bool has_avx512_ifma() {
    static const bool result = []() {
        int info[4];
        __cpuidex(info, 7, 0);
        bool avx512f = (info[1] >> 16) & 1;  // EBX bit 16
        bool ifma    = (info[1] >> 21) & 1;  // EBX bit 21
        return avx512f && ifma;
    }();
    return result;
}

// 64-bit limb → radix-2^29 limb conversion
// dst must hold at least (n64 * 64 + 28) / 29 + 1 elements
inline size_t to_r29(uint64_t* dst, const uint64_t* src, size_t n64) {
    size_t total_bits = n64 * 64;
    size_t n29 = (total_bits + BITS - 1) / BITS;
    const uint8_t* bytes = reinterpret_cast<const uint8_t*>(src);
    size_t total_bytes = n64 * 8;

    for (size_t i = 0; i < n29; ++i) {
        size_t bit_start  = i * BITS;
        size_t byte_start = bit_start / 8;
        size_t bit_off    = bit_start % 8;
        // read 5 bytes (covers up to 40 bits)
        uint64_t val = 0;
        for (size_t b = 0; b < 5 && byte_start + b < total_bytes; ++b)
            val |= static_cast<uint64_t>(bytes[byte_start + b]) << (b * 8);
        dst[i] = (val >> bit_off) & MASK;
    }
    return n29;
}

// radix-2^29 → 64-bit limb conversion
inline void from_r29(uint64_t* dst, const uint64_t* src, size_t n29, size_t n64) {
    std::memset(dst, 0, n64 * 8);
    uint8_t* bytes = reinterpret_cast<uint8_t*>(dst);
    size_t total_bytes = n64 * 8;

    for (size_t i = 0; i < n29; ++i) {
        size_t bit_start  = i * BITS;
        size_t byte_start = bit_start / 8;
        size_t bit_off    = bit_start % 8;
        uint64_t val = (src[i] & MASK) << bit_off;
        for (size_t b = 0; b < 5 && byte_start + b < total_bytes; ++b)
            bytes[byte_start + b] |= static_cast<uint8_t>(val >> (b * 8));
    }
}

// carry propagation (radix-2^29)
inline void carry_prop(uint64_t* t, size_t n) {
    uint64_t carry = 0;
    for (size_t i = 0; i < n; ++i) {
        t[i] += carry;
        carry = t[i] >> BITS;
        t[i] &= MASK;
    }
}

// -m29[0]^{-1} mod 2^29
inline uint64_t neg_inv_r29(uint64_t m0) {
    uint64_t x = 1;
    for (int i = 0; i < 5; ++i)
        x = (x * (2 - ((m0 * x) & MASK))) & MASK;
    return (MASK + 1 - x) & MASK;
}

// AVX2 vectorized addmul_1: t[0..n-1] += a[0..n-1] * b
// Preconditions: a[i], b are 29-bit or less, and t[i] fits in 64 bits
inline void addmul_1(uint64_t* t, const uint64_t* a, size_t n, uint64_t b) {
    __m256i vb = _mm256_set1_epi64x(static_cast<int64_t>(b));
    size_t i = 0;
    for (; i + 4 <= n; i += 4) {
        __m256i va = _mm256_loadu_si256(reinterpret_cast<const __m256i*>(a + i));
        __m256i vt = _mm256_loadu_si256(reinterpret_cast<const __m256i*>(t + i));
        vt = _mm256_add_epi64(vt, _mm256_mul_epu32(va, vb));
        _mm256_storeu_si256(reinterpret_cast<__m256i*>(t + i), vt);
    }
    for (; i < n; ++i)
        t[i] += static_cast<uint32_t>(a[i]) * static_cast<uint32_t>(b);
}

// r29 domain REDC: T[0..2n29+1] → r29_out[0..n29-1]
// T is destroyed
inline void redc_r29(uint64_t* r29_out, uint64_t* T,
                     const uint64_t* m29, size_t n29,
                     uint64_t m_inv29) {
    for (size_t i = 0; i < n29; ++i) {
        uint64_t q = (T[i] * m_inv29) & MASK;
        addmul_1(T + i, m29, n29, q);
        uint64_t carry = T[i] >> BITS;
        T[i] = 0;
        T[i + 1] += carry;
        if (i > 0 && (i % SAFE_ITERS == 0))
            carry_prop(T + i + 1, n29 - 1);
    }
    carry_prop(T + n29, n29 + 1);

    // conditional subtraction
    int cmp = 0;
    for (int j = static_cast<int>(n29) - 1; j >= 0; --j) {
        if (T[n29 + j] > m29[j]) { cmp = 1; break; }
        if (T[n29 + j] < m29[j]) { cmp = -1; break; }
    }
    if (T[2 * n29] != 0 || cmp >= 0) {
        uint64_t borrow = 0;
        for (size_t j = 0; j < n29; ++j) {
            int64_t diff = static_cast<int64_t>(T[n29 + j])
                         - static_cast<int64_t>(m29[j])
                         - static_cast<int64_t>(borrow);
            borrow = (diff < 0) ? 1 : 0;
            T[n29 + j] = static_cast<uint64_t>(diff) & MASK;
        }
    }
    std::memcpy(r29_out, T + n29, n29 * sizeof(uint64_t));
}

// r29 domain Montgomery multiplication: r29_out = r29_a * r29_b * R_r29^{-1} mod m
// All inputs/outputs radix-2^29 (R_r29 = 2^(29*n29))
// T: working buffer of 2*n29+2 words or more
void mont_mul_r29(uint64_t* r29_out,
                  const uint64_t* r29_a, const uint64_t* r29_b,
                  const uint64_t* m29, size_t n29,
                  uint64_t m_inv29,
                  uint64_t* T) {
    size_t T_len = 2 * n29 + 2;
    std::memset(T, 0, T_len * sizeof(uint64_t));
    for (size_t i = 0; i < n29; ++i)
        addmul_1(T + i, r29_a, n29, r29_b[i]);
    carry_prop(T, 2 * n29 + 1);
    redc_r29(r29_out, T, m29, n29, m_inv29);
}

// r29 domain Montgomery squaring: r29_out = r29_a² * R_r29^{-1} mod m
// Uses symmetry: compute the off-diagonal once and double it
// T: working buffer of 2*n29+2 words or more
void mont_sqr_r29(uint64_t* r29_out,
                  const uint64_t* r29_a,
                  const uint64_t* m29, size_t n29,
                  uint64_t m_inv29,
                  uint64_t* T) {
    size_t T_len = 2 * n29 + 2;
    std::memset(T, 0, T_len * sizeof(uint64_t));
    for (size_t i = 0; i < n29; ++i)
        addmul_1(T + 2 * i + 1, r29_a + i + 1, n29 - i - 1, r29_a[i]);
    for (size_t i = 0; i < 2 * n29; ++i)
        T[i] <<= 1;
    for (size_t i = 0; i < n29; ++i)
        T[2 * i] += static_cast<uint32_t>(r29_a[i]) * static_cast<uint32_t>(r29_a[i]);
    carry_prop(T, 2 * n29 + 1);
    redc_r29(r29_out, T, m29, n29, m_inv29);
}

// Montgomery encode (r29 domain): a * R_r29 mod m → r29 form
// R_r29 = 2^(29 * n29), a is normalized (0 <= a < m)
void mont_encode_r29(uint64_t* r29_out, size_t n29,
                     const uint64_t* a, size_t an,
                     const uint64_t* m64, size_t n64) {
    std::vector<uint64_t> awords(a, a + an);
    while (!awords.empty() && awords.back() == 0) awords.pop_back();
    if (awords.empty()) {
        std::memset(r29_out, 0, n29 * sizeof(uint64_t));
        return;
    }
    Int aInt = Int::fromRawWords(awords, 1);
    IntOps::leftShift(aInt, static_cast<int>(29 * n29));
    Int mInt = Int::fromRawWords(std::vector<uint64_t>(m64, m64 + n64), 1);
    Int result = aInt % mInt;
    // 64-bit → r29 conversion
    std::vector<uint64_t> tmp64(n64, 0);
    size_t rn = result.size();
    if (rn > 0)
        std::memcpy(tmp64.data(), result.words().data(),
                    std::min(rn, n64) * sizeof(uint64_t));
    to_r29(r29_out, tmp64.data(), n64);
}

// Final REDC: r29 Montgomery form → 64-bit plain value
// work: 2*n29+2 + n29 words (T + r29_tmp)
void mont_redc_final_r29(uint64_t* r64_out, size_t n64,
                         const uint64_t* r29_in,
                         const uint64_t* m29, size_t n29,
                         uint64_t m_inv29,
                         uint64_t* work) {
    uint64_t* T = work;
    uint64_t* r29_tmp = work + 2 * n29 + 2;
    std::memcpy(T, r29_in, n29 * sizeof(uint64_t));
    std::memset(T + n29, 0, (n29 + 2) * sizeof(uint64_t));
    redc_r29(r29_tmp, T, m29, n29, m_inv29);
    from_r29(r64_out, r29_tmp, n29, n64);
}

} // namespace r29
} // anonymous namespace
#endif // _MSC_VER && _M_X64

// Window width selection (dynamic, based on the exponent bit length)
// Thresholds equivalent to GMP mpn_sec_powm. The optimal width is determined by
// balancing the precomputation table cost (2^(w-1) multiplications) against
// the reduction in exponent-loop multiplications (expBits/w).
inline int choose_window_width(size_t expBits) {
    if (expBits <= 24)   return 1;
    if (expBits <= 64)   return 3;
    if (expBits <= 256)  return 4;
    if (expBits <= 1024) return 5;
    if (expBits <= 4096) return 6;
    return 7;
}

// Montgomery modular exponentiation: base^exp mod m (m is odd)
// Left-to-right sliding window + Montgomery multiplication
// Combined with dedicated squaring (using mpn::square symmetry)
//
// Advantages of sliding window (vs fixed k-ary):
// - precompute odd powers only (half the table)
// - skip runs of 0 bits (squaring only)
// - restricting the window to odd values reduces the multiplication count by ~10%
//
// n=1 (64-bit modulus) fully inlined powerMod
// Montgomery multiplication completes in two _umul128 calls. Zero heap allocation.
static Int mont_power_mod_1(const Int& base, const Int& exp, const Int& m) {
    uint64_t mod = m.data()[0];
    uint64_t m_inv = mont_neg_inv(mod);

    // Montgomery multiplication: r = (a * b * R^{-1}) mod m, where R = 2^64
    auto mont_mul_1 = [mod, m_inv](uint64_t a, uint64_t b) -> uint64_t {
        // t = a * b (128-bit)
        uint64_t t_hi, t_lo;
#if defined(_MSC_VER)
        t_lo = _umul128(a, b, &t_hi);
#else
        __uint128_t t = static_cast<__uint128_t>(a) * b;
        t_hi = static_cast<uint64_t>(t >> 64);
        t_lo = static_cast<uint64_t>(t);
#endif
        // q = t_lo * m_inv (mod R) — lower 64 bits only
        uint64_t q = t_lo * m_inv;
        // u = (t + q * m) >> 64
        uint64_t qm_hi;
#if defined(_MSC_VER)
        uint64_t qm_lo = _umul128(q, mod, &qm_hi);
#else
        __uint128_t qm = static_cast<__uint128_t>(q) * mod;
        qm_hi = static_cast<uint64_t>(qm >> 64);
        uint64_t qm_lo = static_cast<uint64_t>(qm);
#endif
        uint64_t sum_lo = t_lo + qm_lo;
        uint64_t carry = (sum_lo < t_lo) ? 1ULL : 0ULL;
        uint64_t u = t_hi + qm_hi + carry;
        // conditional subtraction
        return (u >= mod) ? (u - mod) : u;
    };

    // compute R² mod m directly
    auto [q_hi, Rmod] = UInt128::divmod_fast(1ULL, 0ULL, mod);  // R mod m = 2^64 mod m
    uint64_t R2_hi;
#if defined(_MSC_VER)
    uint64_t R2_lo = _umul128(Rmod, Rmod, &R2_hi);
#else
    __uint128_t R2_full = static_cast<__uint128_t>(Rmod) * Rmod;
    R2_hi = static_cast<uint64_t>(R2_full >> 64);
    uint64_t R2_lo = static_cast<uint64_t>(R2_full);
#endif
    auto [q2, R2mod] = UInt128::divmod_fast(R2_hi, R2_lo, mod);

    // base mod m → Montgomery encode
    // negative bases or oversized values need normalization via mod()
    uint64_t bval;
    if (base.getSign() < 0 || base.size() > 1 ||
        (base.size() == 1 && base.data()[0] >= mod)) {
        Int bmod = IntModular::mod(base, m);
        bval = (bmod.size() > 0) ? bmod.data()[0] : 0;
    } else {
        bval = (base.size() == 0) ? 0 : base.data()[0];
    }
    uint64_t baseR = mont_mul_1(bval, R2mod);    // base * R mod m
    uint64_t oneR  = mont_mul_1(1, R2mod);        // 1 * R mod m = R mod m

    // Sliding window precomputation
    size_t expBits = exp.bitLength();
    int w = (expBits <= 8) ? 1 : (expBits <= 24) ? 2 : (expBits <= 64) ? 3 : 4;
    int oddTableSize = 1 << (w - 1);

    uint64_t g[8];  // at most oddTableSize=8 (w=4)
    g[0] = baseR;
    if (oddTableSize > 1) {
        uint64_t base2R = mont_mul_1(baseR, baseR);  // base^2 * R
        for (int i = 1; i < oddTableSize; ++i)
            g[i] = mont_mul_1(g[i - 1], base2R);
    }

    // Left-to-right sliding window exponentiation
    uint64_t result = oneR;
    const uint64_t* edata = exp.data();
    size_t esize = exp.size();

    int i = static_cast<int>(expBits) - 1;
    while (i >= 0) {
        size_t word_idx = static_cast<size_t>(i) / 64;
        unsigned bit_idx = static_cast<unsigned>(i) % 64;
        bool bit = (word_idx < esize) && ((edata[word_idx] >> bit_idx) & 1);

        if (!bit) {
            result = mont_mul_1(result, result);  // sqr
            --i;
        } else {
            // window start: extract an odd value of up to w bits
            int wlen = 1;
            unsigned wval = 1;
            for (int j = 1; j < w && (i - j) >= 0; ++j) {
                size_t wj = static_cast<size_t>(i - j) / 64;
                unsigned bj = static_cast<unsigned>(i - j) % 64;
                unsigned b = (wj < esize) ? ((edata[wj] >> bj) & 1) : 0;
                wval = (wval << 1) | b;
                wlen++;
            }
            // strip trailing 0 bits to make it odd
            while ((wval & 1) == 0) {
                wval >>= 1;
                wlen--;
            }
            // wlen squarings
            for (int j = 0; j < wlen; ++j)
                result = mont_mul_1(result, result);
            // table lookup: g[(wval-1)/2]
            result = mont_mul_1(result, g[(wval - 1) / 2]);
            i -= wlen;
        }
    }

    // Final REDC: result * 1 → result * R^{-1} mod m
    uint64_t final_val = mont_mul_1(result, 1);
    if (final_val == 0) return Int::Zero();
    return Int(final_val);
}

// AVX2 radix-2^29 path:
// - unified with R_r29 = 2^(29*n29) (avoids mismatch with R_64)
// - encode/precompute/loop/final REDC all run in the r29 domain
// - throughput improved via VPMULUDQ (4×32→64 parallel multiplication)
// n <= 8 dedicated fast path: zero heap allocation, direct ASM kernel calls
// Completes all processing in fixed-size buffers on the stack
#ifdef SANGI_INT_HAS_ASM
template<size_t N>  // N = number of words (2, 4, 8)
static Int mont_power_mod_fixed(const Int& base, const Int& exp, const Int& m) {
    static_assert(N == 2 || N == 4 || N == 8);

    auto _prof_lap = [](const char*) {};

    const uint64_t* mdata = m.data();
    uint64_t m_inv = mont_neg_inv(mdata[0]);

    size_t expBits = exp.bitLength();
    int w = choose_window_width(expBits);
    int oddTableSize = 1 << (w - 1);

    // place all buffers on the stack (zero heap allocation)
    alignas(64) uint64_t R2[N] = {};
    alignas(64) uint64_t baseR[N] = {};
    alignas(64) uint64_t oneR[N] = {};
    alignas(64) uint64_t result[N], temp_buf[N];
    alignas(64) uint64_t pad[N] = {};
    // g_buf: at most oddTableSize (w=7 → 64) entries.
    // w depends on expBits, but for N<=8, w<=5 (oddTableSize<=16)
    alignas(64) uint64_t g_buf[16 * N] = {};  // at most 16 entries
    alignas(64) uint64_t base2R[N];

    // direct calls to Montgomery multiplication/squaring (no dispatch needed)
    auto do_mul = [&](uint64_t* dst, const uint64_t* a, const uint64_t* b) {
        if constexpr (N == 2) mpn_mont_mul_2(dst, a, b, mdata, m_inv);
        else if constexpr (N == 4) mpn_mont_mul_4(dst, a, b, mdata, m_inv);
        else mpn_mont_mul_8(dst, a, b, mdata, m_inv);
    };
    auto do_sqr = [&](uint64_t* dst, const uint64_t* src) {
        if constexpr (N == 2) mpn_mont_sqr_2(dst, src, mdata, m_inv);
        else if constexpr (N == 4) mpn_mont_sqr_4(dst, src, mdata, m_inv);
        else mpn_mont_sqr_8(dst, src, mdata, m_inv);
    };

    // R² mod m: R = 2^(64*N), R² = R*R mod m
    // Method: R_buf = 2^(64*N) (N+1 words) → R mod m → (R mod m)² mod m
    // Faster than repeated doubling (a 2*64*N loop): done in 2 divisions + 1 multiplication
    _prof_lap("setup");
    {
        alignas(64) uint64_t R_buf[N + 1] = {};
        R_buf[N] = 1;  // R = 2^(64*N)
        alignas(64) uint64_t q_buf[2] = {};
        alignas(64) uint64_t Rmod[N] = {};

        // R mod m (N+1 ÷ N word division)
        size_t div_scratch_sz = mpn::divide_scratch_size(N + 1, N);
        alignas(64) uint64_t div_scratch[512];  // sufficient for N<=8
        mpn::divide(q_buf, Rmod, R_buf, N + 1, mdata, N, div_scratch);

        // (R mod m)² mod m
        alignas(64) uint64_t Rmod_sq[2 * N] = {};
        size_t mul_scratch_sz = mpn::multiply_scratch_size(N, N);
        alignas(64) uint64_t mul_scratch[512];
        size_t rmod_n = N;
        while (rmod_n > 0 && Rmod[rmod_n - 1] == 0) --rmod_n;
        if (rmod_n == 0) rmod_n = 1;
        mpn::multiply(Rmod_sq, Rmod, rmod_n, Rmod, rmod_n, mul_scratch);
        size_t sq_n = 2 * rmod_n;
        while (sq_n > 0 && Rmod_sq[sq_n - 1] == 0) --sq_n;
        if (sq_n == 0) sq_n = 1;

        if (sq_n <= N) {
            std::memcpy(R2, Rmod_sq, sq_n * sizeof(uint64_t));
        } else {
            alignas(64) uint64_t q2[N + 1] = {};
            mpn::divide(q2, R2, Rmod_sq, sq_n, mdata, N, div_scratch);
        }
    }

    _prof_lap("R2_compute");
    // convert base to Montgomery form
    {
        Int bmod = IntModular::mod(base, m);
        _prof_lap("base_mod");
        size_t bn = bmod.size();
        std::memset(pad, 0, N * sizeof(uint64_t));
        if (bn > 0)
            std::memcpy(pad, bmod.data(), std::min(bn, N) * sizeof(uint64_t));
        do_mul(baseR, pad, R2);
    }

    // 1R
    {
        std::memset(pad, 0, N * sizeof(uint64_t));
        pad[0] = 1;
        do_mul(oneR, pad, R2);
    }
    _prof_lap("encode");

    // window table
    uint64_t* g[16];
    for (int ii = 0; ii < oddTableSize && ii < 16; ++ii)
        g[ii] = g_buf + static_cast<size_t>(ii) * N;

    std::memcpy(g[0], baseR, N * sizeof(uint64_t));
    if (oddTableSize > 1) {
        do_sqr(base2R, baseR);
        for (int ii = 1; ii < oddTableSize; ++ii)
            do_mul(g[ii], g[ii - 1], base2R);
    }
    _prof_lap("precomp");

    // Sliding window exponentiation
    uint64_t* rp = result;
    uint64_t* tp = temp_buf;
    std::memcpy(rp, oneR, N * sizeof(uint64_t));

    int i = static_cast<int>(expBits) - 1;
    while (i >= 0) {
        if (!exp.getBit(i)) {
            do_sqr(tp, rp);
            std::swap(rp, tp);
            --i;
        } else {
            int j = (i - w + 1 > 0) ? (i - w + 1) : 0;
            while (!exp.getBit(j)) ++j;
            int wval = 0;
            for (int k = i; k >= j; --k)
                wval = (wval << 1) | (exp.getBit(k) ? 1 : 0);
            for (int s = 0; s < i - j + 1; ++s) {
                do_sqr(tp, rp);
                std::swap(rp, tp);
            }
            do_mul(tp, rp, g[(wval - 1) / 2]);
            std::swap(rp, tp);
            i = j - 1;
        }
    }

    _prof_lap("loop");
    // Montgomery → normal form (REDC with input = result * 1)
    {
        std::memset(pad, 0, N * sizeof(uint64_t));
        pad[0] = 1;
        do_mul(tp, rp, pad);
    }

    size_t rn = N;
    while (rn > 0 && tp[rn - 1] == 0) --rn;
    if (rn == 0) return Int::Zero();
    std::vector<uint64_t> rwords(tp, tp + rn);
    _prof_lap("final_redc");
    return Int::fromRawWords(rwords, 1);
}
#endif // SANGI_INT_HAS_ASM

Int mont_power_mod(const Int& base, const Int& exp, const Int& m) {
    size_t n = m.size();      // number of words in m
    const uint64_t* mdata = m.data();

    // n=1 fully inlined (no ASM needed, only _umul128)
    if (n == 1) return mont_power_mod_1(base, exp, m);

    // n <= 8 dedicated fast path (zero heap allocation, direct ASM kernel calls)
#ifdef SANGI_INT_HAS_ASM
    if (mpn::detail::has_bmi2_adx()) {
        if (n == 2) return mont_power_mod_fixed<2>(base, exp, m);
        if (n == 4) return mont_power_mod_fixed<4>(base, exp, m);
        if (n == 8) return mont_power_mod_fixed<8>(base, exp, m);
    }
#endif

    uint64_t m_inv = mont_neg_inv(mdata[0]);
    uint64_t mip1 = (n >= 2) ? mont_neg_inv2(mdata[0], mdata[1], m_inv) : 0;

    size_t expBits = exp.bitLength();
    int w = choose_window_width(expBits);
    int oddTableSize = 1 << (w - 1);  // odd powers only: 2^(w-1) entries

    // --- AVX2 radix-2^29 domain decision ---
#if defined(_MSC_VER) && defined(_M_X64)
    bool use_avx2 = false;  // schoolbook r29 is slower than scalar MULX+Karatsuba (1024-bit: 2x, 2048-bit: 1.9x)
    size_t n29 = 0;
    uint64_t m_inv29 = 0;
    std::vector<uint64_t> m29_buf;
    // work: 2*n29+2 (T) + n29 (r29_tmp for final REDC) = 3*n29+2
    std::vector<uint64_t> avx2_work;
    if (use_avx2) {
        n29 = (n * 64 + r29::BITS - 1) / r29::BITS;
        m29_buf.resize(n29 + 1, 0);
        r29::to_r29(m29_buf.data(), mdata, n);
        m_inv29 = r29::neg_inv_r29(m29_buf[0]);
        avx2_work.resize(3 * n29 + 4, 0);
    }
#else
    constexpr bool use_avx2 = false;
    constexpr size_t n29 = 0;
#endif

    // limb count (AVX2: r29 limbs, scalar: 64-bit limbs)
    size_t limb_n = use_avx2 ? n29 : n;

    // ========================================================================
    // Single buffer allocation: consolidates all scratch/work areas into one heap allocation
    // Previously 15-18 std::vector allocations occurred; reduced to one
    // ========================================================================
    size_t sq_scratch_sz = 0, mul_scratch_sz = 0, extra_sz = 0;
    if (!use_avx2) {
        sq_scratch_sz = mpn::square_scratch_size(n);
        mul_scratch_sz = mpn::multiply_scratch_size(n, n);
        extra_sz = std::max(sq_scratch_sz, mul_scratch_sz);
    }

    // scratch size for R² computation (when n >= 2)
    size_t r2_scratch_sz = 0;
    if (!use_avx2 && n >= 2) {
        size_t ds1 = mpn::divide_scratch_size(n + 1, n);
        size_t ds2 = mpn::divide_scratch_size(2 * n, n);
        size_t msz = mpn::multiply_scratch_size(n, n);
        r2_scratch_sz = std::max({ds1, ds2, msz});
    }

    // for redc_n: precompute the Hensel inverse when n >= REDC_N_THRESHOLD
    bool use_redc_n = (!use_avx2 && n >= REDC_N_THRESHOLD);
    size_t inv_n_sz = use_redc_n ? n : 0;         // m_inv_n buffer
    size_t inv_scratch_sz = use_redc_n ? 3 * n : 0; // hensel_inverse scratch

    // buffer layout (scalar path):
    //   [0..2n]                scratch (product of mont_mul/sqr + REDC)
    //   [2n+1..2n+extra]       sq_scratch / mul_scratch
    //   [+n]                   R2
    //   [+n]                   baseR
    //   [+n]                   oneR
    //   [+oddTableSize*n]      g_buf (window table)
    //   [+n]                   base2R
    //   [+n]                   result
    //   [+n]                   temp
    //   [+n]                   base_padded (for encode)
    //   [+inv_n_sz]            m_inv_n (Hensel inverse for redc_n)
    //   [+n+1+2+n+2n+r2_scr]  temporary area for R² computation
    //   [+inv_scratch_sz]      hensel_inverse scratch (temporary, not needed after R² computation)
    size_t mont_area = 2 * n + 1 + extra_sz;
    size_t buf_total = mont_area
                     + n                                    // R2
                     + n                                    // baseR
                     + n                                    // oneR
                     + static_cast<size_t>(oddTableSize) * n // g_buf
                     + n                                    // base2R
                     + n                                    // result
                     + n                                    // temp
                     + n                                    // base_padded / one_padded
                     + inv_n_sz                             // m_inv_n
                     + (n >= 2 ? (n + 1) + 2 + n + 2 * n + r2_scratch_sz : 0)
                     + inv_scratch_sz;                      // hensel scratch (temporary)

    std::vector<uint64_t> arena(buf_total, 0);
    uint64_t* ap = arena.data();

    // assign pointers to each area
    uint64_t* scratch_ptr = ap;          ap += mont_area;
    uint64_t* sq_scratch_p  = scratch_ptr + 2 * n + 1;
    uint64_t* mul_scratch_p = scratch_ptr + 2 * n + 1;
    uint64_t* R2_p      = ap;           ap += n;
    uint64_t* baseR_p   = ap;           ap += n;
    uint64_t* oneR_p    = ap;           ap += n;
    uint64_t* g_buf_p   = ap;           ap += static_cast<size_t>(oddTableSize) * n;
    uint64_t* base2R_p  = ap;           ap += n;
    uint64_t* result_p  = ap;           ap += n;
    uint64_t* temp_p    = ap;           ap += n;
    uint64_t* pad_p     = ap;           ap += n;
    // Hensel inverse for redc_n
    uint64_t* m_inv_n_p = use_redc_n ? ap : nullptr;
    if (use_redc_n) ap += inv_n_sz;
    // temporary area for R² computation (n >= 2)
    uint64_t* r2_R_buf  = ap;           // n+1 words
    uint64_t* r2_q      = r2_R_buf + n + 1;  // 2 words
    uint64_t* r2_Rmod   = r2_q + 2;          // n words
    uint64_t* r2_full   = r2_Rmod + n;       // 2n words
    uint64_t* r2_work   = r2_full + 2 * n;   // r2_scratch_sz words

    // precompute R² mod m
    if (!use_avx2) {
        if (n == 1) {
            auto [q_hi, Rmod1] = UInt128::divmod_fast(1ULL, 0ULL, mdata[0]);
            uint64_t hi, lo;
#if defined(_MSC_VER)
            lo = _umul128(Rmod1, Rmod1, &hi);
#else
            __uint128_t prod = static_cast<__uint128_t>(Rmod1) * Rmod1;
            hi = static_cast<uint64_t>(prod >> 64);
            lo = static_cast<uint64_t>(prod);
#endif
            auto [q2, R2mod1] = UInt128::divmod_fast(hi, lo, mdata[0]);
            R2_p[0] = R2mod1;
        } else {
            // R = 2^(64n)
            r2_R_buf[n] = 1;  // the rest is already zero-initialized by arena

            // R mod m
            mpn::divide(r2_q, r2_Rmod, r2_R_buf, n + 1, mdata, n, r2_work);

            // (R mod m)² mod m
            size_t rmod_n = n;
            while (rmod_n > 0 && r2_Rmod[rmod_n - 1] == 0) --rmod_n;
            if (rmod_n == 0) rmod_n = 1;

            mpn::multiply(r2_full, r2_Rmod, rmod_n, r2_Rmod, rmod_n, r2_work);
            size_t r2fn = 2 * rmod_n;
            while (r2fn > 0 && r2_full[r2fn - 1] == 0) --r2fn;
            if (r2fn == 0) r2fn = 1;

            if (r2fn <= n) {
                std::memcpy(R2_p, r2_full, r2fn * sizeof(uint64_t));
            } else {
                mpn::divide(r2_q, R2_p, r2_full, r2fn, mdata, n, r2_work);
            }
        }
    }

    // precompute the Hensel inverse for redc_n
    if (use_redc_n) {
        // use the inv_scratch area at the end of arena as a temporary buffer
        uint64_t* inv_scratch = arena.data() + buf_total - inv_scratch_sz;
        hensel_inverse(m_inv_n_p, mdata, n, m_inv, inv_scratch);
    }

    // convert base to Montgomery form
    {
        Int bmod = IntModular::mod(base, m);
        size_t bn = bmod.size();
        const uint64_t* bdata = bmod.data();
#if defined(_MSC_VER) && defined(_M_X64)
        if (use_avx2) {
            r29::mont_encode_r29(baseR_p, n29, bdata, bn, mdata, n);
        } else
#endif
        {
            std::memset(pad_p, 0, n * sizeof(uint64_t));
            if (bn > 0)
                std::memcpy(pad_p, bdata, std::min(bn, n) * sizeof(uint64_t));
            mont_mul(baseR_p, pad_p, R2_p,
                     mdata, n, m_inv, mip1, scratch_ptr, mul_scratch_p);
        }
    }

    // 1R = Mont(1, R²) = R mod m
    {
#if defined(_MSC_VER) && defined(_M_X64)
        if (use_avx2) {
            uint64_t one = 1;
            r29::mont_encode_r29(oneR_p, n29, &one, 1, mdata, n);
        } else
#endif
        {
            std::memset(pad_p, 0, n * sizeof(uint64_t));
            pad_p[0] = 1;
            mont_mul(oneR_p, pad_p, R2_p,
                     mdata, n, m_inv, mip1, scratch_ptr, mul_scratch_p);
        }
    }

    // === Sliding window: precompute odd powers only ===
    uint64_t** g = static_cast<uint64_t**>(alloca(oddTableSize * sizeof(uint64_t*)));
    for (int ii = 0; ii < oddTableSize; ++ii)
        g[ii] = g_buf_p + static_cast<size_t>(ii) * limb_n;

    // lambda: dispatch for mont_sqr / mont_mul
    auto do_sqr = [&](uint64_t* dst, const uint64_t* src) {
#if defined(_MSC_VER) && defined(_M_X64)
        if (use_avx2) {
            r29::mont_sqr_r29(dst, src, m29_buf.data(), n29,
                              m_inv29, avx2_work.data());
            return;
        }
#endif
        mont_sqr(dst, src, mdata, n, m_inv, mip1, scratch_ptr, sq_scratch_p,
                 m_inv_n_p, mul_scratch_p);
    };

    auto do_mul = [&](uint64_t* dst, const uint64_t* a, const uint64_t* b) {
#if defined(_MSC_VER) && defined(_M_X64)
        if (use_avx2) {
            r29::mont_mul_r29(dst, a, b, m29_buf.data(), n29,
                              m_inv29, avx2_work.data());
            return;
        }
#endif
        mont_mul(dst, a, b, mdata, n, m_inv, mip1, scratch_ptr, mul_scratch_p,
                 m_inv_n_p);
    };

    // g[0] = baseR
    std::memcpy(g[0], baseR_p, limb_n * sizeof(uint64_t));

    if (oddTableSize > 1) {
        do_sqr(base2R_p, baseR_p);
        for (int ii = 1; ii < oddTableSize; ++ii)
            do_mul(g[ii], g[ii - 1], base2R_p);
    }

    // === Left-to-right sliding window exponentiation ===
    std::memcpy(result_p, oneR_p, limb_n * sizeof(uint64_t));

    int i = static_cast<int>(expBits) - 1;
    while (i >= 0) {
        if (!exp.getBit(i)) {
            do_sqr(temp_p, result_p);
            std::swap(result_p, temp_p);
            --i;
        } else {
            int j = (i - w + 1 > 0) ? (i - w + 1) : 0;
            while (!exp.getBit(j)) ++j;

            int wval = 0;
            for (int k = i; k >= j; --k)
                wval = (wval << 1) | (exp.getBit(k) ? 1 : 0);

            int sqr_count = i - j + 1;
            for (int s = 0; s < sqr_count; ++s) {
                do_sqr(temp_p, result_p);
                std::swap(result_p, temp_p);
            }

            do_mul(temp_p, result_p, g[(wval - 1) / 2]);
            std::swap(result_p, temp_p);

            i = j - 1;
        }
    }

    // === convert from Montgomery form back to normal form ===
#if defined(_MSC_VER) && defined(_M_X64)
    if (use_avx2) {
        std::vector<uint64_t> r64(n, 0);
        r29::mont_redc_final_r29(r64.data(), n, result_p,
                                 m29_buf.data(), n29, m_inv29,
                                 avx2_work.data());
        while (n > 0 && r64[n - 1] == 0) --n;
        if (n == 0) return Int::Zero();
        return Int::fromRawWords(r64, 1);
    }
#endif
    {
        std::memset(scratch_ptr, 0, (2 * n + 1) * sizeof(uint64_t));
        std::memcpy(scratch_ptr, result_p, n * sizeof(uint64_t));
        mont_redc_2(temp_p, scratch_ptr, mdata, n, m_inv, mip1);

        size_t rn = n;
        while (rn > 0 && temp_p[rn - 1] == 0) --rn;
        if (rn == 0) return Int::Zero();
        std::vector<uint64_t> rwords(temp_p, temp_p + rn);
        return Int::fromRawWords(rwords, 1);
    }
}

// Constant-time modular exponentiation: base^exp mod m (equivalent to GMP mpz_powm_sec)
// - fixed window width (not a sliding window)
// - reads every table entry each time, with conditional selection (cmov-equivalent)
// - executes sqr+mul regardless of bit 0/1, selecting the result via conditional copy
Int mont_power_mod_sec(const Int& base, const Int& exp, const Int& m) {
    size_t n = m.size();
    const uint64_t* mdata = m.data();
    uint64_t m_inv = mont_neg_inv(mdata[0]);
    uint64_t mip1 = (n >= 2) ? mont_neg_inv2(mdata[0], mdata[1], m_inv) : 0;

    size_t expBits = exp.bitLength();
    int w = choose_window_width(expBits);
    int tableSize = 1 << w;  // all entries (0..2^w-1)

    // scratch buffer
    size_t sq_scratch_sz = mpn::square_scratch_size(n);
    size_t mul_scratch_sz = mpn::multiply_scratch_size(n, n);
    size_t extra_sz = std::max(sq_scratch_sz, mul_scratch_sz);
    std::vector<uint64_t> scratch(2 * n + 1 + extra_sz, 0);
    uint64_t* sq_scratch = scratch.data() + 2 * n + 1;
    uint64_t* mul_scratch_ptr = scratch.data() + 2 * n + 1;

    // precompute R² mod m
    std::vector<uint64_t> R2(n, 0);
    {
        Int mInt = Int::fromRawWords(std::vector<uint64_t>(mdata, mdata + n), 1);
        Int R_int = Int::One();
        IntOps::leftShift(R_int, static_cast<int>(64 * n));
        Int R_mod = R_int % mInt;
        Int R2_mod = (R_mod * R_mod) % mInt;
        size_t r2n = R2_mod.size();
        if (r2n > 0)
            std::memcpy(R2.data(), R2_mod.words().data(),
                        std::min(r2n, n) * sizeof(uint64_t));
    }

    // base → Montgomery form
    std::vector<uint64_t> baseR(n, 0);
    {
        Int bmod = IntModular::mod(base, m);
        std::vector<uint64_t> base_padded(n, 0);
        size_t bn = bmod.size();
        if (bn > 0) {
            auto bw = bmod.words();
            std::memcpy(base_padded.data(), bw.data(),
                        std::min(bn, n) * sizeof(uint64_t));
        }
        mont_mul(baseR.data(), base_padded.data(), R2.data(),
                 mdata, n, m_inv, mip1, scratch.data(), mul_scratch_ptr);
    }

    // 1R = 1 in Montgomery form
    std::vector<uint64_t> oneR(n, 0);
    {
        std::vector<uint64_t> one_padded(n, 0);
        one_padded[0] = 1;
        mont_mul(oneR.data(), one_padded.data(), R2.data(),
                 mdata, n, m_inv, mip1, scratch.data(), mul_scratch_ptr);
    }

    // === fixed window: precompute all powers ===
    // g[i] = base^i * R mod m  (i = 0..2^w-1)
    std::vector<uint64_t> g_buf(static_cast<size_t>(tableSize) * n);
    std::vector<uint64_t*> g(tableSize);
    for (int i = 0; i < tableSize; ++i)
        g[i] = g_buf.data() + static_cast<size_t>(i) * n;

    // g[0] = 1R
    std::memcpy(g[0], oneR.data(), n * sizeof(uint64_t));
    // g[1] = baseR
    std::memcpy(g[1], baseR.data(), n * sizeof(uint64_t));
    // g[i] = g[i-1] * baseR
    for (int i = 2; i < tableSize; ++i)
        mont_mul(g[i], g[i - 1], baseR.data(), mdata, n, m_inv, mip1, scratch.data(), mul_scratch_ptr);

    // === constant-time table lookup ===
    // read every entry and select only the one with idx == target via OR
    auto ct_select = [&](uint64_t* dst, int idx) {
        std::memset(dst, 0, n * sizeof(uint64_t));
        for (int i = 0; i < tableSize; ++i) {
            uint64_t mask = static_cast<uint64_t>(-(static_cast<int64_t>(i == idx)));
            for (size_t j = 0; j < n; ++j)
                dst[j] |= g[i][j] & mask;
        }
    };

    // === constant-time conditional copy (kept for future ct final subtraction) ===
    [[maybe_unused]] auto ct_cond_copy = [&](uint64_t* dst, const uint64_t* src, bool cond) {
        uint64_t mask = static_cast<uint64_t>(-(static_cast<int64_t>(cond)));
        uint64_t nmask = ~mask;
        for (size_t j = 0; j < n; ++j)
            dst[j] = (dst[j] & nmask) | (src[j] & mask);
    };

    // === fixed-window modular exponentiation (left-to-right) ===
    // process the exponent w bits at a time. pad to a multiple of w.
    size_t padded_bits = ((expBits + w - 1) / w) * w;

    std::vector<uint64_t> result(n);
    std::vector<uint64_t> temp(n);
    std::vector<uint64_t> sel(n);
    std::memcpy(result.data(), oneR.data(), n * sizeof(uint64_t));

    for (size_t pos = padded_bits; pos >= static_cast<size_t>(w); pos -= w) {
        // w squarings
        for (int s = 0; s < w; ++s) {
            mont_sqr(temp.data(), result.data(), mdata, n, m_inv, mip1,
                     scratch.data(), sq_scratch);
            std::memcpy(result.data(), temp.data(), n * sizeof(uint64_t));
        }

        // extract the window value
        int wval = 0;
        for (int b = 0; b < w; ++b) {
            int bit_idx = static_cast<int>(pos) - w + b;
            int bit = (bit_idx >= 0 && bit_idx < static_cast<int>(expBits))
                      ? (exp.getBit(bit_idx) ? 1 : 0) : 0;
            wval |= bit << b;
        }

        // constant-time table selection
        ct_select(sel.data(), wval);

        // always perform the multiplication
        mont_mul(temp.data(), result.data(), sel.data(),
                 mdata, n, m_inv, mip1, scratch.data(), mul_scratch_ptr);
        std::memcpy(result.data(), temp.data(), n * sizeof(uint64_t));
    }

    // === Montgomery → normal form ===
    std::memset(scratch.data(), 0, (2 * n + 1) * sizeof(uint64_t));
    std::memcpy(scratch.data(), result.data(), n * sizeof(uint64_t));
    mont_redc_2(temp.data(), scratch.data(), mdata, n, m_inv, mip1);

    std::vector<uint64_t> rwords(temp.data(), temp.data() + n);
    while (!rwords.empty() && rwords.back() == 0) rwords.pop_back();
    if (rwords.empty()) return Int::Zero();
    return Int::fromRawWords(rwords, 1);
}

} // anonymous namespace

// Normalized remainder (same behavior as Mathematica's Mod[])
Int IntModular::mod(const Int& x, const Int& m) {
    // 1. handle special states
    if (x.isNaN() || m.isNaN()) {
        return Int::NaN();
    }
    if (x.isInfinite() || m.isInfinite()) {
        return Int::NaN();
    }
    if (m.isZero()) {
        // division by zero
        return Int::NaN();
    }

    // 2. handle the case where m is negative
    // mod(x, -m) = -mod(-x, m)
    if (m.getSign() < 0) {
        return -mod(-x, -m);
    }

    // 3. handle the case where x is negative
    // mod(-a, m) = m - mod(a, m) (but 0 if mod(a,m)=0)
    if (x.getSign() < 0) {
        Int t = mod(-x, m);
        if (t.isZero()) {
            return t;
        } else {
            return m - t;
        }
    }

    // 4. ordinary remainder of two positive values
    // use x % m (C++ operator%)
    return x % m;
}

// Modular exponentiation: base^exp mod m (Binary exponentiation)
Int IntModular::powerMod(const Int& base, const Int& exp, const Int& m) {
    // 1. handle special states
    if (base.isNaN() || exp.isNaN() || m.isNaN()) {
        return Int::NaN();
    }
    if (base.isInfinite() || exp.isInfinite() || m.isInfinite()) {
        return Int::NaN();
    }
    if (m.isZero()) {
        // division by zero
        return Int::NaN();
    }

    // 2. special case m = 1 (any number mod 1 = 0)
    if (m.isOne()) {
        return Int::Zero();
    }

    // 3. case of exponent 0: base^0 = 1
    if (exp.isZero()) {
        return Int::One();
    }

    // 4. case of base 0: 0^n = 0 (n > 0)
    if (base.isZero()) {
        return Int::Zero();
    }

    // 5. handle negative exponent: base^(-n) = (base^n)^(-1) mod m
    if (exp.getSign() < 0) {
        // first compute base^|exp| mod m
        Int pos_result = powerMod(base, -exp, m);
        // return its inverse
        return inverseMod(pos_result, m);
    }

    // 6. if m is odd, use Montgomery modular exponentiation (fast, no division)
    if (m.getBit(0)) {
        return mont_power_mod(base, exp, m);
    }

    // 7. m even: conventional Binary exponentiation (right-to-left method)
    //    special states already checked at entry → use the Unchecked version
    Int x = mod(base, m);
    Int result = Int::One();
    Int e = exp;
    Int temp;

    while (!e.isZero()) {
        if (e.getBit(0)) {
            IntOps::mulUnchecked(result, x, temp);
            result = mod(temp, m);
        }
        IntOps::squareUnchecked(x, temp);
        x = mod(temp, m);
        IntOps::rightShiftUnchecked(e, 1, e);
    }

    return result;
}

// Multiplicative inverse: a^(-1) mod m
Int IntModular::inverseMod(const Int& a, const Int& m, bool is_prime) {
    // 1. handle special states
    if (a.isNaN() || m.isNaN()) {
        return Int::NaN();
    }
    if (a.isInfinite() || m.isInfinite()) {
        return Int::NaN();
    }
    if (m.isZero()) {
        // division by zero
        return Int::NaN();
    }

    // 2. special case m = 1 (any number mod 1 = 0, no inverse)
    if (m.isOne()) {
        return Int::Zero();
    }

    // 3. if a = 0, no inverse exists
    if (a.isZero()) {
        return Int::Zero();
    }

    // 4. normalize a (into the range 0 <= a < m)
    Int a_norm = mod(a, m);
    if (a_norm.isZero()) {
        return Int::Zero();  // a ≡ 0 mod m, so no inverse
    }

    // 5. if a = 1, the inverse is 1
    if (a_norm.isOne()) {
        return Int::One();
    }

    // 6. optimization when m is prime: a^(-1) ≡ a^(m-2) mod m (Fermat's little theorem)
    if (is_prime) {
        return powerMod(a_norm, m - 2, m);
    }

    // 7. general case: extended Euclidean algorithm
    // solve a * x + m * y = gcd(a, m)
    Int x, y;
    Int g = IntGCD::extendedGcd(a_norm, m, x, y);

    // 8. if gcd(a, m) != 1, no inverse exists
    if (!g.isOne()) {
        return Int::Zero();  // no inverse
    }

    // 9. normalize x and return it (0 <= x < m)
    Int result = mod(x, m);
    return result;
}

// Constant-time modular exponentiation: base^exp mod m (equivalent to GMP mpz_powm_sec)
Int IntModular::powerModSec(const Int& base, const Int& exp, const Int& m) {
    // 1. special states
    if (base.isNaN() || exp.isNaN() || m.isNaN()) return Int::NaN();
    if (base.isInfinite() || exp.isInfinite() || m.isInfinite()) return Int::NaN();
    if (m.isZero()) return Int::NaN();

    // 2. m = 1 → 0
    if (m.isOne()) return Int::Zero();

    // 3. exp = 0 → 1
    if (exp.isZero()) return Int::One();

    // 4. negative exponent not allowed (constant-time constraint)
    if (exp.getSign() < 0) return Int::NaN();

    // 5. m must be an odd positive integer
    if (!m.getBit(0)) {
        // even m → fall back to powerMod
        return powerMod(base, exp, m);
    }

    // 6. Montgomery constant-time modular exponentiation
    return mont_power_mod_sec(base, exp, m);
}

} // namespace sangi