数论算法详解：从CSP-S到IOI#

数论是算法竞赛中的重要分支，从CSP-S到IOI级别的比赛都经常出现数论相关题目。本文将系统地介绍数论算法的核心知识点，从基础到进阶，帮助读者建立完整的数论知识体系。

一、质数与筛法#

1.1 质数的定义与判断#

质数（素数）是指在大于1的自然数中，除了1和它本身以外不再有其他因数的数。

朴素判断法：时间复杂度 O(√n)

1
bool isPrime(int n) {
2
    if (n < 2) return false;
3
    for (int i = 2; i * i <= n; i++) {
4
        if (n % i == 0) return false;
5
    }
6
    return true;
7
}

1.2 埃拉托斯特尼筛法（Eratosthenes筛）#

算法思想：从小到大枚举每个质数，标记其所有倍数为合数。

时间复杂度：O(n log log n)

1
const int MAXN = 1e6 + 5;
2
bool isPrime[MAXN];
3
vector<int> primes;
4

5
void sieveOfEratosthenes(int n) {
6
    fill(isPrime, isPrime + n + 1, true);
7
    isPrime[0] = isPrime[1] = false;
8

9
    for (int i = 2; i <= n; i++) {
10
        if (isPrime[i]) {
11
            primes.push_back(i);
12
            for (long long j = (long long)i * i; j <= n; j += i) {
13
                isPrime[j] = false;
14
            }
15
        }
16
    }
17
}

1.3 线性筛法（欧拉筛）#

算法思想：每个合数只被其最小质因子筛掉一次，确保线性时间复杂度。

时间复杂度：O(n)

1
const int MAXN = 1e7 + 5;
2
bool isPrime[MAXN];
3
int primes[MAXN], cnt = 0;
4

5
void linearSieve(int n) {
6
    fill(isPrime, isPrime + n + 1, true);
7
    isPrime[0] = isPrime[1] = false;
8

9
    for (int i = 2; i <= n; i++) {
10
        if (isPrime[i]) primes[cnt++] = i;
11

12
        for (int j = 0; j < cnt && (long long)i * primes[j] <= n; j++) {
13
            isPrime[i * primes[j]] = false;
14
            if (i % primes[j] == 0) break;  // 关键：保证每个数只被最小质因子筛掉
15
        }
16
    }
17
}

核心理解：当 i % primes[j] == 0 时，primes[j] 是 i 的最小质因子，也是 i * primes[j] 的最小质因子，此时应该停止，避免重复筛选。

二、最大公约数与最小公倍数#

2.1 欧几里得算法（辗转相除法）#

定理：gcd(a, b) = gcd(b, a mod b)

证明：设 d = gcd(a, b)，则 d|a 且 d|b。因为 a = bq + r (其中 r = a mod b)，所以 d|r。反之，若 d|b 且 d|r，则 d|a。因此 gcd(a, b) = gcd(b, r)。

时间复杂度：O(log min(a, b))

1
int gcd(int a, int b) {
2
    return b == 0 ? a : gcd(b, a % b);
3
}
4

5
// 最小公倍数
6
int lcm(int a, int b) {
7
    return a / gcd(a, b) * b;  // 先除后乘防止溢出
8
}

2.2 扩展欧几里得算法#

问题：求解方程 ax + by = gcd(a, b) 的整数解。

递归思路：

当 b = 0 时，gcd(a, 0) = a，此时 x = 1, y = 0
否则求解 bx’ + (a mod b)y’ = gcd(b, a mod b)
因为 a mod b = a - ⌊a/b⌋ × b，代入得：bx’ + (a - ⌊a/b⌋ × b)y’ = gcd(a, b)
整理得：ay’ + b(x’ - ⌊a/b⌋ × y’) = gcd(a, b)
因此 x = y’, y = x’ - ⌊a/b⌋ × y’

1
int exgcd(int a, int b, int &x, int &y) {
2
    if (b == 0) {
3
        x = 1, y = 0;
4
        return a;
5
    }
6
    int d = exgcd(b, a % b, y, x);
7
    y -= a / b * x;
8
    return d;
9
}

应用：求模逆元

若 gcd(a, p) = 1，则 a 在模 p 意义下的逆元为 x (mod p)，其中 x 满足 ax ≡ 1 (mod p)。

1
int modInverse(int a, int p) {
2
    int x, y;
3
    int d = exgcd(a, p, x, y);
4
    if (d != 1) return -1;  // 不存在逆元
5
    return (x % p + p) % p;
6
}

三、快速幂与矩阵快速幂#

3.1 快速幂算法#

问题：计算 a^n mod p

算法思想：将指数 n 二进制分解，利用 a^n = a^(2^k1) × a^(2^k2) × … 的性质。

时间复杂度：O(log n)

1
long long quickPow(long long a, long long n, long long p) {
2
    long long res = 1;
3
    a %= p;
4
    while (n > 0) {
5
        if (n & 1) res = res * a % p;
6
        a = a * a % p;
7
        n >>= 1;
8
    }
9
    return res;
10
}

3.2 矩阵快速幂#

应用：加速递推关系的计算，如斐波那契数列。

1
struct Matrix {
2
    long long a[2][2];
3
    Matrix() { memset(a, 0, sizeof(a)); }
4

5
    Matrix operator * (const Matrix &b) const {
6
        Matrix c;
7
        for (int i = 0; i < 2; i++) {
8
            for (int j = 0; j < 2; j++) {
9
                for (int k = 0; k < 2; k++) {
10
                    c.a[i][j] = (c.a[i][j] + a[i][k] * b.a[k][j]) % MOD;
11
                }
12
            }
13
        }
14
        return c;
15
    }
16
};
17

18
Matrix matPow(Matrix a, long long n) {
19
    Matrix res;
20
    res.a[0][0] = res.a[1][1] = 1;  // 单位矩阵
21
    while (n > 0) {
22
        if (n & 1) res = res * a;
23
        a = a * a;
24
        n >>= 1;
25
    }
26
    return res;
27
}
28

29
// 斐波那契数列第n项
30
long long fibonacci(long long n) {
31
    if (n <= 2) return 1;
32
    Matrix m;
33
    m.a[0][0] = m.a[0][1] = m.a[1][0] = 1;
34
    m.a[1][1] = 0;
35
    Matrix res = matPow(m, n - 2);
36
    return (res.a[0][0] + res.a[0][1]) % MOD;
37
}

四、模运算与逆元#

4.1 模运算的基本性质#

(a + b) mod p = ((a mod p) + (b mod p)) mod p
(a - b) mod p = ((a mod p) - (b mod p) + p) mod p
(a × b) mod p = ((a mod p) × (b mod p)) mod p
(a / b) mod p = (a × b^(-1)) mod p （b^(-1) 为 b 的逆元）

4.2 费马小定理#

定理：若 p 为质数，gcd(a, p) = 1，则 a^(p-1) ≡ 1 (mod p)

推论：a^(p-2) ≡ a^(-1) (mod p)

应用：当模数为质数时，快速求逆元

1
long long modInverse(long long a, long long p) {
2
    return quickPow(a, p - 2, p);
3
}

4.3 线性求逆元#

问题：求 1, 2, …, n 在模 p 下的逆元

递推公式：inv[i] = (p - p/i) × inv[p % i] % p

证明：设 p = ki + r (r = p % i)，则 ki + r ≡ 0 (mod p)，两边同时乘以 i^(-1) × r^(-1)，得 i^(-1) ≡ -k × r^(-1) (mod p)

1
long long inv[MAXN];
2

3
void calcInverse(int n, int p) {
4
    inv[1] = 1;
5
    for (int i = 2; i <= n; i++) {
6
        inv[i] = (p - p / i) * inv[p % i] % p;
7
    }
8
}

时间复杂度：O(n)

4.4 中国剩余定理（CRT）#

问题：求解同余方程组

1
x ≡ a1 (mod m1)
2
x ≡ a2 (mod m2)
3
...
4
x ≡ ak (mod mk)

其中 m1, m2, …, mk 两两互质。

算法：

令 M = m1 × m2 × … × mk
令 Mi = M / mi
求 ti 满足 Mi × ti ≡ 1 (mod mi)
答案 x = Σ(ai × Mi × ti) mod M

1
long long CRT(int n, long long a[], long long m[]) {
2
    long long M = 1, ans = 0;
3
    for (int i = 0; i < n; i++) M *= m[i];
4

5
    for (int i = 0; i < n; i++) {
6
        long long Mi = M / m[i];
7
        long long ti = modInverse(Mi, m[i]);
8
        ans = (ans + a[i] * Mi % M * ti % M) % M;
9
    }
10
    return (ans % M + M) % M;
11
}

五、欧拉函数#

5.1 定义与性质#

定义：φ(n) 表示小于等于 n 且与 n 互质的正整数个数。

性质：

φ(1) = 1
若 p 为质数，φ(p) = p - 1
若 p 为质数，φ(p^k) = p^k - p^(k-1) = p^(k-1) × (p - 1)
若 gcd(m, n) = 1，φ(m × n) = φ(m) × φ(n)（积性函数）

5.2 单个数的欧拉函数#

公式：若 n = p1^k1 × p2^k2 × … × pm^km，则 φ(n) = n × (1 - 1/p1) × (1 - 1/p2) × … × (1 - 1/pm)

1
long long phi(long long n) {
2
    long long res = n;
3
    for (long long i = 2; i * i <= n; i++) {
4
        if (n % i == 0) {
5
            res = res / i * (i - 1);
6
            while (n % i == 0) n /= i;
7
        }
8
    }
9
    if (n > 1) res = res / n * (n - 1);
10
    return res;
11
}

时间复杂度：O(√n)

5.3 线性筛求欧拉函数#

1
long long phi[MAXN];
2
int primes[MAXN], cnt = 0;
3
bool isPrime[MAXN];
4

5
void phiTable(int n) {
6
    phi[1] = 1;
7
    fill(isPrime, isPrime + n + 1, true);
8

9
    for (int i = 2; i <= n; i++) {
10
        if (isPrime[i]) {
11
            primes[cnt++] = i;
12
            phi[i] = i - 1;
13
        }
14

15
        for (int j = 0; j < cnt && (long long)i * primes[j] <= n; j++) {
16
            isPrime[i * primes[j]] = false;
17
            if (i % primes[j] == 0) {
18
                phi[i * primes[j]] = phi[i] * primes[j];
19
                break;
20
            } else {
21
                phi[i * primes[j]] = phi[i] * (primes[j] - 1);
22
            }
23
        }
24
    }
25
}

5.4 欧拉定理#

定理：若 gcd(a, n) = 1，则 a^φ(n) ≡ 1 (mod n)

应用：降幂公式，a^b ≡ a^(b mod φ(n)) (mod n)（当 b ≥ φ(n) 时）

六、组合数学#

6.1 排列组合基础#

排列数：A(n, m) = n!/(n-m)!

组合数：C(n, m) = n!/(m!(n-m)!)

递推公式：C(n, m) = C(n-1, m-1) + C(n-1, m)（杨辉三角）

1
// 预处理组合数（小范围）
2
long long C[MAXN][MAXN];
3

4
void initCombination(int n) {
5
    for (int i = 0; i <= n; i++) {
6
        C[i][0] = 1;
7
        for (int j = 1; j <= i; j++) {
8
            C[i][j] = (C[i-1][j-1] + C[i-1][j]) % MOD;
9
        }
10
    }
11
}
12

13
// Lucas定理（p为质数，n,m很大）
14
long long lucas(long long n, long long m, long long p) {
15
    if (m == 0) return 1;
16
    return C[n % p][m % p] * lucas(n / p, m / p, p) % p;
17
}

6.2 卡特兰数#

定义：第 n 个卡特兰数 Cat(n) = C(2n, n) / (n + 1)

应用场景：

n 对括号的合法匹配方案数
n 个节点构成的二叉树个数
n × n 方格从 (0,0) 走到 (n,n) 且不穿过对角线的路径数

递推公式：Cat(n) = Σ Cat(i) × Cat(n-1-i)，其中 i = 0 to n-1

1
long long catalan[MAXN];
2

3
void initCatalan(int n) {
4
    catalan[0] = catalan[1] = 1;
5
    for (int i = 2; i <= n; i++) {
6
        catalan[i] = catalan[i-1] * (4 * i - 2) / (i + 1);
7
    }
8
}

6.3 斯特林数#

第一类斯特林数 S(n, k)：n 个元素分成 k 个环的方案数

递推：S(n, k) = S(n-1, k-1) + (n-1) × S(n-1, k)

第二类斯特林数 S(n, k)：n 个元素分成 k 个非空集合的方案数

递推：S(n, k) = S(n-1, k-1) + k × S(n-1, k)

七、莫比乌斯反演#

7.1 莫比乌斯函数#

定义：

μ(1) = 1
μ(n) = (-1)^k，若 n 为 k 个不同质数的乘积
μ(n) = 0，若 n 含有平方因子

性质：Σ μ(d) = [n == 1]（d|n）

7.2 莫比乌斯反演#

定理：若 F(n) = Σ f(d)（d|n），则 f(n) = Σ μ(d) × F(n/d)（d|n）

1
int mu[MAXN];
2

3
void initMobius(int n) {
4
    mu[1] = 1;
5
    for (int i = 2; i <= n; i++) {
6
        if (isPrime[i]) {
7
            primes[cnt++] = i;
8
            mu[i] = -1;
9
        }
10
        for (int j = 0; j < cnt && (long long)i * primes[j] <= n; j++) {
11
            if (i % primes[j] == 0) {
12
                mu[i * primes[j]] = 0;
13
                break;
14
            } else {
15
                mu[i * primes[j]] = -mu[i];
16
            }
17
        }
18
    }
19
}

经典应用：求 Σ Σ [gcd(i, j) == 1]（1 ≤ i ≤ n, 1 ≤ j ≤ m）

八、原根与离散对数#

8.1 阶的概念#

定义：满足 a^r ≡ 1 (mod p) 的最小正整数 r 称为 a 模 p 的阶，记作 ord_p(a)。

8.2 原根#

定义：若 ord_p(g) = φ(p)，则称 g 为模 p 的原根。

性质：若 g 是模 p 的原根，则 g^0, g^1, …, g^(φ(p)-1) 模 p 的余数恰好是所有与 p 互质的数。

应用：离散对数、指标运算

8.3 BSGS算法#

问题：求解 a^x ≡ b (mod p)（Baby Step Giant Step）

算法思想：设 x = im - j，其中 m = ⌈√p⌉，枚举 j，用哈希表存储 a^j，然后枚举 i 查找。

时间复杂度：O(√p)

1
map<long long, int> mp;
2

3
long long BSGS(long long a, long long b, long long p) {
4
    mp.clear();
5
    a %= p; b %= p;
6
    if (b == 1) return 0;
7

8
    int m = sqrt(p) + 1;
9
    long long val = b;
10
    for (int j = 0; j < m; j++) {
11
        mp[val] = j;
12
        val = val * a % p;
13
    }
14

15
    long long base = quickPow(a, m, p);
16
    val = 1;
17
    for (int i = 1; i <= m; i++) {
18
        val = val * base % p;
19
        if (mp.count(val)) {
20
            return i * m - mp[val];
21
        }
22
    }
23
    return -1;  // 无解
24
}

九、Miller-Rabin素性测试#

9.1 算法原理#

基于费马小定理的改进：对于奇质数 p，p - 1 = 2^s × d（d为奇数），若 a^d ≢ 1 (mod p)，则必存在 r ∈ [0, s-1] 使得 a^(2^r × d) ≡ -1 (mod p)。

时间复杂度：O(k log^3 n)，k 为测试轮数

1
bool millerRabin(long long n, int times = 10) {
2
    if (n < 2) return false;
3
    if (n == 2 || n == 3) return true;
4
    if (n % 2 == 0) return false;
5

6
    // n - 1 = 2^s × d
7
    long long d = n - 1, s = 0;
8
    while (d % 2 == 0) {
9
        d /= 2;
10
        s++;
11
    }
12

13
    auto check = [&](long long a) {
14
        long long x = quickPow(a, d, n);
15
        if (x == 1 || x == n - 1) return true;
16
        for (int i = 0; i < s - 1; i++) {
17
            x = (__int128)x * x % n;
18
            if (x == n - 1) return true;
19
        }
20
        return false;
21
    };
22

23
    // 测试多个底数
24
    for (int i = 0; i < times; i++) {
25
        long long a = rand() % (n - 3) + 2;
26
        if (!check(a)) return false;
27
    }
28
    return true;
29
}

十、Pollard-Rho大数分解#

10.1 算法思想#

利用生日悖论的概率方法，寻找 n 的非平凡因子。通过伪随机函数 f(x) = (x^2 + c) mod n 生成序列，利用 Floyd 判圈算法找到重复。

时间复杂度：期望 O(n^(1/4))

1
long long pollardRho(long long n) {
2
    if (n % 2 == 0) return 2;
3

4
    long long c = rand() % (n - 1) + 1;
5
    long long x = rand() % n, y = x;
6
    long long d = 1;
7

8
    auto f = [&](long long val) {
9
        return ((__int128)val * val + c) % n;
10
    };
11

12
    while (d == 1) {
13
        x = f(x);
14
        y = f(f(y));
15
        d = gcd(abs(x - y), n);
16
    }
17
    return d == n ? pollardRho(n) : d;
18
}
19

20
// 完全分解
21
void factorize(long long n, vector<long long> &factors) {
22
    if (n == 1) return;
23
    if (millerRabin(n)) {
24
        factors.push_back(n);
25
        return;
26
    }
27
    long long p = pollardRho(n);
28
    factorize(p, factors);
29
    factorize(n / p, factors);
30
}

十一、IOI级别应用题#

11.1 例题1：约数个数问题#

问题：给定 n 个正整数 a1, a2, …, an，求它们的乘积的约数个数模 10^9+7。

思路：分解质因数，利用约数个数公式 (k1+1)(k2+1)…(km+1)，其中 ki 是第 i 个质因子的指数。

11.2 例题2：同余方程求解#

问题：求解 ax ≡ b (mod p)，其中 p 不一定是质数。

思路：

求 d = gcd(a, p)
若 d ∤ b，无解
否则将方程化简为 (a/d)x ≡ (b/d) (mod p/d)
用扩展欧几里得求解

11.3 例题3：线性递推优化#

问题：求 f(n) = c1×f(n-1) + c2×f(n-2) + … + ck×f(n-k) 的第 n 项。

思路：构造转移矩阵，用矩阵快速幂求解，时间复杂度 O(k^3 log n)。

11.4 例题4：多项式模运算#

问题：给定多项式 P(x) 和 Q(x)，求 P(x) mod Q(x)。

思路：利用快速幂思想，结合多项式乘法（FFT/NTT优化），可在 O(n log^2 n) 内完成。

总结#

数论算法是算法竞赛的核心内容之一，从基础的质数筛、GCD到高级的Miller-Rabin、Pollard-Rho，每个知识点都有其独特的应用场景。掌握这些算法需要：

理解数学原理：证明和推导能加深理解
熟练代码实现：模板要能快速默写
注意边界条件：模运算、溢出、特殊情况
灵活应用组合：实际问题往往需要多个知识点结合

建议学习路径：

CSP-S/提高组：质数筛、GCD、快速幂、逆元、欧拉函数
省选/NOI：扩展欧几里得、CRT、组合数学、莫比乌斯反演
IOI/国际赛：原根、BSGS、Miller-Rabin、Pollard-Rho

持续练习，不断总结，相信你能在数论这个美妙的领域取得优异成绩！

参考资料#

《算法竞赛进阶指南》
《具体数学》(Concrete Mathematics)
《初等数论》(Elementary Number Theory)
OI Wiki: https://oi-wiki.org/math/number-theory/
Codeforces 数论专题题单