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

// IntPrime.hpp
// 素数判定（Miller-Rabin、Fermat）for Int

#ifndef CALX_INT_PRIME_HPP
#define CALX_INT_PRIME_HPP

#include "IntBase.hpp"
#include <vector>
#include <utility>
#include <atomic>

namespace calx {

/**
 * @brief 素数判定アルゴリズム
 *
 * 確率的素数判定（Miller-Rabin）と決定的な小素数チェックを提供します。
 * RSA暗号、楕円曲線暗号、素因数分解などで使用されます。
 */
class IntPrime {
public:
    /**
     * @brief 確率的素数判定（Miller-Rabin テスト）
     *
     * Miller-Rabin アルゴリズムによる確率的素数判定を行います。
     * k 回のテストを行い、誤判定確率は最大 4^(-k) です。
     *
     * アルゴリズム:
     *   n-1 = 2^s × d (d は奇数) と分解
     *   k 回のラウンドで、底 a (2 <= a <= n-2) について:
     *   - a^d mod n が 1 なら次のラウンドへ
     *   - a^(2^r × d) mod n が -1 (r = 0..s-1) なら次のラウンドへ
     *   - どちらも満たさなければ合成数と確定
     *
     * @param n 判定する整数
     * @param k テスト回数（デフォルト 20, 誤判定確率 ≤ 2^(-40)）
     * @return true: おそらく素数, false: 確実に合成数
     *
     * 特殊ケース:
     *   isMillerRabinPrime(n < 2) = false
     *   isMillerRabinPrime(2) = true
     *   isMillerRabinPrime(偶数) = false
     *   isMillerRabinPrime(NaN) = false
     *   isMillerRabinPrime(∞) = false
     *
     * 注意:
     *   - k = 20 の場合、誤判定確率は 2^(-40) ≈ 10^(-12)
     *   - 暗号用途では k = 40 以上を推奨（誤判定確率 ≤ 2^(-80)）
     *   - 決定的判定ではないため、確実性が必要な場合は別のアルゴリズムを使用
     */
    static bool isMillerRabinPrime(const Int& n, int k = 20);

    /**
     * @brief フェルマーの小定理による確率的素数判定
     *
     * フェルマーの小定理: p が素数なら a^(p-1) ≡ 1 (mod p) (gcd(a,p)=1)
     * この逆を利用して素数判定を行います。
     *
     * 注意:
     *   - Carmichael 数（擬素数）を素数と誤判定する可能性があります
     *   - Miller-Rabin の方が信頼性が高いため、通常は Miller-Rabin を使用してください
     *   - 教育目的や簡易チェック用
     *
     * @param n 判定する整数
     * @param k テスト回数（デフォルト 5）
     * @return true: おそらく素数, false: 確実に合成数
     *
     * 特殊ケース:
     *   isFermatPrime(n < 2) = false
     *   isFermatPrime(2) = true
     *   isFermatPrime(偶数) = false
     *   isFermatPrime(NaN) = false
     */
    static bool isFermatPrime(const Int& n, int k = 5);

    /**
     * @brief 汎用素数判定（推奨）
     *
     * Miller-Rabin テストを使用した確率的素数判定です。
     * isMillerRabinPrime のエイリアスです。
     *
     * @param n 判定する整数
     * @param k テスト回数（デフォルト 20）
     * @return true: おそらく素数, false: 確実に合成数
     */
    static bool isProbablePrime(const Int& n, int k = 20) {
        return isMillerRabinPrime(n, k);
    }

    /**
     * @brief 決定的素数判定 (APRCL アルゴリズム)
     *
     * Adleman-Pomerance-Rumely-Cohen-Lenstra 法による決定的素数証明。
     * Miller-Rabin と異なり、確率的でなく数学的に正確な結果を返す。
     *
     * アルゴリズム:
     *   1. 円分環 Z[ζ_q]/(n) 上で Jacobi 和テストを実施
     *   2. 全テスト通過 → n の素因数は限定された剰余類に属する
     *   3. 最終因数探索で非自明因数がなければ素数確定
     *
     * @param n 判定する整数
     * @return true: 素数 (確定), false: 合成数 (確定)
     *
     * 注意:
     *   - Miller-Rabin より大幅に遅い (数秒〜数分)
     *   - 暗号鍵生成などの高速性が必要な場面では isProbablePrime を使用
     *   - 数学的証明や検証が必要な場面で使用
     */
    static bool isProvablePrime(const Int& n);

    /**
     * @brief n より大きい最小の素数を返す
     *
     * n の次の素数を探索します。n+1, n+2, ... と順に素数判定を行います。
     *
     * @param n 基準となる整数
     * @param k Miller-Rabin テストの回数（デフォルト 20）
     * @return n より大きい最小の素数
     *
     * 特殊ケース:
     *   nextPrime(n < 2) = 2
     *   nextPrime(2) = 3
     *   nextPrime(NaN) = NaN
     *   nextPrime(∞) = ∞
     *
     * 注意:
     *   - 大きい数では探索に時間がかかる可能性があります
     *   - 素数定理により、n 付近の素数間隔は約 ln(n) です
     */
    static Int nextPrime(const Int& n, int k = 20);

    /**
     * @brief n より小さい最大の素数を返す
     *
     * n の前の素数を探索します。n-1, n-2, ... と順に素数判定を行います。
     *
     * @param n 基準となる整数
     * @param k Miller-Rabin テストの回数（デフォルト 20）
     * @return n より小さい最大の素数
     *
     * 特殊ケース:
     *   prevPrime(n <= 2) = NaN（2 より小さい素数は存在しない）
     *   prevPrime(3) = 2
     *   prevPrime(NaN) = NaN
     *   prevPrime(∞) = NaN
     *
     * 注意:
     *   - 大きい数では探索に時間がかかる可能性があります
     */
    static Int prevPrime(const Int& n, int k = 20);

    /**
     * @brief Fermat 法による因数探索
     *
     * n = a² - b² = (a+b)(a-b) の形を探す。
     * n が 2 つの近い素数の積の場合に高速。
     *
     * @param n 合成数（奇数、n > 1）
     * @param max_iterations 最大反復回数（デフォルト 10000）
     * @return n の非自明な因数。見つからなければ n を返す
     */
    static Int findFactor_Fermat(const Int& n, uint64_t max_iterations = 10000);

    /**
     * @brief Pollard's rho 法による因数探索
     *
     * 疑似乱数列 x_{k+1} = x_k² + c (mod n) を用いた因数探索。
     * Brent の改良版（サイクル検出）を使用。
     * 期待計算量: O(n^{1/4})
     *
     * @param n 合成数（n > 1）
     * @return n の非自明な因数。見つからなければ n を返す
     */
    static Int findFactor_PollardRho(const Int& n);

    /**
     * @brief Pollard p-1 法による因数探索
     *
     * p-1 が B-smooth（全素因数 ≤ B）な素因数 p を見つける。
     * 移植元: lib++20 Factor_P1
     *
     * アルゴリズム:
     *   M = lcm(2, 3, ..., B) を段階的に計算し、
     *   gcd(a^M - 1, n) を求める。p | n かつ (p-1) | M なら
     *   a^M ≡ 1 (mod p) となり gcd > 1 で因数が見つかる。
     *
     * @param n 合成数（n > 1）
     * @param B1 Stage 1 の smooth bound（デフォルト 10000）
     * @return n の非自明な因数。見つからなければ n を返す
     */
    static Int findFactor_PollardP1(const Int& n, uint64_t B1 = 10000);

    /**
     * @brief ECM (楕円曲線法) による因数探索
     *
     * Lenstra の楕円曲線素因数分解法。
     * Montgomery 形式の楕円曲線 By² = x³ + Ax² + x を使用。
     * 20〜60 桁程度の因数の探索に有効。
     * 期待計算量: O(exp(√(2 ln p ln ln p))) (p は最小素因数)
     *
     * @param n 合成数（n > 1）
     * @param curves 試行する曲線の数（デフォルト 100）
     * @param B1 Stage 1 の smooth bound（デフォルト 10000）
     * @param B2 Stage 2 の bound（デフォルト 100*B1）
     * @return n の非自明な因数。見つからなければ n を返す
     */
    static Int findFactor_ECM(const Int& n, int curves = 100,
                              uint64_t B1 = 10000, uint64_t B2 = 0,
                              const std::atomic<bool>* cancel = nullptr);

    /**
     * @brief MPQS (多重多項式二次篩) による因数探索
     *
     * 二次篩法の多重多項式版。50〜100 桁の合成数の素因数分解に有効。
     * Q(x) = (ax + b)² - n の値が factor base 上で smooth になる
     * x の組を集め、線形代数 (GF(2) 上のガウス消去) で因数を導出。
     *
     * @param n 合成数（奇数, n > 1）
     * @return n の非自明な因数。見つからなければ n を返す
     */
    static Int findFactor_MPQS(const Int& n,
                               const std::atomic<bool>* cancel = nullptr);

    /**
     * @brief SIQS (Self-Initializing Quadratic Sieve) による因数分解
     *
     * MPQS の強化版。主な改良点:
     *   1. Knuth-Schroeppel 乗数選択
     *   2. sieve array による高速 smooth 判定 (log₂(p) 加算方式)
     *   3. SIQS 多項式: A = 複数 FB 素数の積, B を CRT で構成
     *   4. Gray code による B 切り替え (self-initialization)
     *   5. 符号列 (-1) の正しい処理
     *
     * @param n 合成数（奇数, n > 1）
     * @return n の非自明な因数。見つからなければ n を返す
     */
    static Int findFactor_SIQS(const Int& n,
                               const std::atomic<bool>* cancel = nullptr);

    /**
     * @brief GNFS (一般数体篩法) による因数探索
     *
     * 大きな合成数（60桁以上）に対する準指数時間アルゴリズム。
     * base-m 多項式選択、有理・代数因数基底、line sieve、
     * GF(2) 線形代数、Hensel lifting による代数平方根を実装。
     *
     * @param n 合成数（奇数, n > 1）
     * @return n の非自明な因数。見つからなければ n を返す
     */
    static Int findFactor_GNFS(const Int& n);

    /**
     * @brief Williams p+1 法による因数探索
     *
     * p+1 が B-smooth な素因数 p を見つける。p-1 法の双対。
     * Lucas 列 V_k(a) mod n を用いて gcd(V_M - 2, n) を計算。
     *
     * アルゴリズム:
     *   1. Lucas 列 V_k: V_0=2, V_1=a, V_{k+1}=a*V_k - V_{k-1}
     *   2. Binary ladder で V_M mod n を高速計算
     *   3. g = gcd(V_M - 2, n) で因数検出
     *
     * @param n 合成数（奇数, n > 1）
     * @param B1 smooth bound（デフォルト 10000）
     * @return n の非自明な因数。見つからなければ n を返す
     */
    static Int findFactor_WilliamsP1(const Int& n, uint64_t B1 = 10000);

    /**
     * @brief SQUFOF (Square Form Factorization) による因数探索
     *
     * Shanks の平方形式因数分解法。連分数展開に基づく。
     * O(n^{1/4}) 時間、O(1) メモリ。小〜中規模 (〜18桁) で高速。
     *
     * @param n 合成数（奇数, n > 1, 完全平方数でない）
     * @return n の非自明な因数。見つからなければ n を返す
     */
    static Int findFactor_SQUFOF(const Int& n);

    /**
     * @brief Pollard rho 法（Montgomery 乗算版）による因数探索
     *
     * Montgomery 剰余乗算を用いた高速な Pollard rho 法。
     * 移植元: lib++20 Factor_Rho_Montgomery
     *
     * アルゴリズム:
     *   R = 2^k > n となる R を取り、Montgomery 表現に変換。
     *   x_{k+1} = x_k² + 1 (mod n) を Montgomery 乗算で計算。
     *   Floyd のサイクル検出で gcd(x-y, n) > 1 となるまで反復。
     *
     * @param n 合成数（奇数, n > 1）
     * @return n の非自明な因数。見つからなければ n を返す
     */
    static Int findFactor_RhoMontgomery(const Int& n);

    /**
     * @brief Goldbach 予想に基づく因数探索
     *
     * N = p * q (半素数) に対し、sum ≈ p + q を二分探索で推定し、
     * 二次方程式 t² - sum*t + N = 0 を解いて因数を求める。
     * 移植元: lib++20 Factor_Goldbach (Dr. Roger G. Doss, PhD)
     *
     * 注意: Factor_Fermat より一般に遅い。教育目的で提供。
     *
     * @param n 合成数（n > 1, 半素数を想定）
     * @return n の非自明な因数。見つからなければ n を返す
     */
    static Int findFactor_Goldbach(const Int& n);

    /**
     * @brief 試し割り法による因数探索
     *
     * SMALL_PRIMES テーブルの素数で順に試し割りを行う。
     * 移植元: lib++20 Factor_Try
     *
     * @param n 合成数（n > 1）
     * @return n の非自明な因数。見つからなければ n を返す
     */
    static Int findFactor_TrialDivision(const Int& n);

    /**
     * @brief n 以下の素数の個数 (素数計数関数 π(n))
     *
     * エラトステネスの篩を用いて n 以下の素数を数えます。
     *
     * @param n 上限値（n >= 0）
     * @return n 以下の素数の個数
     *
     * 特殊ケース:
     *   primePi(0) = 0, primePi(1) = 0, primePi(2) = 1
     */
    static int64_t primePi(int64_t n);

    // 小素数で割り切れる因子をすべて抽出する（1リム・バッチ GCD 版）
    // remaining から小素数の因子を除去し、見つかった因子を factors に追加
    // remaining は小素数で割り切った後の値に更新される
    // 注意: 素因数 2 は事前に countTrailingZeros() で除去済みの前提
    static void extractSmallPrimeFactors(
        Int& remaining,
        std::vector<std::pair<Int, int>>& factors);

private:
    // 小素数リスト（3〜1987 の奇素数 299 個）
    // 素因数 2 はバッチ GCD ではなく LSB で処理するためここには含めない
    static constexpr int SMALL_PRIMES[] = {
        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,1009,1013,
        1019,1021,1031,1033,1039,1049,1051,1061,1063,1069,
        1087,1091,1093,1097,1103,1109,1117,1123,1129,1151,
        1153,1163,1171,1181,1187,1193,1201,1213,1217,1223,
        1229,1231,1237,1249,1259,1277,1279,1283,1289,1291,
        1297,1301,1303,1307,1319,1321,1327,1361,1367,1373,
        1381,1399,1409,1423,1427,1429,1433,1439,1447,1451,
        1453,1459,1471,1481,1483,1487,1489,1493,1499,1511,
        1523,1531,1543,1549,1553,1559,1567,1571,1579,1583,
        1597,1601,1607,1609,1613,1619,1621,1627,1637,1657,
        1663,1667,1669,1693,1697,1699,1709,1721,1723,1733,
        1741,1747,1753,1759,1777,1783,1787,1789,1801,1811,
        1823,1831,1847,1861,1867,1871,1873,1877,1879,1889,
        1901,1907,1913,1931,1933,1949,1951,1973,1979,1987
    };
    static constexpr int SMALL_PRIMES_COUNT = sizeof(SMALL_PRIMES) / sizeof(SMALL_PRIMES[0]);

    // 1リム・バッチ GCD 用テーブル
    // 連続する奇素数の積を 64bit に詰めたグループ
    struct PrimeBatch {
        uint64_t product;       // バッチ内の素数の積
        int      start_index;   // SMALL_PRIMES[] 内の開始インデックス
        int      count;         // バッチ内の素数の個数
    };

    // バッチテーブルを取得（初回呼び出し時に構築）
    static const std::vector<PrimeBatch>& getBatches();

    // 小素数での割り切れチェック（1リム・バッチ GCD 版）
    static bool isDivisibleBySmallPrime(const Int& n);
};

} // namespace calx

#endif // CALX_INT_PRIME_HPP