一、整除、素数与最大公约数

这是数论大厦的地基，看似简单，却蕴含着后续一切知识的基础。

1. 整除性 (Divisibility)

如果整数 a 除以非零整数 b，商为整数 q 且余数为 0，我们就说 b 整除 a，或者 a 是 b 的倍数，b 是 a 的约数（或因子）。记作 b | a。

性质：

• 自反性：a | a
• 传递性：若 a | b 且 b | c，则 a | c
• 若 a | b 且 a | c，则对于任意整数 x, y，有 a | (bx + cy)

这些性质在推导和证明中非常有用。

2. 素数与合数 (Primes and Composites)

• 素数 (Prime): 一个大于 1 的正整数 p，如果它的正约数只有 1 和它本身，则称 p 为素数。例如：2, 3, 5, 7, 11, ...
• 合数 (Composite): 一个大于 1 的正整数，如果不是素数，则称为合数。例如：4, 6, 8, 9, 10, ...
• 1: 既不是素数，也不是合数。

算术基本定理 (Fundamental Theorem of Arithmetic): 任何一个大于 1 的整数 n 都可以唯一地分解成若干个素数的乘积，形式如下：

$$n = p_1^{a_1} p_2^{a_2} \cdots p_k^{a_k}$$

其中 $p_1 < p_2 < \cdots < p_k$ 是素数， $a_i \ge 1$ 是整数。这被称为 n 的标准分解式。

这个定理是数论的基石之一。很多数论问题最终都归结为对素数幂的操作。

如何判断一个数 n 是否为素数？ 最朴素的方法是试除法 (Trial Division)：检查从 2 到 $\sqrt{n}$ 的所有整数是否能整除 n。如果都不能，则 n 是素数。

#include <bits/stdc++.h>

using namespace std;
using ll = long long; // 使用 ll 代替 long long

// 判断 n 是否为素数 (试除法)
// 时间复杂度: O(sqrt(n))
bool is_prime(int n) {
    if (n < 2) return false;
    // 注意 i*i 可能溢出 int，但这里 n 是 int，所以 i*i <= n <= INT_MAX，通常不会溢出
    // 如果 n 是 ll 类型，则应写为 for (ll i = 2; i * i <= n; ++i)
    for (int i = 2; i * i <= n; ++i) {
        if (n % i == 0) {
            return false;
        }
    }
    return true;
}

int main() {
    ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
    int n;
    cin >> n;
    if (is_prime(n)) {
        cout << n << " is prime\n";
    } else {
        cout << n << " is not prime\n";
    }
    return 0;
}

如何分解质因数？ 同样使用试除法，从 2 开始尝试除 n，如果 i 能整除 n，就不断除以 i 并记录 i 的次数，直到 n 不能再被 i 整除。然后尝试下一个数 i+1。持续这个过程直到 i * i > n。如果此时 n 仍然大于 1，那么剩下的 n 本身也是一个素因子。

#include <bits/stdc++.h>

using namespace std;
using ll = long long; // 使用 ll 代替 long long

// 分解 n 的质因数
// 返回一个 map，键是质因子，值是其指数
// 时间复杂度: O(sqrt(n))
map<int, int> prime_factorize(int n) {
    map<int, int> factors;
    if (n < 2) return factors;
    // 同样，若 n 是 ll 类型，i 也应为 ll
    for (int i = 2; i * i <= n; ++i) {
        if (n % i == 0) {
            int count = 0;
            while (n % i == 0) {
                n /= i;
                count++;
            }
            factors[i] = count;
        }
    }
    if (n > 1) { // 如果 n 还有一个大于 sqrt(n) 的质因子
        factors[n] = 1;
    }
    return factors;
}

int main() {
    ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
    int n;
    cin >> n;
    map<int, int> factors = prime_factorize(n);
    cout << n << " = ";
    bool first = true;
    for (auto const& [p, a] : factors) {
        if (!first) {
            cout << " * ";
        }
        cout << p;
        if (a > 1) {
            cout << "^" << a;
        }
        first = false;
    }
    cout << "\n";
    return 0;
}

3. 最大公约数 (GCD) 与 最小公倍数 (LCM)

• 最大公约数 (Greatest Common Divisor, GCD): 两个或多个整数共有的约数中最大的一个。记作 gcd(a, b) 或 (a, b)。
• 最小公倍数 (Least Common Multiple, LCM): 两个或多个整数共有的倍数中最小的一个正整数。记作 lcm(a, b) 或 [a, b]。

重要关系: 对于两个正整数 a, b，有：

$$a \times b = \gcd(a, b) \times \lcm(a, b)$$

这个关系非常有用，通常我们计算出 GCD 后，就可以直接用它来计算 LCM，避免了大数乘法溢出的风险（先除后乘）：lcm(a, b) = (a / gcd(a, b)) * b。

如何计算 GCD？ 欧几里得算法 (Euclidean Algorithm)，也叫辗转相除法。这是数论中最古老、最优美、也是最重要的算法之一。

原理: 基于定理 gcd(a, b) = gcd(b, a % b) (假设 a >= b > 0)。递归下去，直到 b = 0，此时 gcd(a, 0) = a。

#include <bits/stdc++.h>

using namespace std;
using ll = long long; // 使用 ll 代替 long long

// 欧几里得算法计算 gcd(a, b)
// 时间复杂度: O(log(min(a, b)))
ll gcd(ll a, ll b) {
    return b == 0 ? a : gcd(b, a % b);
}

// 计算 lcm(a, b)
// 注意可能溢出，使用 ll，并且先除后乘
ll lcm(ll a, ll b) {
    if (a == 0 || b == 0) return 0; // 处理 0 的情况
    // 防止 a/gcd(a,b) * b 溢出，即使 a, b, gcd 都是 ll
    // 可以写成 a * (b / gcd(a, b))，或者确保结果不超 ll 范围
    // 在竞赛环境下，通常假设结果在 ll 内，直接写 (a / gcd(a, b)) * b
    return (a / gcd(a, b)) * b;
}

int main() {
    ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
    ll a, b;
    cin >> a >> b;
    cout << "gcd(" << a << ", " << b << ") = " << gcd(a, b) << "\n";
    cout << "lcm(" << a << ", " << b << ") = " << lcm(a, b) << "\n";
    return 0;
}

复杂度分析: 欧几里得算法的复杂度是对数级别的，非常高效。每一次递归，至少会使其中一个数减半（斐波那契数列是其最坏情况）。

扩展欧几里得算法 (Extended Euclidean Algorithm - ExGCD)

欧几里得算法不仅能求出 gcd(a, b)，还能找到一对整数 x, y，使得：

$$ax + by = \gcd(a, b)$$

这就是贝祖定理 (Bézout's Identity)。ExGCD 就是求解这对 x, y 的算法。

原理: 基于 gcd(a, b) = gcd(b, a % b) 的递归过程。 假设我们已经求得 b 和 a % b 的一组解 x' 和 y'，满足：

$$bx' + (a \pmod b)y' = \gcd(b, a \pmod b)$$

我们知道 a % b = a - floor(a / b) * b，代入上式：

$$bx' + (a - \lfloor a/b \rfloor \times b) y' = \gcd(b, a \pmod b) $$

整理得：

$$ay' + b(x' - \lfloor a/b \rfloor y') = \gcd(b, a \pmod b) $$

由于 gcd(a, b) = gcd(b, a % b)，我们令：

$$x = y'$$ $$y = x' - \lfloor a/b \rfloor y'$$

这就找到了满足 ax + by = gcd(a, b) 的一组解 (x, y)。递归的边界是当 b = 0 时，gcd(a, 0) = a，此时 ax + 0y = a，解为 x = 1, y = 0。

#include <bits/stdc++.h>

using namespace std;
using ll = long long; // 使用 ll 代替 long long

// 扩展欧几里得算法
// 求解 ax + by = gcd(a, b) 的一组解 (x, y)
// 返回 gcd(a, b)，并将解存储在 x, y 引用的变量中
// 时间复杂度: O(log(min(a, b)))
ll exgcd(ll a, ll b, ll &x, ll &y) {
    if (b == 0) {
        x = 1;
        y = 0;
        return a;
    }
    ll d = exgcd(b, a % b, x, y);
    ll temp = x;
    x = y;
    y = temp - (a / b) * y;
    return d;
}

int main() {
    ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
    ll a, b, x, y;
    cin >> a >> b;
    ll d = exgcd(a, b, x, y);
    cout << "gcd(" << a << ", " << b << ") = " << d << "\n";
    cout << "A solution for " << a << "x + " << b << "y = " << d << " is: ";
    cout << "x = " << x << ", y = " << y << "\n";
    // 验证
    cout << "Verification: " << a << "*" << x << " + " << b << "*" << y << " = " << a * x + b * y << "\n";

    // 求解线性丢番图方程 ax + by = c
    ll c;
    cout << "Enter c for equation " << a << "x + " << b << "y = c: ";
    cin >> c;
    if (c % d != 0) {
        cout << "No integer solution exists.\n";
    } else {
        // 找到 ax + by = d 的一组解 x0, y0 (即上面求出的 x, y)
        // 则 ax + by = c 的一组特解是 x1 = x0 * (c/d), y1 = y0 * (c/d)
        // 注意 x * (c / d) 可能溢出 ll，如果 a, b, c 很大
        // 竞赛中通常假设不会爆 ll，或者题目会有说明
        ll factor = c / d;
        ll x1 = x * factor;
        ll y1 = y * factor;
        cout << "One particular solution for " << a << "x + " << b << "y = " << c << " is: ";
        cout << "x = " << x1 << ", y = " << y1 << "\n";

        // 通解形式： x = x1 + k * (b/d), y = y1 - k * (a/d) for integer k
        ll b_d = b / d;
        ll a_d = a / d;
        cout << "General solution form: x = " << x1 << " + k*(" << b_d << "), y = " << y1 << " - k*(" << a_d << ")\n";
        // 求最小正整数解 x
        // k = (x1 % b_d - x1) / b_d; // 这样计算 k 不太稳健
        // 最小非负 x 解为 (x1 % b_d + b_d) % b_d
        ll min_non_neg_x = (x1 % b_d + b_d) % b_d;
        // 如果要求最小正整数解
        ll min_pos_x = min_non_neg_x;
        if (min_pos_x == 0) {
            min_pos_x += abs(b_d); // 加 b/d 的绝对值，因为 b/d 可能为负
        }
        cout << "Smallest positive integer solution for x (if exists): " << min_pos_x << "\n";
    }

    return 0;
}

ExGCD 的主要应用：

1. 求解线性丢番图方程 (Linear Diophantine Equation): 形如 ax + by = c 的方程。该方程有整数解当且仅当 gcd(a, b) | c。若有解，可用 ExGCD 求出 ax + by = gcd(a, b) 的一组解 (x_0, y_0)，则 ax + by = c 的一组特解为 (x_1, y_1) = (x_0 * (c/d), y_0 * (c/d))，其中 d = gcd(a, b)。通解为 x = x_1 + k * (b/d), y = y_1 - k * (a/d)，其中 k 为任意整数。
2. 求解模线性方程 (Modular Linear Equation) / 求解模逆元 (Modular Inverse): 这是 ExGCD 在竞赛中最常见的应用，我们将在下一节详细讨论。

二、模运算的艺术：同余与模逆元

模运算 (Modulo Operation) 是数论在算法竞赛中的核心舞台。很多问题需要处理非常大的数，但我们往往只关心它们对某个数取模后的结果。

1. 同余 (Congruence)

如果两个整数 a 和 b 除以一个正整数 m 所得的余数相同，则称 a 和 b 模 m 同余。记作：

$$a \equiv b \pmod m$$

这等价于 m | (a - b)。m 称为模数 (Modulus)。

同余的性质 (模 m 下)：

• 自反性：a ≡ a
• 对称性：若 a ≡ b，则 b ≡ a
• 传递性：若 a ≡ b 且 b ≡ c，则 a ≡ c
• 和差性： 若 a ≡ b 且 c ≡ d，则 a + c ≡ b + d 且 a - c ≡ b - d
• 积性： 若 a ≡ b 且 c ≡ d，则 ac ≡ bd
• 幂性： 若 a ≡ b，则 a^k ≡ b^k (k 为非负整数)
• 若 ac ≡ bc \pmod m 且 gcd(c, m) = d，则 a ≡ b \pmod{m/d}。特别地，若 gcd(c, m) = 1，则 a ≡ b \pmod m (模意义下的除法需要小心)。

这些性质意味着我们可以在运算过程中随时取模，而不影响最终模 m 的结果（除法除外）。这对于防止中间结果溢出至关重要。

// C++ 中的取模运算注意事项
#include <iostream>
using ll = long long; // 使用 ll 代替 long long

int main() {
    ll a = -5;
    ll m = 3;
    // C++ 的 % 运算符对于负数的结果可能也是负数 (-5 % 3 = -2)
    ll cpp_remainder = a % m;
    std::cout << "C++ %: " << cpp_remainder << std::endl; // Output: -2 (或依赖编译器)

    // 在数论中，我们通常希望余数是非负的 (介于 0 和 m-1 之间)
    ll positive_remainder = (a % m + m) % m; // 结果为 1
    std::cout << "Positive remainder: " << positive_remainder << std::endl; // Output: 1
    return 0;
}

2. 模逆元 (Modular Inverse)

在实数运算中，a 的逆元是 1/a，因为 a * (1/a) = 1。模运算中也有类似的概念。 如果存在一个整数 x 使得：

$$ax \equiv 1 \pmod m$$

则称 x 是 a 模 m 的乘法逆元。通常记作 $a^{-1} \pmod m$。

模逆元存在的条件： a 模 m 的逆元存在 当且仅当 gcd(a, m) = 1。即 a 和 m 互素。

为什么需要模逆元？ 模运算对于加、减、乘是封闭的，但没有直接的“模除法”。如果你想计算 (a / b) % m，你不能简单地计算 (a % m) / (b % m) % m。正确的方法是找到 b 模 m 的逆元 b^{-1}，然后计算：

$$(a / b) \pmod m \equiv a \times b^{-1} \pmod m$$

这在处理组合计数问题（如计算组合数 $C(n, k) \pmod m$）时非常常见，因为组合数公式涉及除法。

如何求解模逆元？

方法一：使用扩展欧几里得算法 (ExGCD) 求解 ax \equiv 1 \pmod m 等价于求解线性丢番图方程 ax + my = 1。 根据 ExGCD，这个方程有解当且仅当 gcd(a, m) | 1，即 gcd(a, m) = 1。 当 gcd(a, m) = 1 时，我们使用 ExGCD 找到 x, y 使得 ax + my = 1。那么这个 x 就是 a 模 m 的一个逆元。因为我们通常需要一个 [0, m-1] 范围内的逆元，所以最终结果是 (x % m + m) % m。

#include <bits/stdc++.h>

using namespace std;
using ll = long long; // 使用 ll 代替 long long

// ExGCD (复用之前的代码)
ll exgcd(ll a, ll b, ll &x, ll &y) {
    if (b == 0) { x = 1; y = 0; return a; }
    ll d = exgcd(b, a % b, x, y);
    ll temp = x; x = y; y = temp - (a / b) * y;
    return d;
}

// 求解 a 模 m 的逆元
// 返回 -1 表示逆元不存在
// 时间复杂度: O(log(min(a, m)))
ll modInverse_exgcd(ll a, ll m) {
    ll x, y;
    ll d = exgcd(a, m, x, y);
    if (d != 1) {
        return -1; // 逆元不存在
    }
    // 确保逆元是正数且在 [0, m-1] 范围内
    return (x % m + m) % m;
}

int main() {
    ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
    ll a = 7;
    ll m = 26; // gcd(7, 26) = 1
    ll inv = modInverse_exgcd(a, m);
    if (inv != -1) {
        cout << "Inverse of " << a << " mod " << m << " is " << inv << "\n"; // inv = 15
        cout << "Verification: (" << a << " * " << inv << ") % " << m << " = " << (a * inv) % m << "\n"; // 7 * 15 = 105. 105 % 26 = 1.
    } else {
        cout << "Inverse of " << a << " mod " << m << " does not exist.\n";
    }

    a = 6;
    m = 26; // gcd(6, 26) = 2
    inv = modInverse_exgcd(a, m);
     if (inv != -1) {
        cout << "Inverse of " << a << " mod " << m << " is " << inv << "\n";
    } else {
        cout << "Inverse of " << a << " mod " << m << " does not exist.\n"; // Correct
    }
    return 0;
}

方法二：费马小定理 (Fermat's Little Theorem) 条件： 模数 m 必须是素数，并且 a 不能是 m 的倍数 (a % m != 0)。 定理内容： 若 p 是素数，且 p 不能整除 a，则：

$$a^{p-1} \equiv 1 \pmod p$$

推论求逆元： 将上式两边同乘以 $a^{-1}$ (假设存在)，得到：

$$a^{p-2} \equiv a^{-1} \pmod p$$

所以，当 p 是素数时，a 模 p 的逆元就是 $a^{p-2} \pmod p$。我们可以使用快速幂 (Modular Exponentiation) 算法来高效计算 $a^{p-2} \pmod p$。

快速幂算法 (Modular Exponentiation / Exponentiation by Squaring) 用于快速计算 $a^b \pmod m$。其思想是将指数 b 表示为二进制形式，然后利用 $a^{2k} = (a^k)^2$ 和 $a^{k+1} = a^k \times a$ 来加速计算。

#include <bits/stdc++.h>

using namespace std;
using ll = long long; // 使用 ll 代替 long long

// 快速幂算法 计算 (base^exp) % mod
// 时间复杂度: O(log(exp))
ll qpow(ll base, ll exp, ll mod) {
    ll res = 1;
    base %= mod; // 预处理，防止 base 过大或为负
    if (base < 0) base += mod; // 确保 base 非负
    while (exp > 0) {
        if (exp % 2 == 1) { // 如果当前位是 1
            res = (res * base) % mod;
        }
        base = (base * base) % mod; // base 平方
        exp /= 2; // 指数右移一位
    }
    return res;
}

// 使用费马小定理求逆元 (要求 mod 是素数)
// 时间复杂度: O(log(mod))
ll modInverse_fermat(ll a, ll p) {
    if (p <= 1) return -1; // 素数必须大于 1
    // 这里可以加一个素性测试，但通常题目会保证 p 是素数
    a %= p;
    if (a < 0) a += p;
    if (a == 0) return -1; // a 是 p 的倍数，无逆元
    // gcd(a, p) != 1 时无逆元，这已包含在 a=0 的情况里，因为 p 是素数

    // 指数是 p - 2
    if (p - 2 < 0) { // 处理 p=2 的情况
        if (p == 2 && a == 1) return 1; // 1 mod 2 的逆元是 1
        else return -1; // 其他情况 p=2 时指数为 0 或负数，或 a=0 mod 2
    }
    return qpow(a, p - 2, p);
}

int main() {
    ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
    ll a = 7;
    ll p = 29; // 29 is prime
    ll inv = modInverse_fermat(a, p);
     if (inv != -1) {
        cout << "Inverse of " << a << " mod " << p << " is " << inv << "\n"; // inv = 25
        cout << "Verification: (" << a << " * " << inv << ") % " << p << " = " << (a * inv) % p << "\n"; // 7 * 25 = 175. 175 = 6*29 + 1. 175 % 29 = 1.
    } else {
        cout << "Inverse of " << a << " mod " << p << " does not exist or p is not prime.\n";
    }

    a = 3;
    p = 2; // 2 is prime
    inv = modInverse_fermat(a, p); // a=3 % 2 = 1. p-2=0. 1^0 = 1.
     if (inv != -1) {
        cout << "Inverse of " << a << " mod " << p << " is " << inv << "\n"; // inv = 1
        cout << "Verification: (" << a << " * " << inv << ") % " << p << " = " << (a * inv) % p << "\n"; // (3*1)%2 = 1
    } else {
         cout << "Inverse of " << a << " mod " << p << " does not exist or p is not prime.\n";
    }

    return 0;
}

ExGCD vs. 费马小定理：

• ExGCD:
  • 优点：适用范围更广，只要 gcd(a, m) = 1 即可，m 不必是素数。
  • 缺点：实现稍复杂一点点。
• 费马小定理:
  • 优点：实现简单，只需要快速幂。
  • 缺点：强烈依赖 m 是素数 这个条件。如果 m 不是素数，用费马小定理求逆元是错误的！

竞赛中如何选择？ 如果题目明确说明模数 p 是一个大素数 (常见的如 10^9 + 7, 998244353)，并且你需要求逆元，使用费马小定理 + 快速幂通常更方便。 如果模数 m 不确定是否为素数，或者可能不是素数，必须使用 ExGCD 来求解逆元。

方法三：线性递推求逆元 (Linear Method for Inverses) 有时我们需要预处理出 1 到 n 所有数模 p (p 通常是素数) 的逆元。如果对每个数都用 ExGCD 或快速幂，复杂度是 O(n log p)。存在一种 O(n) 的线性递推方法。

令 inv[i] 表示 i 模 p 的逆元。我们要求 inv[i]。 设 p = k * i + r，其中 $0 \le r < i$ (这是带余除法的定义)。 则 k = floor(p / i)，r = p % i。 将上式置于模 p 意义下：

$$k \times i + r \equiv 0 \pmod p$$ $$k \times i \equiv -r \pmod p$$

两边同乘以 $i^{-1} \times r^{-1}$ (假设它们都存在)：

$$k \times r^{-1} \equiv -i^{-1} \pmod p$$ $$i^{-1} \equiv -k \times r^{-1} \pmod p$$

将 k 和 r 的表达式代回：

$$\text{inv}[i] \equiv -\lfloor p/i \rfloor \times \text{inv}[p \pmod i] \pmod p $$

由于 p % i < i，inv[p % i] 在计算 inv[i] 之前已经被计算出来了。这就是递推关系。 边界是 inv[1] = 1。

#include <bits/stdc++.h>

using namespace std;
using ll = long long; // 使用 ll 代替 long long

const int MAXN = 200005; // 需要求逆元的上限
ll inv[MAXN];
ll p = 1e9 + 7; // 假设模数为素数

// 线性预处理 1 到 n 的逆元
// 时间复杂度: O(n)
void precompute_inverses(int n) {
    inv[0] = 0; // 0 没有逆元，或者根据题目定义
    inv[1] = 1;
    for (int i = 2; i <= n; ++i) {
        // inv[i] = (p - (p / i) * inv[p % i] % p) % p;
        // 更安全的写法，防止 p/i * inv[p%i] 本身是 p 的倍数导致结果为 0
        ll k = p / i;
        ll r = p % i;
        inv[i] = (p - k) * inv[r] % p;
        // 确保结果是非负的 (虽然按理说 k < p，inv[r] > 0，结果模 p 不会负)
        // inv[i] = (inv[i] % p + p) % p; // 如果不确定中间结果，加一层保险
    }
}

int main() {
    ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
    int n = 10;
    if (n >= MAXN) { /* Handle error or resize */ }
    precompute_inverses(n);
    cout << "Inverses mod " << p << " up to " << n << ":\n";
    for (int i = 1; i <= n; ++i) {
        cout << "inv(" << i << ") = " << inv[i] << "\n";
        // 验证: (i * inv[i]) % p 应该等于 1
        if (i > 0) {
             cout << "  Check: (" << i << " * " << inv[i] << ") % " << p << " = " << (1ll * i * inv[i]) % p << "\n";
        }
    }
    return 0;
}

注意: 线性求逆元方法也要求模数 p 是素数，因为它隐含地假设了 p % i (即 r) 的逆元存在 (r 不为 0 且与 p 互素)。当 p 是素数且 0 < r < p 时，此条件自然满足。

三、筛法：高效寻找素数

在很多问题中，我们需要快速得到一定范围内的所有素数，或者计算某些数论函数的值。筛法 (Sieve) 就是为此而生的。

1. 埃拉托斯特尼筛法 (Sieve of Eratosthenes)

最古老、最著名的筛法。

思想: 从 2 开始，将每个素数的倍数都标记为合数。

1. 创建一个布尔数组 is_composite[]，初始化所有值为 false。is_composite[0] 和 is_composite[1] 设为 true。
2. 从 p = 2 开始，如果 is_composite[p] 为 false (p 是素数)，则将 p 的所有倍数 2p, 3p, 4p, ... (直到上限 n) 标记为 true (合数)。
3. 继续找下一个 is_composite 仍为 false 的数 p，重复步骤 2。
4. 直到 p * p > n 时停止。

#include <bits/stdc++.h>

using namespace std;
using ll = long long; // 使用 ll 代替 long long

const int MAXN = 1000005;
bool is_composite[MAXN]; // is_composite[i] = true 表示 i 是合数
vector<int> primes;      // 存储找到的素数

// 埃氏筛法，找出 [2, n] 范围内的所有素数
// 时间复杂度: O(n log log n)
void sieve_eratosthenes(int n) {
    if (n >= MAXN) { /* Handle error or resize */ return; }
    fill(is_composite, is_composite + n + 1, false);
    is_composite[0] = is_composite[1] = true;
    // 注意 p*p 溢出 int 的问题
    for (int p = 2; (ll)p * p <= n; ++p) { // 强制转 ll 判断
        if (!is_composite[p]) { // p 是素数
            // 将 p 的倍数标记为合数
            // 从 p*p 开始标记，因为小于 p*p 的合数已被更小的素数标记过
            // 注意 i 溢出 int 的问题
            for (ll i = (ll)p * p; i <= n; i += p) { // 使用 ll 迭代
                 if (i < MAXN) { // 检查数组越界
                    is_composite[i] = true;
                 } else {
                    // 如果 i 超过了数组大小，理论上不应发生，因为 i<=n 且 n < MAXN
                    break;
                 }
            }
        }
    }
    // 收集素数
    primes.clear();
    for (int p = 2; p <= n; ++p) {
        if (!is_composite[p]) {
            primes.push_back(p);
        }
    }
}

int main() {
    ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
    int n = 100;
    sieve_eratosthenes(n);
    cout << "Primes up to " << n << ":\n";
    for (int p : primes) {
        cout << p << " ";
    }
    cout << "\n";
    return 0;
}

复杂度分析: 埃氏筛法的复杂度是 $O(N \log \log N)$，非常接近线性。直观理解，每个数 k 被其素因子标记。对于素数 p，它会标记 N/p 个数。所有素数加起来，$\sum_{p \le N} \frac{N}{p} \approx N \log \log N$。

2. 线性筛法 (Linear Sieve / Euler Sieve)

埃氏筛法的一个缺点是，一个合数可能会被它的多个素因子重复标记（例如 12 会被 2 和 3 都标记）。线性筛法通过巧妙的设计，确保每个合数只被其最小的素因子筛掉一次，从而达到严格的 O(N) 复杂度。

思想: 维护一个当前找到的素数列表 primes[]。 遍历 i 从 2 到 n：

1. 如果 i 没有被标记过（即 is_composite[i] 为 false），说明 i 是素数，将其加入 primes 列表。
2. 遍历已找到的素数 p (从 primes[0] 开始):
   • 将 i * p 标记为合数 (is_composite[i * p] = true)。
   • 关键一步: 如果 i 是 p 的倍数 (i % p == 0)，则停止用当前 i 去筛后续的素数。
     • 原因: 如果 i % p == 0，说明 p 是 i 的最小素因子（因为 p 是从 primes 列表里从小到大取的）。那么对于任何比 p 更大的素数 p' (即 primes 列表中 p 后面的素数)，i * p' 这个合数的最小素因子仍然是 p (因为 i 包含因子 p)，而不是 p'。根据线性筛的原则（每个合数只被其最小素因子筛掉），i * p' 应该在后面遍历到 (i * p' / p) 这个数时，被 p 筛掉，而不是现在被 p' 筛掉。因此，当 i % p == 0 时，我们必须 break，停止用当前的 i 和更大的素数 p' 去筛，以保证每个合数只被其最小素因子筛一次。

#include <bits/stdc++.h>

using namespace std;
using ll = long long; // 使用 ll 代替 long long

const int MAXN = 1000005;
bool is_composite[MAXN];
vector<int> primes;
int min_prime_factor[MAXN]; // 可选：记录每个数的最小素因子

// 线性筛法 (欧拉筛)
// 时间复杂度: O(n)
void sieve_linear(int n) {
    if (n >= MAXN) { /* Handle error or resize */ return; }
    fill(is_composite, is_composite + n + 1, false);
    // min_prime_factor[0] = min_prime_factor[1] = 1; // or 0, based on def
    is_composite[0] = is_composite[1] = true;
    primes.clear();

    for (int i = 2; i <= n; ++i) {
        if (!is_composite[i]) {
            primes.push_back(i);
            min_prime_factor[i] = i; // i 的最小素因子是它自己
        }
        // 尝试用已找到的素数 primes[j] (即 pj) 去筛掉 i * pj
        for (int pj : primes) {
            // 如果 i * pj 超过范围 n，停止
            // 注意溢出
            if (1ll * i * pj > n) { // 使用 1ll 保证计算不溢出
                break;
            }
            // i * pj 肯定在 MAXN 范围内，因为 i <= n < MAXN, pj <= n < MAXN
            // 但如果 MAXN 很大，i*pj 仍然可能超过 MAXN，需要检查
            // if (i * pj >= MAXN) break; // 更严谨的检查，如果MAXN不是n+1的话

            is_composite[i * pj] = true;
            min_prime_factor[i * pj] = pj; // pj 是 i*pj 的一个素因子，且是最小的那个

            // 关键: 保证每个合数只被其最小素因子筛一次
            if (i % pj == 0) {
                // 如果 i 能被 pj 整除, 说明 pj 是 i 的最小素因子
                // 那么 i * primes[k] (k>j) 的最小素因子也是 pj
                // 这些数应该在后面 i' = i * primes[k] / pj 时被 pj 筛掉
                // 所以这里 break, 不再用 i 去筛更大的素数
                break;
            }
        }
    }
}

int main() {
    ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
    int n = 100;
    sieve_linear(n);
    cout << "Primes up to " << n << " using linear sieve:\n";
    for (int p : primes) {
        cout << p << " ";
    }
    cout << "\n";

    // Optional: Print minimum prime factors
    // cout << "\nMinimum prime factor for numbers up to " << n << ":\n";
    // for (int i = 2; i <= n; ++i) {
    //     cout << "min_prime_factor(" << i << ") = " << min_prime_factor[i] << "\n";
    // }

    return 0;
}

线性筛的优势:

• 严格 O(N) 复杂度，在 N 达到 $10^7$ 或 $10^8$ 级别时比埃氏筛有明显优势。
• 可以在筛素数的同时，线性地预处理出一些重要的积性函数 (Multiplicative Function) 的值，例如欧拉函数 phi、莫比乌斯函数 mu、约数个数 d、约数和 sigma 等。这是线性筛在竞赛中应用更广泛的原因。

什么是积性函数？ 一个数论函数 f(n)，如果对于任意一对互素的正整数 a 和 b (即 gcd(a, b) = 1)，都有 f(a * b) = f(a) * f(b)，则称 f 为积性函数。 如果对于任意正整数 a 和 b，都有 f(a * b) = f(a) * f(b)，则称 f 为完全积性函数。

线性筛利用 n = i * p (其中 p 是 n 的最小素因子) 这个关系，结合积性函数的性质，可以在 O(1) 的时间内由 f(i) 和 f(p)（或相关值）计算出 f(n)。

四、数论函数与定理

除了基础运算，一些重要的数论函数和定理在解题中扮演着关键角色。

1. 欧拉函数 (Euler's Totient Function)

定义: 对于正整数 n，欧拉函数 $\phi(n)$ (读作 phi of n) 定义为小于或等于 n 的正整数中，与 n 互素 (即 gcd(k, n) = 1) 的数的个数。 例如：$\phi(1)=1$, $\phi(6)=2$ (与 6 互素的有 1, 5)，$\phi(7)=6$ (1, 2, 3, 4, 5, 6 都与素数 7 互素)。

性质：

1. 若 p 是素数，则 $\phi(p) = p - 1$。 (因为 1 到 p-1 都与 p 互素)
2. 若 p 是素数，k >= 1，则 $\phi(p^k) = p^k - p^{k-1} = p^k(1 - 1/p)$。 (在 $1$ 到 $p^k$ 中，不与 $p^k$ 互素的数是 $p$ 的倍数，共有 $p^{k-1}$ 个，即 $p, 2p, ..., p^k$)
3. 欧拉函数是积性函数。 即若 gcd(a, b) = 1，则 $\phi(a \times b) = \phi(a) \times \phi(b)$。

计算公式: 基于算术基本定理和积性性质，如果 n 的标准分解式为 $n = p_1^{a_1} p_2^{a_2} \cdots p_k^{a_k}$，则：

$$\phi(n) = \phi(p_1^{a_1}) \phi(p_2^{a_2}) \cdots \phi(p_k^{a_k}) $$$$= (p_1^{a_1} - p_1^{a_1-1}) \cdots (p_k^{a_k} - p_k^{a_k-1}) $$$$= p_1^{a_1}(1 - \frac{1}{p_1}) \cdots p_k^{a_k}(1 - \frac{1}{p_k}) $$$$= (p_1^{a_1} \cdots p_k^{a_k}) \times (1 - \frac{1}{p_1}) \cdots (1 - \frac{1}{p_k}) $$$$\phi(n) = n \prod_{i=1}^{k} (1 - \frac{1}{p_i})$$

其中 $p_1, ..., p_k$ 是 n 的所有不同素因子。

如何计算 $\phi(n)$？

• 单个计算： 先对 n 进行质因数分解，然后套用上面的公式。复杂度瓶颈在质因数分解，即 $O(\sqrt{n})$。

#include <bits/stdc++.h>

using namespace std;
using ll = long long; // 使用 ll 代替 long long

// 计算单个 n 的欧拉函数 phi(n)
// 时间复杂度: O(sqrt(n))
ll phi(ll n) {
    if (n <= 0) return 0; // 处理非正整数输入
    ll res = n;
    // 使用 ll 迭代
    for (ll i = 2; i * i <= n; ++i) {
        if (n % i == 0) {
            res = res / i * (i - 1); // 等价于 res = res * (1 - 1/i)
            while (n % i == 0) {
                n /= i;
            }
        }
    }
    if (n > 1) { // 如果 n 还有一个大于 sqrt(n) 的质因子
        res = res / n * (n - 1);
    }
    return res;
}

int main() {
    ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
    ll n;
    cin >> n;
    cout << "phi(" << n << ") = " << phi(n) << "\n";
    return 0;
}

• 线性筛预处理： 利用线性筛，可以在 O(N) 时间内预处理出 1 到 N 所有数的欧拉函数值。

  在线性筛的过程中，当筛到合数 n = i * p (其中 p 是 n 的最小素因子 primes[j]) 时：

  1. 如果 i % p == 0 (即 p | i)： 这意味着 n 和 i 含有相同的素因子 p，只是 p 的指数在 n 中比在 i 中多 1。根据 $\phi(p^k) = p^k - p^{k-1} = p \times (p^{k-1} - p^{k-2}) = p \times \phi(p^{k-1})$，以及积性函数的性质，可以推导出 $\phi(n) = \phi(i \times p) = \phi(i) \times p$。
  2. 如果 i % p != 0 (即 p 不是 i 的因子)： 这意味着 i 和 p 互素。根据 $\phi$ 是积性函数，$\phi(n) = \phi(i \times p) = \phi(i) \times \phi(p)$。因为 p 是素数，$\phi(p) = p - 1$。所以 $\phi(n) = \phi(i) \times (p - 1)$。

#include <bits/stdc++.h>

using namespace std;
using ll = long long; // 使用 ll 代替 long long

const int MAXN = 1000005;
bool is_composite[MAXN];
vector<int> primes;
ll phi_val[MAXN]; // 存储欧拉函数值，用 ll

// 线性筛法预处理欧拉函数 phi
// 时间复杂度: O(n)
void sieve_phi(int n) {
    if (n >= MAXN) { /* Handle error or resize */ return; }
    fill(is_composite, is_composite + n + 1, false);
    phi_val[1] = 1;
    is_composite[0] = is_composite[1] = true;
    primes.clear();

    for (int i = 2; i <= n; ++i) {
        if (!is_composite[i]) {
            primes.push_back(i);
            phi_val[i] = i - 1; // p 是素数, phi(p) = p - 1
        }
        for (int pj : primes) {
            if (1ll * i * pj > n) break; // 注意溢出
            is_composite[i * pj] = true;
            if (i % pj == 0) {
                // pj 是 i 的一个素因子 (也是 i*pj 的最小素因子)
                // phi(i*pj) = phi(i) * pj
                phi_val[i * pj] = phi_val[i] * pj; // pj 是 int, phi_val[i] 是 ll
                break; // 保证线性时间
            } else {
                // i 和 pj 互素
                // phi(i*pj) = phi(i) * phi(pj) = phi(i) * (pj - 1)
                phi_val[i * pj] = phi_val[i] * (pj - 1);
            }
        }
    }
}

int main() {
    ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
    int n = 20;
    sieve_phi(n);
    cout << "Euler's totient function phi(k) for k from 1 to " << n << ":\n";
    for (int i = 1; i <= n; ++i) {
        cout << "phi(" << i << ") = " << phi_val[i] << "\n";
    }
    return 0;
}

2. 欧拉定理 (Euler's Theorem)

这是费马小定理的推广。

定理内容: 若 a 和 m 是互素的正整数 (即 gcd(a, m) = 1)，则：

$$a^{\phi(m)} \equiv 1 \pmod m$$

证明思路: 考虑所有小于 m 且与 m 互素的 $\phi(m)$ 个数 $r_1, r_2, ..., r_{\phi(m)}$。因为 gcd(a, m) = 1，所以集合 $\{ar_1 \pmod m, ar_2 \pmod m, ..., ar_{\phi(m)} \pmod m\}$ 恰好是集合 $\{r_1, r_2, ..., r_{\phi(m)}\}$ 的一个排列。将两组数分别相乘并模 m，得到：

$$(ar_1)(ar_2)\cdots(ar_{\phi(m)}) \equiv r_1 r_2 \cdots r_{\phi(m)} \pmod m $$$$a^{\phi(m)} (r_1 r_2 \cdots r_{\phi(m)}) \equiv (r_1 r_2 \cdots r_{\phi(m)}) \pmod m $$

令 $R = r_1 r_2 \cdots r_{\phi(m)}$。因为所有 $r_i$ 都与 m 互素，所以它们的乘积 R 也与 m 互素，即 gcd(R, m) = 1。根据同余的性质，我们可以约掉 R：

$$a^{\phi(m)} \equiv 1 \pmod m$$

应用:

1. 降幂: 在计算 $a^b \pmod m$ 时，如果 gcd(a, m) = 1，我们可以将指数 b 对 $\phi(m)$ 取模，即：

   $$a^b \equiv a^{b \pmod{\phi(m)}} \pmod m \quad (\text{如果 } b \ge \phi(m)) $$

   注意： 这个降幂公式严格来说要求 $b \ge \phi(m)$。对于 $b < \phi(m)$，直接计算 $a^b \pmod m$。 扩展欧拉定理 (Extended Euler's Theorem): 处理 gcd(a, m) != 1 或 b < phi(m) 的情况。

   $$a^b \equiv \begin{cases} a^b \pmod m & \text{if } b < \phi(m) \\ a^{b \pmod{\phi(m)} + \phi(m)} \pmod m & \text{if } b \ge \phi(m) \end{cases} $$

   这个形式可以统一处理所有情况（即使 gcd(a, m) != 1）。在竞赛中，如果你需要计算 $a^b \pmod m$，可以直接用扩展欧拉定理的公式来降幂：先计算 $b' = b \pmod{\phi(m)}$，然后令新的指数为 $exp = b' + \phi(m)$。但是，需要判断原始的 $b$ 是否小于 $\phi(m)$，如果 $b < \phi(m)$，则指数应该直接是 $b$。

   #include <bits/stdc++.h>
   using namespace std;
   using ll = long long; // 使用 ll 代替 long long
   
   // 快速幂 (复用)
   ll qpow(ll base, ll exp, ll mod) {
       ll res = 1;
       base %= mod; if(base < 0) base += mod;
       while (exp > 0) {
           if (exp % 2 == 1) res = (res * base) % mod;
           base = (base * base) % mod;
           exp /= 2;
       }
       return res;
   }
   
   // 计算单个 phi(n) (复用)
   ll phi(ll n) {
       if (n <= 0) return 0;
       ll res = n;
       for (ll i = 2; i * i <= n; ++i) {
           if (n % i == 0) {
               res = res / i * (i - 1);
               while (n % i == 0) n /= i;
           }
       }
       if (n > 1) res = res / n * (n - 1);
       return res;
   }
   
   // 使用扩展欧拉定理计算 a^b mod m
   // 注意: b 可能是非常大的数
   ll power_extended_euler(ll a, ll b, ll m) {
       if (m == 1) return 0; // 对 1 取模结果总是 0
       if (b == 0) return 1 % m; // a^0 = 1
   
       ll phm = phi(m);
       ll exp;
       bool b_ge_phm = (b >= phm); // 判断 b 是否大于等于 phi(m)
   
       if (b_ge_phm) {
            exp = b % phm + phm; // 核心：指数变为 b % phi(m) + phi(m)
       } else {
           exp = b; // b < phi(m)，指数不变
       }
   
       // 使用快速幂计算 a^exp mod m
       return qpow(a, exp, m);
   }
   
   // 处理大指数 b (以字符串形式给出)
   ll power_large_b(ll a, string b_str, ll m) {
        if (m == 1) return 0;
        // 处理 b=0 的情况
        if (b_str == "0") return 1 % m;
   
        ll phm = phi(m);
        ll b_mod_phm = 0;
        bool b_ge_phm = false; // 判断 b 是否 >= phm
   
        // 计算 b mod phi(m)，同时判断 b 是否 >= phi(m)
        for (char digit : b_str) {
            b_mod_phm = (b_mod_phm * 10 + (digit - '0'));
            if (b_mod_phm >= phm) { // 只要中间结果超过 phi(m) 就说明 b >= phi(m)
                b_ge_phm = true;
                // 取模要在比较之后，否则无法正确判断 b >= phm
                b_mod_phm %= phm;
            }
        }
        // 如果 b_ge_phm 为 true，b_mod_phm 存储的是 b % phm
        // 如果 b_ge_phm 为 false，b_mod_phm 存储的是 b 的实际值 (因为 b < phm)
   
        ll exp;
        if (b_ge_phm) {
            exp = b_mod_phm + phm;
        } else {
            // 如果 b_ge_phm 是 false，说明 b < phm
            exp = b_mod_phm; // 指数就是 b
        }
        return qpow(a, exp, m);
   }
   

   int main() { ios::sync_with_stdio(0); cin.tie(0); cout.tie(0); ll a = 3; ll b = 100; ll m = 10; // gcd(3, 10) = 1, phi(10) = 4 cout << a << "^" << b << " mod " << m << " = " << power_extended_euler(a, b, m) << endl; // b=100>=4. exp = 100%4 + 4 = 0+4=4. 3^4 = 81. 81 mod 10 = 1.

    a = 2;
    b = 10;
    m = 6; // gcd(2, 6) = 2 != 1. phi(6)=2
    cout << a << "^" << b << " mod " << m << " = " << power_extended_euler(a, b, m) << endl; // b=10 >= phi(6)=2. exp = 10%2 + 2 = 0+2 = 2. 2^2 = 4. 4 mod 6 = 4.
   
    a = 7;
    string b_large = "123456789123456789";
    m = 100; // phi(100) = 40
    cout << a << "^(" << b_large << ") mod " << m << " = " << power_large_b(a, b_large, m) << endl;
    // b mod 40: 123...789 mod 40 = 9 mod 40. b >= 40.
    // exp = 9 + 40 = 49.
    // Need 7^49 mod 100.
    // 7^1=7, 7^2=49, 7^3=343=43, 7^4=43*7=301=1 (mod 100)
    // 7^49 = 7^(4*12 + 1) = (7^4)^12 * 7^1 = 1^12 * 7 = 7 (mod 100).
    // Let's check the code output.
   
    return 0;
   

   }

   **注意:** 计算 $b \pmod{\phi(m)}$ 时，如果 `b` 本身就非常大 (比如指数塔 $a^{b^c}$ 或者用字符串表示)，你需要先处理大数 `b` 对 $\phi(m)$ 的取模，并正确判断 $b$ 是否大于等于 $\phi(m)$。
   
   

2. 求解模逆元: 如果 m 不是素数，但 gcd(a, m) = 1，可以用欧拉定理求逆元： $a^{\phi(m)-1} \pmod m$ 就是 a 模 m 的逆元。这需要先计算 $\phi(m)$，然后做一次快速幂。复杂度是 $O(\sqrt{m} + \log(\phi(m)))$ 或 O(筛法复杂度 + $\log(\phi(m)))$。通常 ExGCD 更直接。

3. 费马小定理 (Fermat's Little Theorem)

前面已经提过，它是欧拉定理在模数 p 是素数时的特殊情况。 定理内容: 若 p 是素数，a 不是 p 的倍数，则 $a^{p-1} \equiv 1 \pmod p$。 另一种形式：对于任意整数 a 和素数 p，有 $a^p \equiv a \pmod p$。

主要应用：

1. 素性测试 (Primality Testing): Miller-Rabin 测试的基础。如果 $a^{p-1} \not\equiv 1 \pmod p$ 对某个 a 成立，则 p 一定是合数。但反之不一定，存在伪素数 (Carmichael numbers)。
2. 求模逆元 (模数为素数时): $a^{-1} \equiv a^{p-2} \pmod p$。

4. 中国剩余定理 (Chinese Remainder Theorem, CRT)

用于求解一元线性同余方程组 (System of Linear Congruences)。

标准形式: 求解 x 满足以下方程组：

$$\begin{cases} x \equiv a_1 \pmod{m_1} \\ x \equiv a_2 \pmod{m_2} \\ \vdots \\ x \equiv a_k \pmod{m_k} \end{cases} $$

条件： 模数 $m_1, m_2, ..., m_k$ 两两互素 (pairwise coprime)。

解法： 令 $M = m_1 \times m_2 \times \cdots \times m_k$。 对每个 $i$（从 1 到 k），令 $M_i = M / m_i$。 因为 $m_i$ 与 $M_i$ 互素 ($m_i$ 的因子只在 $m_i$ 里，$M_i$ 的因子是除了 $m_i$ 之外的所有 $m_j$)，所以 $M_i$ 模 $m_i$ 的逆元存在。令 $t_i$ 是 $M_i$ 模 $m_i$ 的逆元，即 $M_i t_i \equiv 1 \pmod{m_i}$。可以用 ExGCD 求解 $t_i$。

则方程组的一个特解为：

$$x_0 = \sum_{i=1}^{k} a_i M_i t_i$$

验证: 对于任意第 j 个方程 $x \equiv a_j \pmod{m_j}$： 当 $i = j$ 时，$a_j M_j t_j \equiv a_j \times 1 \equiv a_j \pmod{m_j}$。 当 $i \neq j$ 时，$M_i = M / m_i$ 包含因子 $m_j$ (因为模数两两互素)，所以 $M_i \equiv 0 \pmod{m_j}$。因此 $a_i M_i t_i \equiv 0 \pmod{m_j}$。 所以 $x_0 = \sum_{i=1}^{k} a_i M_i t_i \equiv 0 + \dots + a_j + \dots + 0 \equiv a_j \pmod{m_j}$。 $x_0$ 满足所有方程。

方程组的通解为：

$$x \equiv x_0 \pmod M$$

即所有解的形式为 $x = x_0 + kM$，其中 k 是任意整数。 最小非负整数解为 $x_0 \pmod M$。

#include <bits/stdc++.h>

using namespace std;
using ll = long long; // 使用 ll 代替 long long

// ExGCD (复用)
ll exgcd(ll a, ll b, ll &x, ll &y) {
    if (b == 0) { x = 1; y = 0; return a; }
    ll d = exgcd(b, a % b, x, y);
    ll temp = x; x = y; y = temp - (a / b) * y;
    return d;
}

// 模逆元 (复用 ExGCD)
ll modInverse(ll a, ll m) {
    ll x, y;
    ll d = exgcd(a, m, x, y);
    if (d != 1) return -1; // 逆元不存在
    return (x % m + m) % m;
}

// 安全乘法取模，防止 a*b 溢出 ll
ll safe_mul(ll a, ll b, ll m) {
    a %= m; b %= m;
    ll res = 0;
    while (b > 0) {
        if (b & 1) res = (res + a) % m;
        a = (a + a) % m;
        b >>= 1;
    }
    return res;
    // 或者如果 __int128 可用: return (ll)((__int128)a * b % m);
}


// 中国剩余定理 (CRT)
// 求解 x = a_i (mod m_i), 其中 m_i 两两互素
// a: 存储余数的数组 (vector<ll>)
// m: 存储模数的数组 (vector<ll>)
// 返回最小非负整数解 x 和总模数 M
// 时间复杂度: O(k * log(M)), k 是方程个数, M 是模数乘积
pair<ll, ll> crt(const vector<ll>& a, const vector<ll>& m) {
    int k = a.size();
    if (k == 0) return {0, 1};

    ll M = 1; // 所有模数的乘积
    for (ll mi : m) {
        // 需要检查 M * mi 是否溢出 ll
        if ((double)M * mi > 2e18) { // 粗略检查溢出，更精确需要 __int128
             // Handle overflow error, maybe return {-2, -2}
             // 在竞赛中，有时 M 会爆 ll，题目可能要求对另一个大数取模，或者保证 M 不爆
             // 这里假设 M 不会爆 ll
        }
        M *= mi;
    }
    if (M == 0) { // 如果某个 mi 是 0，处理一下，这里简单返回错误
        return {-1, -1};
    }


    ll x0 = 0;
    for (int i = 0; i < k; ++i) {
        if (m[i] == 0) return {-1, -1}; // 模数不能为 0
        ll Mi = M / m[i];
        ll ti = modInverse(Mi, m[i]);
        if (ti == -1) {
            // m_i 不互素时会发生，但 CRT 假设互素
            return {-1, -1}; // 表示无解或输入错误
        }

        // 计算项 (a[i] * Mi * ti) % M
        // 需要使用安全乘法防止溢出
        ll term = safe_mul(a[i], Mi, M);
        term = safe_mul(term, ti, M);
        // 如果用 __int128:
        // __int128 term128 = (__int128)a[i] * Mi % M * ti % M;
        // x0 = (x0 + (ll)term128) % M;

        x0 = (x0 + term) % M;
    }

    // 结果可能是负数，调整为最小非负解
    return {(x0 + M) % M, M}; // 返回解 x0 和总模数 M
}

int main() {
    ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);

    // 例子:
    // x = 2 (mod 3)
    // x = 3 (mod 5)
    // x = 2 (mod 7)
    vector<ll> a = {2, 3, 2};
    vector<ll> m = {3, 5, 7}; // 3, 5, 7 两两互素

    pair<ll, ll> result = crt(a, m);

    if (result.first != -1) {
        cout << "System:\n";
        for (size_t i = 0; i < a.size(); ++i) {
            cout << "x = " << a[i] << " (mod " << m[i] << ")\n";
        }
        cout << "Smallest non-negative solution x = " << result.first << "\n"; // x = 23
        cout << "General solution x = " << result.first << " + k * " << result.second << "\n"; // M = 105
    } else {
        cout << "Error in CRT calculation (e.g., moduli not pairwise coprime or overflow).\n";
    }

    return 0;
}

扩展中国剩余定理 (Extended CRT / Excrt) 用于解决模数 $m_i$ 不一定两两互素 的情况。

思想： 增量法，一次合并两个同余方程。 假设我们已经求得前 i-1 个方程的解，满足 $x \equiv ans \pmod M$，其中 $M = \text{lcm}(m_1, ..., m_{i-1})$。 现在要加入第 i 个方程 $x \equiv a_i \pmod{m_i}$。 这意味着当前的解 x 必须满足 $x = ans + kM$ (对某个整数 k)，并且 $x \equiv a_i \pmod{m_i}$。 代入得到：

$$ans + kM \equiv a_i \pmod{m_i}$$ $$kM \equiv a_i - ans \pmod{m_i}$$

这是一个关于 k 的模线性方程 (M)k \equiv (a_i - ans) \pmod{m_i}。 令 $A = M$, $B = m_i$, $C = (a_i - ans) \pmod{m_i}$ (确保 C 在 $[0, B-1]$ 范围内)。方程变为 $Ak \equiv C \pmod B$。 这个方程有解当且仅当 $d = \gcd(A, B)$ 能整除 $C$。 如果 $C \% d \neq 0$，则整个方程组无解。 如果有解，用 ExGCD 求解 $Ax' + By' = d$。得到一组解 $(x', y')$。 那么 $A(x' \times C/d) + B(y' \times C/d) = C$。 所以 $A(x' \times C/d) \equiv C \pmod B$。 因此， $k_0 = x' \times (C/d)$ 是方程 $Ak \equiv C \pmod B$ 的一个特解。 模 $B/d$ 的通解为 $k \equiv k_0 \pmod{B/d}$。 即 $k = k_0 + t (B/d)$，其中 t 是任意整数。

我们只需要一个特解 k，比如取最小非负解 $k = (k_0 \pmod{B/d} + B/d) \pmod{B/d}$。 将这个 k 代回 $x = ans + k M$，得到满足前 i 个方程的一个新解 $ans_{new} = ans + k M$。 新的总模数变为 $M_{new} = \text{lcm}(M, m_i) = (M / \gcd(M, m_i)) \times m_i$ (注意计算 LCM 防溢出)。 重复这个过程直到合并完所有方程。

#include <bits/stdc++.h>

using namespace std;
using ll = long long; // 使用 ll 代替 long long

// ExGCD (复用)
ll exgcd(ll a, ll b, ll &x, ll &y) {
    if (b == 0) { x = 1; y = 0; return a; }
    ll d = exgcd(b, a % b, x, y);
    ll temp = x; x = y; y = temp - (a / b) * y;
    return d;
}

// GCD (复用)
ll gcd(ll a, ll b) { return b == 0 ? abs(a) : gcd(b, a % b); } // 处理负数，返回非负

// 最小公倍数 (复用 gcd)
ll lcm(ll a, ll b) {
    if (a == 0 || b == 0) return 0;
    ll common_divisor = gcd(a, b);
    // 防止溢出: (a / gcd) * b
    ll res_num = a / common_divisor;
    // 检查 res_num * b 是否溢出
    if (abs(b) > 0 && abs(res_num) > numeric_limits<ll>::max() / abs(b)) {
        return -2; // 返回特殊值表示溢出
    }
    return res_num * b;
}

// 安全乘法取模 (复用)
ll safe_mul(ll a, ll b, ll m) {
    a %= m; b %= m;
    if (a < 0) a += m;
    if (b < 0) b += m;
    ll res = 0;
    while (b > 0) {
        if (b & 1) res = (res + a);
        if (res >= m) res -=m; // 手动取模防止 res+a 溢出一点点
        a = (a + a);
        if (a >= m) a -= m;
        b >>= 1;
    }
    return res;
    // Or use __int128 if available and needed
    // return (ll)((__int128)a * b % m);
}


// 扩展中国剩余定理 (Excrt)
// 求解 x = a_i (mod m_i), m_i 不一定两两互素
// 返回 pair<ll, ll>: {最小非负解 x, 最终模数 M = lcm(m_i)}
// 如果无解，返回 {-1, -1}
// 如果 LCM 计算溢出，返回 {-2, -2}
pair<ll, ll> excrt(const vector<ll>& a, const vector<ll>& m) {
    int k = a.size();
    if (k == 0) return {0, 1};

    ll ans = a[0]; // 初始解满足第一个方程
    ll M = m[0];   // 初始模数
    if (M <= 0) return {-1, -1}; // 模数必须为正

    ans = (ans % M + M) % M; // 保证初始解非负

    for (int i = 1; i < k; ++i) {
        ll mi = m[i];
        ll ai = a[i];
        if (mi <= 0) return {-1, -1}; // 模数必须为正

        // 合并方程 x = ans (mod M) 和 x = ai (mod mi)
        // 求解 kM = ai - ans (mod mi)
        ll A = M;
        ll B = mi;
        ll C = (ai - ans % B + B) % B; // C = (ai - ans) mod B, 保证 C 非负

        ll x, y;
        ll d = exgcd(A, B, x, y); // d = gcd(M, mi)

        if (C % d != 0) {
            return {-1, -1}; // 方程组无解
        }

        // 求解 k = x * (C / d) (mod B / d)
        ll mod_k = B / d;
        // x * (C / d) 可能溢出，需要安全乘法
        // ll k0 = safe_mul(x, C / d, mod_k);
        // 使用 __int128 通常更方便且不易出错（如果允许）
        ll k0 = (ll)((__int128)x * (C / d) % mod_k);


        k0 = (k0 + mod_k) % mod_k; // 确保 k0 非负

        // 更新解和模数
        // ans = ans + k0 * M; // 可能溢出
        // M_new = lcm(M, mi); // 可能溢出

        ll M_new = lcm(M, mi);
        if (M_new == -2) return {-2, -2}; // LCM 溢出

        // 更新 ans: ans = ans + k0 * M
        // __int128 next_ans = (__int128)ans + (__int128)k0 * M;
        // ans = (ll)(next_ans % M_new);
        // 不用 __int128 的话，计算 k0 * M 可能会溢出
        // ans = (ans + safe_mul(k0, M, M_new)) % M_new; // 这种计算也不完全对

        // 标准做法：
        ans = ans + k0 * M; // 计算新的特解，此时 ans 可能很大
        M = M_new;          // 更新模数
        ans = (ans % M + M) % M; // 将解归约到 [0, M-1]
    }

    return {(ans + M) % M, M};
}

int main() {
    ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);

    // 例子: x = 1 (mod 4), x = 3 (mod 6)
    vector<ll> a1 = {1, 3};
    vector<ll> m1 = {4, 6};
    pair<ll, ll> result1 = excrt(a1, m1); // Expected: 9 mod 12
    cout << "System 1:\n";
    if (result1.first == -1) cout << "No solution.\n";
    else if (result1.first == -2) cout << "LCM Overflow.\n";
    else cout << "x = " << result1.first << " (mod " << result1.second << ")\n";

    // 例子: x = 2 (mod 3), x = 3 (mod 4), x = 1 (mod 5)
    vector<ll> a2 = {2, 3, 1};
    vector<ll> m2 = {3, 4, 5};
    pair<ll, ll> result2 = excrt(a2, m2); // Expected: 11 mod 60
    cout << "System 2:\n";
     if (result2.first == -1) cout << "No solution.\n";
    else if (result2.first == -2) cout << "LCM Overflow.\n";
    else cout << "x = " << result2.first << " (mod " << result2.second << ")\n";


    // 例子: 无解 x = 1 (mod 2), x = 0 (mod 4)
    vector<ll> a3 = {1, 0};
    vector<ll> m3 = {2, 4};
    pair<ll, ll> result3 = excrt(a3, m3); // Expected: -1
    cout << "System 3:\n";
    if (result3.first == -1) cout << "No solution.\n";
    else if (result3.first == -2) cout << "LCM Overflow.\n";
    else cout << "x = " << result3.first << " (mod " << result3.second << ")\n";


    return 0;
}

注意: Excrt 实现中，计算 k0 = x * (C / d)、ans = ans + k0 * M 以及 M = lcm(M, m[i]) 时要特别小心 ll 溢出。如果模数 $m_i$ 很大，导致 LCM 或中间乘积超过 $10^{18}$ 范围，必须使用 __int128 (如果编译器支持且题目允许) 或者实现更复杂的 ll 安全乘法和 LCM 计算。

五、组合计数与卢卡斯定理

数论与组合计数常常结合出现，尤其是在需要对组合数取模时。

1. 组合数取模

计算 $C(n, k) = \binom{n}{k} = \frac{n!}{k!(n-k)!} \pmod p$。

方法一：利用模逆元 (当 p 是素数) 如果模数 p 是素数，我们可以预处理阶乘 fact[i] = i! \pmod p$，以及阶乘的逆元 invfact[i] = (i!)^{-1} \pmod p$。$C(n, k) \equiv n! \times (k!)^{-1} \times ((n-k)! )^{-1} \pmod p$

$$\binom{n}{k} \equiv \text{fact}[n] \times \text{invfact}[k] \times \text{invfact}[n-k] \pmod p $$

预处理阶乘是 O(N)。预处理阶乘逆元：可以先用快速幂或 ExGCD 求出 fact[N] 的逆元 invfact[N]，然后利用 invfact[i-1] = invfact[i] * i % p 递推 O(N) 求出所有阶乘逆元。总复杂度 O(N + log p) 或 O(N)。

#include <bits/stdc++.h>

using namespace std;
using ll = long long; // 使用 ll 代替 long long

const int MAXN = 1000005;
ll p = 1e9 + 7; // 假设 p 是素数
ll fact[MAXN];
ll invfact[MAXN];

// 快速幂 (复用)
ll qpow(ll base, ll exp, ll mod) {
    ll res = 1; base %= mod; if(base < 0) base += mod;
    while (exp > 0) {
        if (exp % 2 == 1) res = (res * base) % mod;
        base = (base * base) % mod; exp /= 2;
    }
    return res;
}

// 预处理阶乘和阶乘逆元
// 时间复杂度: O(N + log p) (用快速幂求第一个逆元) 或 O(N) (如果 p 是素数)
void precompute_combinations(int n) {
    if (n >= MAXN) { /* Handle error or resize */ return; }
    fact[0] = 1;
    invfact[0] = 1;
    for (int i = 1; i <= n; ++i) {
        fact[i] = (fact[i - 1] * i) % p;
    }
    // 计算 fact[n] 的逆元
    // 如果 p 不是素数，这里必须用 ExGCD 求逆元
    // invfact[n] = modInverse_exgcd(fact[n], p);
    invfact[n] = qpow(fact[n], p - 2, p); // 使用费马小定理 (仅当 p 是素数)
    if (invfact[n] == -1 && fact[n] != 0) { /* Handle error: inverse doesn't exist */ return;}

    // 递推计算其他阶乘逆元
    for (int i = n - 1; i >= 1; --i) {
        // invfact[i] = invfact[i+1] * (i+1) mod p
        // 注意 (i+1) 可能大于 p，所以要取模
        invfact[i] = (invfact[i + 1] * ((i + 1) % p)) % p;
    }
}

// 计算 C(n, k) mod p
// 时间复杂度: O(1) (在预处理之后)
ll combinations(int n, int k) {
    if (k < 0 || k > n || n >= MAXN) { // 检查 k 范围和 n 是否超出预处理范围
        return 0;
    }
    // C(n, k) = n! / (k! * (n-k)!) mod p
    // C(n, k) = fact[n] * invfact[k] * invfact[n-k] mod p
    ll res = fact[n];
    res = (res * invfact[k]) % p;
    res = (res * invfact[n - k]) % p;
    return (res % p + p) % p; // 确保结果非负
}

int main() {
    ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
    int max_n = 1000; // 预处理到 1000
    precompute_combinations(max_n);

    int n = 10, k = 3;
    cout << "C(" << n << ", " << k << ") mod " << p << " = " << combinations(n, k) << endl; // C(10, 3) = 120

    n = 5, k = 5;
    cout << "C(" << n << ", " << k << ") mod " << p << " = " << combinations(n, k) << endl; // C(5, 5) = 1

    n = 1000, k = 500;
    cout << "C(" << n << ", " << k << ") mod " << p << " = " << combinations(n, k) << endl;

    return 0;
}

方法二：卢卡斯定理 (Lucas's Theorem) 当 n 和 k 很大，但模数 p 是素数且不大时，直接预处理阶乘不可行。卢卡斯定理提供了一种拆解计算的方法。

定理内容: 对于非负整数 n, k 和素数 p，令 n = n_0 + n_1 p + n_2 p^2 + ... 和 k = k_0 + k_1 p + k_2 p^2 + ... 分别是 n 和 k 的 p 进制表示 ($0 \le n_i, k_i < p$)。则：

$$\binom{n}{k} \equiv \prod_{i=0} \binom{n_i}{k_i} \pmod p $$

其中，如果 $k_i > n_i$，则 $\binom{n_i}{k_i} = 0$。

理解: 该定理将计算大组合数模小素数 p 的问题，分解为计算若干个小于 p 的组合数 $\binom{n_i}{k_i} \pmod p$ 的乘积。 计算 $\binom{n_i}{k_i} \pmod p$ 时，因为 $n_i, k_i < p$，可以使用方法一（预处理阶乘和逆元到 p-1）来快速计算。

实现: 我们可以递归或者迭代地实现 Lucas 定理。 核心是：

$$\text{Lucas}(n, k, p) = \binom{n \pmod p}{k \pmod p} \times \text{Lucas}(\lfloor n/p \rfloor, \lfloor k/p \rfloor, p) \pmod p $$

递归边界是当 k = 0 时，返回 1。

#include <bits/stdc++.h>

using namespace std;
using ll = long long; // 使用 ll 代替 long long

// --- 需要组合数计算 C(n, k) mod p (n, k < p) ---
const int MAXP = 1005; // 假设 p 不会很大，预处理阶乘到 p-1
ll fact_lucas[MAXP];
ll invfact_lucas[MAXP];
ll current_p = -1; // 当前 Lucas 定理使用的素数模数

// 快速幂 (复用)
ll qpow(ll base, ll exp, ll mod) {
    ll res = 1; base %= mod; if(base < 0) base += mod;
    while (exp > 0) {
        if (exp % 2 == 1) res = (res * base) % mod;
        base = (base * base) % mod; exp /= 2;
    }
    return res;
}
// 模逆元 (复用费马小定理)
ll modInverse_fermat(ll a, ll p) {
    if (p <= 1) return -1; a %= p; if(a<0) a+=p; if (a == 0) return -1;
    if (p - 2 < 0) return (p == 2 && a == 1) ? 1 : -1;
    return qpow(a, p - 2, p);
}

// 预处理 Lucas 定理所需的阶乘和逆元 (模 p)
void precompute_lucas(ll p) {
    if (p == current_p) return; // 如果模数没变，无需重新计算
    if (p >= MAXP) { /* Handle error: p is too large for precomputation array */ return; }
    current_p = p;
    fact_lucas[0] = 1;
    invfact_lucas[0] = 1;
    for (int i = 1; i < p; ++i) {
        fact_lucas[i] = (fact_lucas[i - 1] * i) % p;
    }
    invfact_lucas[p - 1] = modInverse_fermat(fact_lucas[p - 1], p);
    if (invfact_lucas[p - 1] == -1 && fact_lucas[p-1] != 0) { /* Error */ return; }
    for (int i = p - 2; i >= 1; --i) {
        invfact_lucas[i] = (invfact_lucas[i + 1] * ((i + 1) % p)) % p;
    }
}

// 计算 C(n, k) mod p (n, k < p)
ll combinations_small(ll n, ll k, ll p) {
    if (k < 0 || k > n) return 0;
    if (k == 0 || k == n) return 1;
    if (k > n / 2) k = n - k; // 优化 C(n, k) = C(n, n-k)
    if (n >= p || k >= p || p != current_p) {
        // 如果 n, k >= p 或模数变化，需要报错或重新预处理
         if (p != current_p) precompute_lucas(p); // 尝试重新预处理
         if (n >= p || k >= p || p!= current_p) return -1; // 预处理失败或参数错误
    }
    // 使用预处理结果
    ll res = fact_lucas[n];
    res = (res * invfact_lucas[k]) % p;
    res = (res * invfact_lucas[n - k]) % p;
    return (res % p + p) % p;
}

// 卢卡斯定理 计算 C(n, k) mod p (p 是素数)
// 时间复杂度: O(p + log_p(n) * log p) (预处理O(p), 每次递归O(log p)求组合数, 递归log_p(n)次)
// 如果组合数查询是 O(1), 则复杂度是 O(p + log_p(n))
ll lucas(ll n, ll k, ll p) {
    if (k < 0 || k > n) return 0;
    if (k == 0) return 1;
    // 确保预处理完成
    if (p != current_p) {
        precompute_lucas(p);
        if (p != current_p) return -1; // 预处理失败
    }

    ll ni = n % p;
    ll ki = k % p;

    // 计算 C(ni, ki) mod p
    ll c_nk = combinations_small(ni, ki, p);
    if (c_nk == 0) return 0; // 优化：如果某一项是 0，总结果就是 0

    // 递归计算 Lucas(n/p, k/p, p)
    ll remaining = lucas(n / p, k / p, p);
    if (remaining == -1) return -1; // 错误传递

    return (c_nk * remaining) % p;
}

int main() {
    ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);

    ll n = 10, k = 3, p = 7; // C(10, 3) mod 7
    // Lucas(10, 3, 7) = C(10%7, 3%7) * Lucas(10/7, 3/7, 7)
    //                 = C(3, 3) * Lucas(1, 0, 7)
    //                 = 1 * 1 = 1
    // C(10, 3) = 120. 120 mod 7 = 1.
    cout << "Lucas C(" << n << ", " << k << ") mod " << p << " = " << lucas(n, k, p) << endl;

    n = 100, k = 30, p = 13;
    cout << "Lucas C(" << n << ", " << k << ") mod " << p << " = " << lucas(n, k, p) << endl;
    // n=100 = 7*13 + 9  => n1=7, n0=9
    // k=30  = 2*13 + 4  => k1=2, k0=4
    // C(100, 30) = C(9, 4) * C(7, 2) (mod 13)
    // C(9, 4) = (9*8*7*6)/(4*3*2*1) = 9*2*7 = 126. 126 = 9*13 + 9. C(9,4)=9 (mod 13).
    // C(7, 2) = (7*6)/2 = 21. 21 = 1*13 + 8. C(7,2)=8 (mod 13).
    // Result = 9 * 8 = 72. 72 = 5*13 + 7. Result = 7 (mod 13).
    // Let's check the code output.

    return 0;
}

方法三：扩展卢卡斯定理 (Extended Lucas Theorem) 当模数 m 不是素数时，不能直接用卢卡斯定理或简单的逆元法。 扩展卢卡斯定理的思想是：

1. 将模数 m 进行质因数分解：$m = p_1^{a_1} p_2^{a_2} \cdots p_k^{a_k}$。
2. 分别计算 $C(n, k) \pmod{p_i^{a_i}}$ 对于每个 $i$ 的值。
3. 利用中国剩余定理 (CRT) 将这些结果合并，得到最终的 $C(n, k) \pmod m$。

核心在于如何计算 $C(n, k) \pmod{p^a}$ (p 是素数)。 $C(n, k) = \frac{n!}{k!(n-k)!}$。我们不能直接计算模 $p^a$ 的逆元，因为 $k!$ 或 $(n-k)!$ 可能含有因子 p，与 $p^a$ 不互素。

我们需要计算 $n! \pmod{p^a}$。可以将 $n!$ 中所有 p 的因子提出来： $n! = p^{\lfloor n/p \rfloor} (\lfloor n/p \rfloor!) \times (\text{其他因子})$ 这个 "其他因子" 指的是 $1 \times 2 \times \dots \times (p-1) \times (p+1) \times \dots \times (2p-1) \times \dots$ 可以发现，这个 "其他因子" 模 $p^a$ 具有周期性 (周期为 $p^a$)。 我们可以计算出 $n! = p^v \times \text{Rest} \pmod{p^a}$，其中 $v$ 是 $n!$ 中因子 $p$ 的总次数（可以用勒让德定理 Legendre's Formula 快速计算：$v = \sum_{i=1}^{\infty} \lfloor n/p^i \rfloor$），Rest 是除去所有因子 $p$ 后剩余部分的乘积模 $p^a$。

计算 $C(n, k) \pmod{p^a}$ 的步骤：

1. 计算 $n!, k!, (n-k)!$ 中因子 $p$ 的次数 $v_n, v_k, v_{n-k}$。
2. 计算 $n!, k!, (n-k)!$ 除去因子 $p$ 后剩余部分的乘积模 $p^a$，记为 $R_n, R_k, R_{n-k}$。
3. $C(n, k) \pmod{p^a} \equiv p^{v_n - v_k - v_{n-k}} \times R_n \times (R_k)^{-1} \times (R_{n-k})^{-1} \pmod{p^a}$。 这里的逆元 $(R_k)^{-1}$ 和 $(R_{n-k})^{-1}$ 是存在的，因为 $R_k, R_{n-k}$ 与 $p$ 互素，所以它们与 $p^a$ 也互素。可以用 ExGCD 求解模 $p^a$ 的逆元。

计算 $R_n \pmod{p^a}$ (即 $n! / p^{v_n} \pmod{p^a}$) 是 ExLucas 的难点，通常需要递归和快速幂结合周期性来完成。

#include <bits/stdc++.h>

using namespace std;
using ll = long long; // 使用 ll 代替 long long

// --- 需要 ExGCD, qpow, safe_mul, CRT ---

// ExGCD (复用)
ll exgcd(ll a, ll b, ll &x, ll &y) {
    if (b == 0) { x = 1; y = 0; return a; }
    ll d = exgcd(b, a % b, x, y);
    ll temp = x; x = y; y = temp - (a / b) * y;
    return d;
}

// 模逆元 (复用 ExGCD)
ll modInverse(ll a, ll m) {
    ll x, y;
    ll d = exgcd(a, m, x, y);
    if (d != 1) return -1; // 逆元不存在
    return (x % m + m) % m;
}

// 快速幂 (复用)
ll qpow(ll base, ll exp, ll mod) {
    ll res = 1; base %= mod; if(base < 0) base += mod;
    while (exp > 0) {
        if (exp % 2 == 1) res = (res * base) % mod;
        base = (base * base) % mod; exp /= 2;
    }
    return res;
}

// 计算 n! 中因子 p 的个数 (勒让德定理)
ll count_prime_factors(ll n, ll p) {
    ll count = 0;
    while (n > 0) {
        count += n / p;
        n /= p;
    }
    return count;
}

// 计算 n! / p^v mod p^a，即 R_n
// 其中 v = count_prime_factors(n, p)
// p 是素数, pe = p^a
ll factorial_mod_pk(ll n, ll p, ll pe) {
    if (n == 0) return 1;
    ll res = 1;
    // 计算 [1, pe] 中不含因子 p 的数的乘积模 pe
    for (ll i = 1; i <= pe; ++i) {
        if (i % p != 0) {
            res = (res * i) % pe;
        }
    }
    // [1, pe] 的结果可以被快速幂应用到 n / pe 次
    res = qpow(res, n / pe, pe);

    // 计算剩余部分 [1, n % pe] 中不含因子 p 的数的乘积
    for (ll i = 1; i <= n % pe; ++i) {
        if (i % p != 0) {
            res = (res * i) % pe;
        }
    }
    // 递归计算 (n/p)! / p^v' mod pe
    res = (res * factorial_mod_pk(n / p, p, pe)) % pe;
    return res;
}


// 计算 C(n, k) mod p^a (pe = p^a)
ll combinations_mod_pk(ll n, ll k, ll p, ll pe) {
    if (k < 0 || k > n) return 0;
    if (k == 0 || k == n) return 1;

    // 计算 n!, k!, (n-k)! 中 p 的因子数
    ll vn = count_prime_factors(n, p);
    ll vk = count_prime_factors(k, p);
    ll vnk = count_prime_factors(n - k, p);
    ll vp = vn - vk - vnk; // C(n, k) 中 p 的因子数

    // 计算 n!, k!, (n-k)! 除去 p 后的部分模 pe
    ll Rn = factorial_mod_pk(n, p, pe);
    ll Rk = factorial_mod_pk(k, p, pe);
    ll Rnk = factorial_mod_pk(n - k, p, pe);

    // 计算 Rk 和 Rnk 的逆元 mod pe
    ll invRk = modInverse(Rk, pe);
    ll invRnk = modInverse(Rnk, pe);

    if (invRk == -1 || invRnk == -1) {
        // 理论上不应发生，因为 Rk, Rnk 与 p 互素，因此与 pe 互素
        return -1; // Error
    }

    // 计算结果
    ll ans = qpow(p, vp, pe); // p^(vn-vk-vnk) mod pe
    ans = (ans * Rn) % pe;
    ans = (ans * invRk) % pe;
    ans = (ans * invRnk) % pe;

    return (ans % pe + pe) % pe;
}

// 扩展卢卡斯定理
// 计算 C(n, k) mod m
ll extended_lucas(ll n, ll k, ll m) {
     if (m <= 0) return -1; // 模数需为正
     if (k < 0 || k > n) return 0;
     if (m == 1) return 0; // mod 1 结果是 0

     ll temp_m = m;
     vector<ll> a_crt, m_crt; // For CRT

     for (ll p = 2; p * p <= temp_m; ++p) {
         if (temp_m % p == 0) {
             ll pe = 1; // p^a
             while (temp_m % p == 0) {
                 pe *= p;
                 temp_m /= p;
             }
             // 计算 C(n, k) mod pe
             ll res_pk = combinations_mod_pk(n, k, p, pe);
             if (res_pk == -1) return -1; // Error occurred
             a_crt.push_back(res_pk);
             m_crt.push_back(pe);
         }
     }
     // 如果 temp_m 还剩下 > 1，说明它是一个素数
     if (temp_m > 1) {
         ll p = temp_m;
         ll pe = p;
         ll res_pk = combinations_mod_pk(n, k, p, pe);
          if (res_pk == -1) return -1; // Error occurred
         a_crt.push_back(res_pk);
         m_crt.push_back(pe);
     }

     // 使用 CRT 合并结果
     // 注意：这里的 CRT 需要处理模数不互素的情况吗？不需要，因为分解出的 p_i^a_i 是互素的。
     // 可以使用之前的 CRT 或 Excrt 实现。Excrt 更通用。
     // 但这里用 CRT 即可，因为 m_i = p_i^a_i 两两互素。

     // --- CRT (复用之前的 CRT 代码) ---
     // 需要 safe_mul 和 modInverse
     auto crt_result = crt(a_crt, m_crt); // 使用前面定义的 CRT 函数
     if(crt_result.first == -1 || crt_result.first == -2) {
         return -1; // CRT failed
     }
     return crt_result.first; // 返回最小非负解
}


// CRT (复用，确保它在 ExLucas 前或包含在内)
pair<ll, ll> crt(const vector<ll>& a, const vector<ll>& m) {
    int k = a.size();
    if (k == 0) return {0, 1};
    ll M = 1;
    for (ll mi : m) {
        // 简单检查溢出
        if (mi <= 0) return {-1, -1};
        if (M > numeric_limits<ll>::max() / mi) return {-2, -2}; // Overflow check
        M *= mi;
    }
    if (M == 0) return {-1,-1};

    ll x0 = 0;
    for (int i = 0; i < k; ++i) {
        ll Mi = M / m[i];
        ll ti = modInverse(Mi, m[i]);
        if (ti == -1) return {-1, -1};

        // 使用 __int128 或 safe_mul 计算项
         __int128 term128 = (__int128)a[i] * Mi % M * ti % M;
         x0 = (x0 + (ll)term128) % M;
        // 或者用 safe_mul
        // ll term = safe_mul(a[i], Mi, M);
        // term = safe_mul(term, ti, M);
        // x0 = (x0 + term) % M;
    }
    return {(x0 + M) % M, M};
}


int main() {
    ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);

    // Example: C(10, 3) mod 12
    // m = 12 = 2^2 * 3^1
    // p=2, pe=4: C(10, 3) mod 4
    //   v_n(10, 2)=5+2+1=8. v_k(3, 2)=1. v_nk(7, 2)=3+1=4. vp=8-1-4=3.
    //   R_n(10, 2, 4) = fact(10)/2^8 mod 4. fact_mod_pk(10, 2, 4) = 1.
    //   R_k(3, 2, 4) = fact(3)/2^1 mod 4. fact_mod_pk(3, 2, 4) = 3. inv(3,4)=3.
    //   R_nk(7, 2, 4) = fact(7)/2^4 mod 4. fact_mod_pk(7, 2, 4) = 1. inv(1,4)=1.
    //   C(10,3) = 2^3 * 1 * 3 * 1 = 8 * 3 = 24 = 0 (mod 4).
    // p=3, pe=3: C(10, 3) mod 3 (Lucas)
    //   Lucas(10, 3, 3) = C(1, 0) * Lucas(3, 1, 3) = 1 * C(0, 1) * Lucas(1, 0, 3) = 1 * 0 * 1 = 0.
    // CRT: x=0 (mod 4), x=0 (mod 3) => x=0 (mod 12).
    // C(10, 3) = 120. 120 mod 12 = 0. Correct.
    ll n = 10, k = 3, m = 12;
    cout << "ExLucas C(" << n << ", " << k << ") mod " << m << " = " << extended_lucas(n, k, m) << endl;

    // Example: C(6, 3) mod 10
    // m = 10 = 2 * 5
    // p=2, pe=2: C(6, 3) mod 2 (Lucas)
    //   Lucas(6, 3, 2) = C(0,1)*Luc(3,1) = 0 * ... = 0.
    // p=5, pe=5: C(6, 3) mod 5 (Lucas)
    //   Lucas(6, 3, 5) = C(1, 3)*Luc(1,0) = 0 * ... = 0.
    // CRT: x=0 (mod 2), x=0 (mod 5) => x=0 (mod 10).
    // C(6, 3) = 20. 20 mod 10 = 0. Correct.
    n = 6, k = 3, m = 10;
    cout << "ExLucas C(" << n << ", " << k << ") mod " << m << " = " << extended_lucas(n, k, m) << endl;

    return 0;
}

复杂度: ExLucas 的复杂度较高。分解模数 $m$ 需要 $O(\sqrt{m})$。计算 $C(n, k) \pmod{p^a}$ 的 factorial_mod_pk 函数，其递归深度是 $O(\log_p n)$，每层递归涉及 $O(p^a)$ 或 $O(p^a \log n)$ （取决于快速幂实现），总复杂度可能达到 $O(\sum p_i^{a_i} \log n)$ 或更高。CRT 合并需要 $O(k \log M)$，其中 $k$ 是 $m$ 的不同质因子个数，$M=m$。对于竞赛来说，如果 $m$ 的质因子 $p$ 或其指数 $a$ 很大，ExLucas 可能会超时。

六、实战

理论学完了，来看看赛场上的一些实用建议：

1. 模数看清楚！

   • 是 $10^9+7$, $998244353$ (素数) 还是其他数 (可能不是素数)？这决定了求逆元、用 Lucas 还是 ExLucas 等方法的选择。
   • 运算过程中随时取模，尤其是乘法。a * b % mod 是标准操作，但如果 a * b 本身会爆 ll，需要用 safe_mul 或 __int128。
   • 负数取模要变成正数：(a % mod + mod) % mod。

2. GCD 和 ExGCD 是瑞士军刀。

   • 熟练掌握它们的代码模板。
   • 理解 ExGCD 如何解线性丢番图方程 ax + by = c 和模线性方程 ax = c (mod m)。

3. 筛法很重要。

   • 线性筛 O(N) 预处理素数、欧拉函数、莫比乌斯函数（后面会讲）等非常有用。
   • 理解线性筛的核心：每个合数只被其最小素因子筛一次。

4. 欧拉定理/费马小定理降幂。

   • 记住扩展欧拉定理 $a^b \equiv a^{b \pmod{\phi(m)} + \phi(m)} \pmod m$ (当 $b \ge \phi(m)$ 时)。这是处理大指数模非素数的标准方法。
   • 计算 $b \pmod{\phi(m)}$ 时，如果 b 本身就很大 (如字符串)，要小心处理。

5. CRT / Excrt 合并同余方程。

   • 分清 CRT (模数互素) 和 Excrt (模数不互素) 的应用场景和代码。
   • Excrt 合并过程中的模线性方程求解是关键。
   • 极其小心溢出，尤其是 Excrt 中计算 LCM 和更新解的时候。

6. 组合数取模。

   • 模数是素数：
     • N, K 较小 (e.g., $10^6$)：预处理阶乘+逆元，O(1) 查询。
     • N, K 很大，P 较小：Lucas 定理。
   • 模数不是素数：ExLucas 定理。实现复杂，易错，确保理解每个步骤。

7. 数据范围和时间限制。

   • $O(\sqrt{n})$ 的算法通常适用于 $n \le 10^{12} \sim 10^{14}$。
   • $O(N \log \log N)$ 或 $O(N)$ 的筛法适用于 $N \le 10^7 \sim 10^8$。
   • $O(\log n)$ 的算法 (GCD, 快速幂, ExGCD) 几乎适用于所有范围。
   • 注意常数因子，线性筛常数比埃氏筛大。

8. 练习，练习，再练习！

   • 找各种类型的数论题目来做。
   • 尝试自己推导公式，加深理解。
   • 多写代码，调试，形成自己的模板库。

更高阶的主题，如莫比乌斯反演、杜教筛、原根、二次剩余、Pell 方程等，则是在你冲击更高目标时可能需要涉猎的。

但无论走多远，今天我们打下的基础都是至关重要的。掌握好整除、素数、GCD、模运算、同余、逆元、筛法、欧拉函数/定理、CRT、组合数取模这些基本功，你就能解决相当一部分竞赛中的数论问题，也能为学习更高级的算法打下坚实的基础。

数论题往往代码不长，但思维含量高。多思考，多推导，多总结。