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

// IntCombinatorics.cpp
// Implementation of combinatorial functions
//
// Factorial algorithms:
//   References:
//     http://www.luschny.de/math/factorial/FastFactorialFunctions.htm
//     http://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/1138-22.pdf (tournament method)
//     http://www.apfloat.org/ (binary splitting method)

#include "math/core/mp/Int/IntCombinatorics.hpp"
#include "math/core/mp/Int/IntOps.hpp"
#include "math/core/mp/Int/IntSequence.hpp"
#include <algorithm>
#include <vector>
#include <cmath>
#include <bit>
#include <cassert>

namespace sangi {

// ============================================================================
// Factorial table (0! to 100!)
// ============================================================================

static const Int& getFactorialTable(int n) {
    assert(n >= 0 && n <= 100);
    static const Int table[101] = {
        Int("1"),
        Int("1"),
        Int("2"),
        Int("6"),
        Int("24"),
        Int("120"),
        Int("720"),
        Int("5040"),
        Int("40320"),
        Int("362880"),
        Int("3628800"),
        Int("39916800"),
        Int("479001600"),
        Int("6227020800"),
        Int("87178291200"),
        Int("1307674368000"),
        Int("20922789888000"),
        Int("355687428096000"),
        Int("6402373705728000"),
        Int("121645100408832000"),
        Int("2432902008176640000"),
        Int("51090942171709440000"),
        Int("1124000727777607680000"),
        Int("25852016738884976640000"),
        Int("620448401733239439360000"),
        Int("15511210043330985984000000"),
        Int("403291461126605635584000000"),
        Int("10888869450418352160768000000"),
        Int("304888344611713860501504000000"),
        Int("8841761993739701954543616000000"),
        Int("265252859812191058636308480000000"),
        Int("8222838654177922817725562880000000"),
        Int("263130836933693530167218012160000000"),
        Int("8683317618811886495518194401280000000"),
        Int("295232799039604140847618609643520000000"),
        Int("10333147966386144929666651337523200000000"),
        Int("371993326789901217467999448150835200000000"),
        Int("13763753091226345046315979581580902400000000"),
        Int("523022617466601111760007224100074291200000000"),
        Int("20397882081197443358640281739902897356800000000"),
        Int("815915283247897734345611269596115894272000000000"),
        Int("33452526613163807108170062053440751665152000000000"),
        Int("1405006117752879898543142606244511569936384000000000"),
        Int("60415263063373835637355132068513997507264512000000000"),
        Int("2658271574788448768043625811014615890319638528000000000"),
        Int("119622220865480194561963161495657715064383733760000000000"),
        Int("5502622159812088949850305428800254892961651752960000000000"),
        Int("258623241511168180642964355153611979969197632389120000000000"),
        Int("12413915592536072670862289047373375038521486354677760000000000"),
        Int("608281864034267560872252163321295376887552831379210240000000000"),
        Int("30414093201713378043612608166064768844377641568960512000000000000"),
        Int("1551118753287382280224243016469303211063259720016986112000000000000"),
        Int("80658175170943878571660636856403766975289505440883277824000000000000"),
        Int("4274883284060025564298013753389399649690343788366813724672000000000000"),
        Int("230843697339241380472092742683027581083278564571807941132288000000000000"),
        Int("12696403353658275925965100847566516959580321051449436762275840000000000000"),
        Int("710998587804863451854045647463724949736497978881168458687447040000000000000"),
        Int("40526919504877216755680601905432322134980384796226602145184481280000000000000"),
        Int("2350561331282878571829474910515074683828862318181142924420699914240000000000000"),
        Int("138683118545689835737939019720389406345902876772687432540821294940160000000000000"),
        Int("8320987112741390144276341183223364380754172606361245952449277696409600000000000000"),
        Int("507580213877224798800856812176625227226004528988036003099405939480985600000000000000"),
        Int("31469973260387937525653122354950764088012280797258232192163168247821107200000000000000"),
        Int("1982608315404440064116146708361898137544773690227268628106279599612729753600000000000000"),
        Int("126886932185884164103433389335161480802865516174545192198801894375214704230400000000000000"),
        Int("8247650592082470666723170306785496252186258551345437492922123134388955774976000000000000000"),
        Int("544344939077443064003729240247842752644293064388798874532860126869671081148416000000000000000"),
        Int("36471110918188685288249859096605464427167635314049524593701628500267962436943872000000000000000"),
        Int("2480035542436830599600990418569171581047399201355367672371710738018221445712183296000000000000000"),
        Int("171122452428141311372468338881272839092270544893520369393648040923257279754140647424000000000000000"),
        Int("11978571669969891796072783721689098736458938142546425857555362864628009582789845319680000000000000000"),
        Int("850478588567862317521167644239926010288584608120796235886430763388588680378079017697280000000000000000"),
        Int("61234458376886086861524070385274672740778091784697328983823014963978384987221689274204160000000000000000"),
        Int("4470115461512684340891257138125051110076800700282905015819080092370422104067183317016903680000000000000000"),
        Int("330788544151938641225953028221253782145683251820934971170611926835411235700971565459250872320000000000000000"),
        Int("24809140811395398091946477116594033660926243886570122837795894512655842677572867409443815424000000000000000000"),
        Int("1885494701666050254987932260861146558230394535379329335672487982961844043495537923117729972224000000000000000000"),
        Int("145183092028285869634070784086308284983740379224208358846781574688061991349156420080065207861248000000000000000000"),
        Int("11324281178206297831457521158732046228731749579488251990048962825668835325234200766245086213177344000000000000000000"),
        Int("894618213078297528685144171539831652069808216779571907213868063227837990693501860533361810841010176000000000000000000"),
        Int("71569457046263802294811533723186532165584657342365752577109445058227039255480148842668944867280814080000000000000000000"),
        Int("5797126020747367985879734231578109105412357244731625958745865049716390179693892056256184534249745940480000000000000000000"),
        Int("475364333701284174842138206989404946643813294067993328617160934076743994734899148613007131808479167119360000000000000000000"),
        Int("39455239697206586511897471180120610571436503407643446275224357528369751562996629334879591940103770870906880000000000000000000"),
        Int("3314240134565353266999387579130131288000666286242049487118846032383059131291716864129885722968716753156177920000000000000000000"),
        Int("281710411438055027694947944226061159480056634330574206405101912752560026159795933451040286452340924018275123200000000000000000000"),
        Int("24227095383672732381765523203441259715284870552429381750838764496720162249742450276789464634901319465571660595200000000000000000000"),
        Int("2107757298379527717213600518699389595229783738061356212322972511214654115727593174080683423236414793504734471782400000000000000000000"),
        Int("185482642257398439114796845645546284380220968949399346684421580986889562184028199319100141244804501828416633516851200000000000000000000"),
        Int("16507955160908461081216919262453619309839666236496541854913520707833171034378509739399912570787600662729080382999756800000000000000000000"),
        Int("1485715964481761497309522733620825737885569961284688766942216863704985393094065876545992131370884059645617234469978112000000000000000000000"),
        Int("135200152767840296255166568759495142147586866476906677791741734597153670771559994765685283954750449427751168336768008192000000000000000000000"),
        Int("12438414054641307255475324325873553077577991715875414356840239582938137710983519518443046123837041347353107486982656753664000000000000000000000"),
        Int("1156772507081641574759205162306240436214753229576413535186142281213246807121467315215203289516844845303838996289387078090752000000000000000000000"),
        Int("108736615665674308027365285256786601004186803580182872307497374434045199869417927630229109214583415458560865651202385340530688000000000000000000000"),
        Int("10329978488239059262599702099394727095397746340117372869212250571234293987594703124871765375385424468563282236864226607350415360000000000000000000000"),
        Int("991677934870949689209571401541893801158183648651267795444376054838492222809091499987689476037000748982075094738965754305639874560000000000000000000000"),
        Int("96192759682482119853328425949563698712343813919172976158104477319333745612481875498805879175589072651261284189679678167647067832320000000000000000000000"),
        Int("9426890448883247745626185743057242473809693764078951663494238777294707070023223798882976159207729119823605850588608460429412647567360000000000000000000000"),
        Int("933262154439441526816992388562667004907159682643816214685929638952175999932299156089414639761565182862536979208272237582511852109168640000000000000000000000"),
        Int("93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000")
    };
    return table[n];
}

// ============================================================================
// factorialTable: table lookup (N <= 100)
// ============================================================================

const Int& IntCombinatorics::factorialTable(int n) {
    assert(n >= 0 && n <= TABLE_MAX);
    return getFactorialTable(n);
}

// ============================================================================
// factorialSimple: simple loop (sequential multiplication based on 100!)
// ============================================================================

Int IntCombinatorics::factorialSimple(int n) {
    if (n <= TABLE_MAX) return getFactorialTable(n);

    Int result = getFactorialTable(TABLE_MAX);
    for (int i = TABLE_MAX + 1; i <= n; i++) {
        result *= Int(i);
    }
    return result;
}

// ============================================================================
// factorialTournament: tournament method
// ============================================================================
//   Separate powers of 2 by LSB, then multiply the rest in a tournament fashion.
//   Multiply the ends together to balance the digit counts.
//   Finally apply 2^shift.

Int IntCombinatorics::factorialTournament(int n) {
    if (n <= 1) return Int::One();

    // Store, for each value 2..N, the value with powers of 2 removed
    std::vector<Int> x(n - 1);
    int shift = 0;
    for (int i = 2; i <= n; i++) {
        int p = std::countr_zero(static_cast<unsigned int>(i));
        x[i - 2] = Int(i >> p);
        shift += p;
    }

    int count = n - 1;

    // Multiply the ends together, halving the element count tournament-style
    while (count > 1) {
        int half = count / 2;
        for (int i = 0; i < half; i++) {
            x[i] *= x[count - 1 - i];
        }
        count = (count + 1) / 2;
    }

    return x[0] << shift;
}

// ============================================================================
// factorialSieve: sieve method (Eratosthenes-style prime factorization + tournament multiplication)
// ============================================================================
//   1. Separate the factor of 2 from each value 2..N (accumulate into shift)
//   2. From the remaining odd parts, extract prime powers Eratosthenes-sieve style
//   3. Compute each prime^exp
//   4. Combine via a balanced tree tournament
//   5. Finally apply the powers of 2 with << shift

Int IntCombinatorics::factorialSieve(int n) {
    if (n <= 1) return Int::One();

    // Step 1: separate powers of 2
    int shift = 0;
    std::vector<int> sieve(n + 1);
    for (int i = 2; i <= n; i++) {
        int val = i;
        int s = std::countr_zero(static_cast<unsigned int>(val));
        shift += s;
        sieve[i] = val >> s;
    }

    // Step 2: extract prime powers
    // Determine the array size from an approximation of PrimePi
    int nprime_est = (n < 55) ? 16 : static_cast<int>(n / (std::log(static_cast<double>(n)) - 4.0)) + 8;
    std::vector<int> primes;
    std::vector<int> exps;
    primes.reserve(nprime_est);
    exps.reserve(nprime_est);

    for (int i = 3; i <= n; i++) {
        if (sieve[i] != 1) {
            int prime = i;
            int count = 0;
            for (int j = i; j <= n; j += prime) {
                while (sieve[j] % prime == 0) {
                    sieve[j] /= prime;
                    count++;
                }
            }
            primes.push_back(prime);
            exps.push_back(count);
        }
    }

    if (primes.empty()) {
        return Int(1) << shift;
    }

    // Step 3: compute each prime^exp
    int pcount = static_cast<int>(primes.size());
    std::vector<Int> t(pcount);
    for (int i = 0; i < pcount; i++) {
        if (exps[i] == 1) {
            t[i] = Int(primes[i]);
        } else {
            t[i] = pow(Int(primes[i]), static_cast<unsigned int>(exps[i]));
        }
    }

    // Step 4: tournament multiplication (balance with end×end)
    while (pcount > 1) {
        int half = pcount / 2;
        for (int i = 0; i < half; i++) {
            t[i] *= t[pcount - 1 - i];
        }
        pcount = (pcount + 1) / 2;
    }

    return t[0] << shift;
}

// ============================================================================
// factorial: main entry point (automatic algorithm selection)
// ============================================================================

Int IntCombinatorics::factorial(const Int& n) {
    // Check for special states
    if (n.isNaN()) return Int::NaN();
    if (n.isInfinite()) {
        return (n.getSign() > 0) ? Int::PositiveInfinity() : Int::NaN();
    }
    if (n.isNegative()) return Int::NaN();
    if (n.isZero() || n.isOne()) return Int::One();

    // Convert to int
    if (!n.fitsInt64()) {
        // Very large n is impractical, but attempt it (sieve method)
        // However, n! large enough not to fit in int will run out of memory
        return Int::NaN();
    }
    int N = n.toInt();

    // Algorithm selection
    if (N <= TABLE_MAX) {
        return getFactorialTable(N);
    } else if (N <= SIMPLE_MAX) {
        return factorialSimple(N);
    } else {
        return factorialSieve(N);
    }
}

// ============================================================================
// Common divide-and-conquer product helper (shared by doubleFactorial / binomial / primorial)
// ============================================================================
namespace {
// Array version of balanced split product (assumes T = uint16_t / uint32_t / uint64_t):
// Keep Int construction cost low; at leaf nodes apply multiplyWord sequentially
template<typename T>
Int binarySplitProductArr(const T* vals, size_t lo, size_t hi) {
    if (lo >= hi) return Int::One();
    size_t len = hi - lo;
    if (len <= 8) {
        Int p(static_cast<uint64_t>(vals[lo]));
        for (size_t i = lo + 1; i < hi; i++) {
            IntOps::multiplyWord(p, static_cast<uint64_t>(vals[i]));
        }
        return p;
    }
    size_t mid = lo + len / 2;
    return binarySplitProductArr(vals, lo, mid) * binarySplitProductArr(vals, mid, hi);
}

// Product of consecutive integers: compute lo * (lo+1) * ... * hi (inclusive) by divide-and-conquer
// Even with single-precision inputs, the deeper recursion becomes multiplication of large Ints.
// Balancing maximizes the gain of Toom/NTT multiplication.
Int productRangeUL(uint64_t lo, uint64_t hi) {
    if (lo > hi) return Int::One();
    uint64_t len = hi - lo + 1;
    if (len <= 8) {
        Int p(lo);
        for (uint64_t i = lo + 1; i <= hi; i++) {
            IntOps::multiplyWord(p, i);
        }
        return p;
    }
    uint64_t mid = lo + len / 2;
    Int left = productRangeUL(lo, mid - 1);
    Int right = productRangeUL(mid, hi);
    return left * right;
}

// Product of odd numbers: lo * (lo+2) * (lo+4) * ... * hi (lo, hi are odd, lo <= hi)
Int productOdd(uint64_t lo, uint64_t hi) {
    if (lo > hi) return Int::One();
    uint64_t count = (hi - lo) / 2 + 1;
    if (count <= 8) {
        Int p(lo);
        for (uint64_t i = lo + 2; i <= hi; i += 2) {
            IntOps::multiplyWord(p, i);
        }
        return p;
    }
    uint64_t mid_count = count / 2;
    uint64_t mid_lo = lo + 2 * mid_count;          // smallest odd from the midpoint on
    uint64_t mid_hi = mid_lo - 2;                  // largest odd in the first half
    Int left = productOdd(lo, mid_hi);
    Int right = productOdd(mid_lo, hi);
    return left * right;
}
} // namespace

// ============================================================================
// doubleFactorial: double factorial n!!
// ============================================================================
// Even n: n!! = 2*4*6*...*n = 2^(n/2) * (n/2)!
// Odd n:  n!! = 1*3*5*...*n → product tree over odd numbers only

Int IntCombinatorics::doubleFactorial(const Int& n) {
    if (n.isNaN()) return Int::NaN();
    if (n.isInfinite()) {
        return (n.getSign() > 0) ? Int::PositiveInfinity() : Int::NaN();
    }
    if (n == -1) return Int::One();
    if (n < -1) return Int::NaN();
    if (n.isZero() || n.isOne()) return Int::One();

    // For huge values of n that do not fit in uint64, fall back to the naive loop
    if (n.bitLength() > 62) {
        Int result = Int::One();
        Int i = n;
        while (i > Int::One()) {
            result = result * i;
            i -= 2;
        }
        return result;
    }

    const uint64_t n_u = n.toUInt64();

    if ((n_u & 1) == 0) {
        // Even: n!! = 2^(n/2) * (n/2)!
        const uint64_t half = n_u / 2;
        Int fact;
        if (half <= 100) {
            fact = factorialTable(static_cast<int>(half));
        } else if (half <= 5000) {
            // Medium scale: divide-and-conquer product via productRangeUL. Avoids the
            // prime-power-extraction overhead of factorialSieve (vector<int> sieve, prime^exp
            // computation, tournament) and directly leverages the leaf multiplyWord chain.
            fact = productRangeUL(2, half);
        } else {
            // Large scale: factorialSieve (prime powers + tournament) is asymptotically favorable
            fact = factorialSieve(static_cast<int>(half));
        }
        return fact << static_cast<int>(half);
    }

    // Odd: 1*3*5*...*n by divide-and-conquer
    return productOdd(1, n_u);
}

// ============================================================================
// binomial: Kummer/Legendre prime-factorization helper
// ============================================================================
// v_p(n!) = Σ_{i≥1} floor(n/p^i)
// exp_p(C(n,k)) = v_p(n!) - v_p(k!) - v_p((n-k)!)
//
// By Kummer's theorem exp_p ≤ log_p(n), hence p^exp ≤ n, which always fits in uint64.
namespace {
inline uint64_t v_p_factorial(uint64_t n, uint64_t p) {
    uint64_t result = 0;
    uint64_t pk = p;
    while (pk <= n) {
        result += n / pk;
        if (pk > n / p) break;  // prevent overflow / exceeding n at the next pk*p
        pk *= p;
    }
    return result;
}

// Compute C(n, k) via Kummer/Legendre factorization.
// Substantially faster than a product tree for large n (e.g. n=100000).
// Avoids building huge intermediate values for numerator/denominator (~Stirling size); only needs
// prime-power tournament multiplication of the result size (~n bits at k=n/2).
Int binomialKummer(uint64_t n, uint64_t k_opt) {
    uint64_t nk = n - k_opt;  // = max(k, n-k)

    // Sieve of Eratosthenes (assumes n ≤ ~10^9; a bool vector is sufficient)
    std::vector<bool> is_prime(static_cast<size_t>(n + 1), true);
    is_prime[0] = is_prime[1] = false;
    for (uint64_t i = 2; i * i <= n; i++) {
        if (is_prime[i]) {
            for (uint64_t j = i * i; j <= n; j += i)
                is_prime[j] = false;
        }
    }

    // Collect prime powers p^exp_p(C(n,k))
    std::vector<uint64_t> prime_powers;
    prime_powers.reserve(static_cast<size_t>(static_cast<double>(n) /
                                              std::log(static_cast<double>(n) + 2.0)) + 16);

    for (uint64_t p = 2; p <= n; p++) {
        if (!is_prime[p]) continue;

        uint64_t exp = v_p_factorial(n, p) - v_p_factorial(k_opt, p) - v_p_factorial(nk, p);
        if (exp == 0) continue;

        if (exp == 1 || p > n / p) {
            // exp == 1 (most p > sqrt(n)), or p² > n always gives exp == 1
            prime_powers.push_back(p);
        } else {
            // Compute p^exp (exp ≤ log_p(n) gives p^exp ≤ n, uint64-safe)
            uint64_t val = 1;
            for (uint64_t i = 0; i < exp; i++) val *= p;
            prime_powers.push_back(val);
        }
    }

    if (prime_powers.empty()) return Int::One();

    // Balanced split tournament multiplication
    return binarySplitProductArr(prime_powers.data(), 0, prime_powers.size());
}
} // namespace

Int IntCombinatorics::binomial(const Int& n, const Int& k) {
    if (n.isNaN() || k.isNaN()) return Int::NaN();
    if (n.isInfinite() || k.isInfinite()) return Int::NaN();
    if (n.isNegative() || k.isNegative()) return Int::NaN();
    if (k > n) return Int::Zero();
    if (k.isZero()) return Int::One();
    if (k == n) return Int::One();

    Int k_opt = k;
    if ((k << 1) > n) {
        k_opt = n - k;
    }
    if (k_opt.isOne()) return n;

    // For huge inputs where n, k_opt do not fit in uint64, fall back to the naive loop
    // (in practice C(2^64, 2^63) is uncomputable on modern machines, so this is a guard)
    if (n.bitLength() > 62 || k_opt.bitLength() > 62) {
        Int result = Int::One();
        Int i = Int::One();
        while (i <= k_opt) {
            result = result * (n - i + Int::One());
            result = result / i;
            i = i + Int::One();
        }
        return result;
    }

    const uint64_t n_u = n.toUInt64();
    const uint64_t k_u = k_opt.toUInt64();

    // For small k_opt an in-place loop is fastest (reduces Int construction cost)
    // Boundary at 32: linear loop for 32 or below, Kummer factorization beyond
    // (measured at k=50: simple loop 1069ns vs Kummer 759ns, Kummer wins)
    if (k_u <= 32) {
        Int result(n_u);  // result = n
        for (uint64_t i = 2; i <= k_u; i++) {
            IntOps::multiplyWord(result, n_u - i + 1);
            IntOps::divExactWord(result, i);
        }
        return result;
    }

    // Large k_opt: compute directly via Kummer/Legendre prime factorization.
    // The old product-tree (building both numerator and denominator at Stirling size and dividing)
    // had intermediate values too large relative to the result; improved 52.9x → 2.23x at n=100000 (2026-05-08).
    return binomialKummer(n_u, k_u);
}

// =============================================================================
// Fibonacci / Lucas delegate to IntSequence::fibonacci2 (GMP-style square-only doubling).
// Uses the square-only form (2 sqr + add/subtract / iter).
// =============================================================================

Int IntCombinatorics::fibonacci(const Int& n) {
    if (n.isNaN()) return Int::NaN();
    if (n.isNegative()) return Int::NaN();
    if (n.isInfinite()) return Int::PositiveInfinity();
    if (n.isZero()) return Int(0);

    // GMP-style square-only doubling (via IntSequence): 2 sqr + add/subtract / iteration.
    // Theoretically ~1.8x faster than the old fibonacci_mpn (1 mul + 2 sqr / iter).
    return IntSequence::fibonacci(n.toUInt64());
}

// =============================================================================
// Lucas — L(n) = 2*F(n+1) - F(n)
// =============================================================================

Int IntCombinatorics::lucas(const Int& n) {
    if (n.isNaN()) return Int::NaN();
    if (n.isNegative()) return Int::NaN();
    if (n.isInfinite()) return Int::PositiveInfinity();
    if (n.isZero()) return Int(2);

    // Via GMP-style square-only doubling (IntSequence::lucas).
    // L(n) = F(n) + 2*F(n-1). Since fibonacci2 is square-only, this avoids general mul.
    return IntSequence::lucas(n.toUInt64());
}

// ============================================================================
// primorial: product of primes up to n
// ============================================================================

// Binary split product: recursively compute the product of vals[lo..hi)
// Balanced multiplication keeps size imbalance in check
static Int binarySplitProduct(const std::vector<Int>& vals, size_t lo, size_t hi) {
    if (lo >= hi) return Int::One();
    if (hi - lo == 1) return vals[lo];
    if (hi - lo == 2) return vals[lo] * vals[lo + 1];
    size_t mid = lo + (hi - lo) / 2;
    return binarySplitProduct(vals, lo, mid) * binarySplitProduct(vals, mid, hi);
}

// Small prime table (≤ 1000, 168 primes).
// Avoids sieve + vector allocation for primorial of small n.
namespace {
constexpr uint16_t SMALL_PRIMES[] = {
       2,    3,    5,    7,   11,   13,   17,   19,   23,   29,
      31,   37,   41,   43,   47,   53,   59,   61,   67,   71,
      73,   79,   83,   89,   97,  101,  103,  107,  109,  113,
     127,  131,  137,  139,  149,  151,  157,  163,  167,  173,
     179,  181,  191,  193,  197,  199,  211,  223,  227,  229,
     233,  239,  241,  251,  257,  263,  269,  271,  277,  281,
     283,  293,  307,  311,  313,  317,  331,  337,  347,  349,
     353,  359,  367,  373,  379,  383,  389,  397,  401,  409,
     419,  421,  431,  433,  439,  443,  449,  457,  461,  463,
     467,  479,  487,  491,  499,  503,  509,  521,  523,  541,
     547,  557,  563,  569,  571,  577,  587,  593,  599,  601,
     607,  613,  617,  619,  631,  641,  643,  647,  653,  659,
     661,  673,  677,  683,  691,  701,  709,  719,  727,  733,
     739,  743,  751,  757,  761,  769,  773,  787,  797,  809,
     811,  821,  823,  827,  829,  839,  853,  857,  859,  863,
     877,  881,  883,  887,  907,  911,  919,  929,  937,  941,
     947,  953,  967,  971,  977,  983,  991,  997
};
constexpr int SMALL_PRIMES_COUNT = static_cast<int>(sizeof(SMALL_PRIMES) / sizeof(SMALL_PRIMES[0]));
constexpr int SMALL_PRIMES_MAX = 1000;
} // namespace

Int IntCombinatorics::primorial(int n) {
    if (n < 2) return Int::One();

    // Small n: use the static prime table to completely avoid sieve + vector allocation
    if (n <= SMALL_PRIMES_MAX) {
        // Get the count of primes up to n by binary search
        int count = static_cast<int>(
            std::upper_bound(SMALL_PRIMES, SMALL_PRIMES + SMALL_PRIMES_COUNT,
                             static_cast<uint16_t>(n)) - SMALL_PRIMES);
        if (count == 0) return Int::One();
        return binarySplitProductArr(SMALL_PRIMES, 0, static_cast<size_t>(count));
    }

    // Large n: Sieve of Eratosthenes
    std::vector<bool> sieve(static_cast<size_t>(n + 1), true);
    sieve[0] = sieve[1] = false;
    for (int i = 2; static_cast<long long>(i) * i <= n; ++i) {
        if (sieve[i]) {
            for (int j = i * i; j <= n; j += i)
                sieve[j] = false;
        }
    }

    // Collect primes as uint32 (n is an int argument, so uint32 is sufficient)
    std::vector<uint32_t> primes;
    primes.reserve(static_cast<size_t>(static_cast<double>(n) / std::log(static_cast<double>(n) + 2.0)) + 8);
    for (int i = 2; i <= n; ++i) {
        if (sieve[i]) primes.push_back(static_cast<uint32_t>(i));
    }

    if (primes.empty()) return Int::One();

    // Binary split product (uint32 → Int converted in bulk at the leaf nodes)
    return binarySplitProductArr(primes.data(), 0, primes.size());
}

// ============================================================================
// mfac: m-factorial n·(n-m)·(n-2m)·...
// ============================================================================

Int IntCombinatorics::mfac(int n, int m) {
    if (m <= 0) return Int::NaN();
    if (n < 0) return Int::NaN();
    if (n <= 1) return Int::One();

    // Delegate special cases
    if (m == 1) return factorial(Int(n));
    if (m == 2) return doubleFactorial(Int(n));

    // Collect factors
    std::vector<Int> factors;
    for (int i = n; i >= 1; i -= m) {
        factors.emplace_back(i);
    }

    if (factors.empty()) return Int::One();

    // Binary split product
    return binarySplitProduct(factors, 0, factors.size());
}

} // namespace sangi