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

// IntCombinatorics.cpp
// 組合せ論的関数の実装
//
// 階乗アルゴリズム:
//   lib++20 からの移植 (4種のアルゴリズム)
//   参考文献:
//     http://www.luschny.de/math/factorial/FastFactorialFunctions.htm
//     http://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/1138-22.pdf (トーナメント法)
//     http://www.apfloat.org/ (バイナリ分割法)

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

namespace calx {

// ============================================================================
// 階乗テーブル (0! 〜 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: テーブル参照 (N <= 100)
// ============================================================================

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

// ============================================================================
// factorialSimple: 単純ループ (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: トーナメント法
// ============================================================================
// lib++20 からの移植:
//   2 の冪を LSB で分離し、残りをトーナメント方式で乗算。
//   端と端を掛け合わせて桁数のバランスを取る。
//   最後に 2^shift を適用。

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

    // 配列に 2..N の各値から 2 の冪を除去した値を設定
    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;

    // 端と端を掛けてトーナメント方式で要素数を半分にしてゆく
    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: 篩法 (エラトステネス風素因数分解 + トーナメント乗算)
// ============================================================================
// lib++20 からの移植:
//   1. 2..N の各値から 2 の因子を分離 (shift に積算)
//   2. 残りの奇数部分から、エラトステネスの篩風に素数の冪を抽出
//   3. 各 prime^exp を計算
//   4. バランス木トーナメントで合算
//   5. 最後に << shift で 2 の冪を適用

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

    // ステップ1: 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;
    }

    // ステップ2: 素数の冪を抽出
    // 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;
    }

    // ステップ3: 各 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]));
        }
    }

    // ステップ4: トーナメント乗算 (端×端でバランス)
    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: メインエントリーポイント (アルゴリズム自動選択)
// ============================================================================

Int IntCombinatorics::factorial(const Int& n) {
    // 特殊状態のチェック
    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();

    // int に変換
    if (!n.fitsInt64()) {
        // 非常に大きい n は実用的でないが、試みる (篩法)
        // ただし int に収まらないほど大きい n! はメモリ不足になる
        return Int::NaN();
    }
    int N = n.toInt();

    // アルゴリズム選択
    if (N <= TABLE_MAX) {
        return getFactorialTable(N);
    } else if (N <= SIMPLE_MAX) {
        return factorialSimple(N);
    } else {
        return factorialSieve(N);
    }
}

// ============================================================================
// doubleFactorial: 二重階乗 n!!
// ============================================================================

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 == Int(-1)) return Int::One();
    if (n < Int(-1)) return Int::NaN();
    if (n.isZero() || n.isOne()) return Int::One();

    Int result = Int::One();
    Int i = n;
    while (i > Int::One()) {
        result = result * i;
        i -= 2;
    }
    return result;
}

// ============================================================================
// binomial: 二項係数 C(n, k)
// ============================================================================

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 * Int(2) > n) {
        k_opt = n - k;
    }

    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;
}

// =============================================================================
// Fibonacci — 高速倍加法 (fast doubling)
// =============================================================================

// 内部: (F(n), F(n+1)) を返す
// mpn レベル fast doubling (反復版, Int オブジェクト構築なし)
// F(n) のビット長 ≈ 0.694 * n bits
// F(2k)   = F(k) * (2*F(k+1) - F(k))
// F(2k+1) = F(k)^2 + F(k+1)^2
static Int fibonacci_mpn(uint64_t n) {
    // F(93) = 12200160415121876738 は uint64_t に収まる最大の Fibonacci 数
    if (n <= 93) {
        static constexpr uint64_t table[] = {
            0,1,1,2,3,5,8,13,21,34,55,89,144,233,377,610,987,
            1597,2584,4181,6765,10946,17711,28657,46368,75025,
            121393,196418,317811,514229,832040,1346269,2178309,
            3524578,5702887,9227465,14930352,24157817,39088169,
            63245986,102334155,165580141,267914296,433494437,
            701408733,1134903170,1836311903,2971215073ULL,
            4807526976ULL,7778742049ULL,12586269025ULL,
            20365011074ULL,32951280099ULL,53316291173ULL,
            86267571272ULL,139583862445ULL,225851433717ULL,
            365435296162ULL,591286729879ULL,956722026041ULL,
            1548008755920ULL,2504730781961ULL,4052739537881ULL,
            6557470319842ULL,10610209857723ULL,17167680177565ULL,
            27777890035288ULL,44945570212853ULL,72723460248141ULL,
            117669030460994ULL,190392490709135ULL,308061521170129ULL,
            498454011879264ULL,806515533049393ULL,1304969544928657ULL,
            2111485077978050ULL,3416454622906707ULL,5527939700884757ULL,
            8944394323791464ULL,14472334024676221ULL,23416728348467685ULL,
            37889062373143906ULL,61305790721611591ULL,99194853094755497ULL,
            160500643816367088ULL,259695496911122585ULL,420196140727489673ULL,
            679891637638612258ULL,1100087778366101931ULL,
            1779979416004714189ULL,2880067194370816120ULL,
            4660046610375530309ULL,7540113804746346429ULL,
            12200160415121876738ULL,
        };
        return Int(table[n]);
    }
    // n > 93: mpn レベル fast doubling

    // F(n) の最大 limb 数を推定: ceil(0.694 * n / 64) + 余裕
    size_t max_limbs = static_cast<size_t>(n * 0.011 + 8); // 0.694/64 ≈ 0.0108, +8 for safety

    // バッファ確保: a=F(k), b=F(k+1), c=F(2k), d=F(2k+1), tmp (乗算結果)
    size_t buf_sz = max_limbs + 2; // carry/overflow 余裕
    size_t prod_sz = 2 * buf_sz + 4;
    size_t sq_scratch_sz = mpn::square_scratch_size(buf_sz);
    size_t mul_scratch_sz = mpn::multiply_scratch_size(buf_sz, buf_sz);
    size_t scratch_sz = std::max(sq_scratch_sz, mul_scratch_sz);
    // 安全余裕: scratch_sz に 16*buf_sz を最低保証
    if (scratch_sz < 16 * buf_sz) scratch_sz = 16 * buf_sz;

    // 全バッファを 1 回の確保で統合
    size_t total = buf_sz * 6 + prod_sz * 2 + scratch_sz;
    std::vector<uint64_t> arena(total, 0);
    uint64_t* ap = arena.data();

    uint64_t* a = ap; ap += buf_sz;      // F(k)
    uint64_t* b = ap; ap += buf_sz;      // F(k+1)
    uint64_t* c = ap; ap += buf_sz;      // temp for F(2k) or swap
    uint64_t* d = ap; ap += buf_sz;      // temp for F(2k+1) or swap
    uint64_t* t1 = ap; ap += buf_sz;     // temp
    uint64_t* prod = ap; ap += prod_sz;  // 乗算結果
    uint64_t* prod2 = ap; ap += prod_sz; // 乗算結果 2
    uint64_t* scratch = ap;              // multiply/square scratch

    // 初期値: F(0)=0, F(1)=1
    a[0] = 0; size_t an = 0; // F(k) = 0
    b[0] = 1; size_t bn = 1; // F(k+1) = 1

    // n のビットを MSB から走査
    int bits = 63 - std::countl_zero(n); // MSB のビット位置 (1-indexed skip)

    for (int i = bits; i >= 0; --i) {
        // doubling step:
        // c = F(2k)   = a * (2*b - a)
        // d = F(2k+1) = a^2 + b^2

        if (an == 0) {
            // F(k) = 0 → F(2k) = 0, F(2k+1) = F(k+1)^2
            // c = 0, d = b^2
            size_t cn = 0;
            std::memset(c, 0, buf_sz * sizeof(uint64_t));

            mpn::square(prod, b, bn, scratch);
            size_t dn = mpn::normalized_size(prod, 2 * bn);
            std::memcpy(d, prod, dn * sizeof(uint64_t));

            if (n & (1ULL << i)) {
                // F(2k+1), F(2k+2) = F(2k+1) + F(2k)
                // a = d, b = c + d = d (since c=0)
                std::memcpy(a, d, dn * sizeof(uint64_t)); an = dn;
                std::memcpy(b, d, dn * sizeof(uint64_t)); bn = dn;
            } else {
                // a = c = 0, b = d
                an = 0;
                std::memcpy(b, d, dn * sizeof(uint64_t)); bn = dn;
            }
            continue;
        }

        // t1 = 2*b (lshift 1)
        uint64_t overflow = mpn::lshift(t1, b, bn, 1);
        size_t t1n = bn;
        if (overflow) { t1[t1n] = overflow; t1n++; }

        // t1 = 2*b - a (t1 >= a は fast doubling の性質で保証)
        mpn::sub(t1, t1, t1n, a, an);
        t1n = mpn::normalized_size(t1, t1n);

        // c = a * t1 = F(k) * (2*F(k+1) - F(k)) = F(2k)
        if (an > 0 && t1n > 0) {
            // scratch_sz が 0 なら nullptr を渡す (NTT は内部確保)
            uint64_t* sc = scratch_sz > 0 ? scratch : nullptr;

            if (an >= t1n)
                mpn::multiply(prod, a, an, t1, t1n, sc);
            else
                mpn::multiply(prod, t1, t1n, a, an, sc);
            size_t cn = mpn::normalized_size(prod, an + t1n);
            std::memcpy(c, prod, cn * sizeof(uint64_t));

            // d = a^2 + b^2 = F(2k+1)
            mpn::square(prod, a, an, sc);
            size_t a2n = mpn::normalized_size(prod, 2 * an);
            mpn::square(prod2, b, bn, sc);
            size_t b2n = mpn::normalized_size(prod2, 2 * bn);

            size_t dn;
            if (a2n >= b2n) {
                mpn::add(d, prod, a2n, prod2, b2n);
                dn = a2n;
            } else {
                mpn::add(d, prod2, b2n, prod, a2n);
                dn = b2n;
            }
            if (d[dn]) dn++; // carry
            dn = mpn::normalized_size(d, dn + 1);

            if (n & (1ULL << i)) {
                // bit=1: a = d, b = c + d
                std::memcpy(a, d, dn * sizeof(uint64_t)); an = dn;
                // b = c + d
                if (cn >= dn) {
                    mpn::add(b, c, cn, d, dn);
                    bn = cn;
                } else {
                    mpn::add(b, d, dn, c, cn);
                    bn = dn;
                }
                if (b[bn]) bn++;
                bn = mpn::normalized_size(b, bn + 1);
            } else {
                // bit=0: a = c, b = d
                std::memcpy(a, c, cn * sizeof(uint64_t)); an = cn;
                std::memcpy(b, d, dn * sizeof(uint64_t)); bn = dn;
            }
        } else {
            // an == 0: handled above
        }
    }

    if (an == 0) return Int::Zero();
    return Int::fromRawWords(std::span<const uint64_t>(a, an), 1);
}

// 旧 fib_pair (互換性のため残す — lucas が使用)
static std::pair<Int, Int> fib_pair(uint64_t n) {
    if (n == 0) return {Int(0), Int(1)};

    auto [a, b] = fib_pair(n >> 1);
    Int c = a * ((b << 1) - a);
    Int d = a * a + b * b;

    if (n & 1) {
        return {d, c + d};
    } else {
        return {c, d};
    }
}

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);

    uint64_t idx = n.toUInt64();
    return fibonacci_mpn(idx);
}

// =============================================================================
// 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);

    uint64_t idx = n.toUInt64();
    auto [f, f1] = fib_pair(idx);
    // L(n) = 2*F(n+1) - F(n)
    return (f1 << 1) - f;
}

// ============================================================================
// primorial: n 以下の素数の積
// ============================================================================

// バイナリ分割積: vals[lo..hi) の積を再帰的に計算
// バランスの良い乗算で大きさの偏りを抑える
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);
}

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

    // エラトステネス篩
    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;
        }
    }

    // 素数を収集
    std::vector<Int> primes;
    for (int i = 2; i <= n; ++i) {
        if (sieve[i]) primes.emplace_back(i);
    }

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

    // バイナリ分割積
    return binarySplitProduct(primes, 0, primes.size());
}

// ============================================================================
// mfac: m-階乗 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();

    // 特殊ケース委譲
    if (m == 1) return factorial(Int(n));
    if (m == 2) return doubleFactorial(Int(n));

    // 因子を収集
    std::vector<Int> factors;
    for (int i = n; i >= 1; i -= m) {
        factors.emplace_back(i);
    }

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

    // バイナリ分割積
    return binarySplitProduct(factors, 0, factors.size());
}

} // namespace calx