📖 第 8.2 章 ⏱️ 约 60 分钟 🥇 Gold

第 8.2 章:素数筛法与积性函数

8.2.0 从一道题出发:数一数有多少素数

先看一道典型的 Gold 风格题目:

🌾 题目:素数普查 农夫 John 想知道 1 到 N(N ≤ 10⁷)之间一共有多少个素数。他还想知道其中每个数的最小质因子,以便快速对任意 a ≤ N 做质因数分解。

输入样例:

N = 20

期望输出:

20 以内有 8 个素数:2 3 5 7 11 13 17 19

先试试你会的方法

思路 1:对每个数 i ∈ [2, N],用试除法判素? 每次 O(√i),总计约 O(N√N)。N = 10⁷ 时 ≈ 3 × 10¹⁰ 次操作——太慢,需要 30 秒以上。

思路 2:换个方向——不去「判定」每个数,而是把合数「筛掉」? 从 2 开始,每遇到一个素数就把它的所有倍数标记为合数。没被标记的就是素数。这就是筛法,也是本章的主角。

而「同时求出最小质因子」「同时求出欧拉 φ 函数」这些进阶需求,会把我们引向欧拉线性筛——每个合数只被它的最小质因子筛掉一次,严格 O(N)。

拆解这道题需要回答的问题

  1. 单个数的素数判定怎么做最省?→ 试除法 O(√n)
  2. 「筛掉合数」的朴素做法及其复杂度?→ 埃拉托斯特尼筛 O(n log log n)
  3. 能否让每个合数只被筛一次?→ 欧拉线性筛 O(n)
  4. 顺手能求出哪些数论函数?→ 欧拉 φ、莫比乌斯 μ 等积性函数
  5. 「x 除以 3 余 2,除以 5 余 3」这类方程组怎么解?→ 中国剩余定理

本章按这条线逐一展开:

小节回答的问题
8.2.1试除法判素:朴素起点
8.2.2埃拉托斯特尼筛
8.2.3欧拉线性筛 + 逐步追踪
8.2.4欧拉 φ 函数与线性筛求 φ
8.2.5积性函数:为什么线性筛能「顺手」求它们
8.2.6中国剩余定理

学习目标:

  • 用试除法在 O(√n) 内判定单个数的素性
  • 实现埃筛与欧拉线性筛,理解后者为何严格 O(n)
  • 用线性筛同时求出 1~n 的欧拉 φ 值
  • 识别积性函数,并用 CRT 合并同余方程组

🔗 前置知识: 模运算、快速幂与模逆元见第 8.1 章。


8.2.1 试除法判定素数(O(√n))

判定单个数 n 是否为素数,最朴素的方法是枚举 [2, n−1] 看有没有因子——O(n)。但有一个关键观察:

若 n 是合数,它必然有一个 ≤ √n 的因子。 因为合数 n = a × b 中,a 和 b 不可能都大于 √n(否则乘积超过 n)。

所以只需枚举到 √n:

// 试除法判素 — O(√n)
bool isPrime(long long n) {
    if (n < 2) return false;
    if (n == 2) return true;
    if (n % 2 == 0) return false;
    for (long long i = 3; i * i <= n; i += 2) {  // 只查奇数,次数减半
        if (n % i == 0) return false;
    }
    return true;
}

两个细节:

  1. 循环条件用 i * i <= n 而不是 i <= sqrt(n)——浮点 sqrt 有精度误差,且每次调用更慢。注意 i 用 long long,否则 i * i 在 i ≈ 46341 时溢出 int。
  2. 先特判 2,再只枚举奇数——把循环次数直接砍半,这在 n ≈ 10¹⁸ 时省的不是一点半点。

试除法在「少量查询」时完全够用。但开篇题要处理 10⁷ 个数,每个都试除就是 3 × 10¹⁰ 次操作。我们需要换武器。


8.2.2 埃拉托斯特尼筛(O(n log log n))

核心思想反转:不问「i 是不是素数」,而是主动「把合数标记出来」。

  • 从 2 开始,2 没被标记 → 2 是素数 → 标记 4, 6, 8, 10, ...
  • 3 没被标记 → 3 是素数 → 标记 9, 12, 15, ...
  • 4 已被标记 → 跳过
  • ……以此类推
const int MAXN = 1e7 + 5;
bool is_composite[MAXN];   // false = 素数(默认值)
vector<int> primes;

void sieve(int n) {
    for (int i = 2; i <= n; i++) {
        if (!is_composite[i]) {           // i 没被筛过 → i 是素数
            primes.push_back(i);
            for (long long j = (long long)i * i; j <= n; j += i)
                is_composite[j] = true;   // 筛掉 i 的倍数
        }
    }
}

三个必问的问题:

Q1:内层为什么从 i * i 开始,而不是 2 * i?

小于 i² 的倍数 2i, 3i, ..., (i−1)·i 已经被更小的素数筛过了。例如 i = 5 时,10、15、20 分别已被 2、3、2 标记,25 才是第一个「首次被 5 筛掉」的数。

Q2:复杂度为什么是 O(n log log n)?

总标记次数 = n/2 + n/3 + n/5 + n/7 + ...(对每个素数 p 标记 n/p 个倍数)≈ n · ln(ln(n))。由素数求和的调和级数估计得出。log log(10⁷) ≈ 2.9,实际跑起来接近线性。

Q3:j 为什么要声明成 long long?

i * i 当 i ≈ 10⁵ 时已达 10¹⁰,用 int 直接溢出成负数,循环行为未定义。

埃筛简单、常数小,是绝大多数题目(「求 10⁶ 以内的素数」)的首选。但它有一个无法回避的缺陷——

一个合数会被重复标记。 比如 12 = 2×6 = 3×4,先被 2 筛一次,又被 3 筛一次;60 = 2²×3×5 会被 2、3、5 各筛一次。合数的质因子越多,被浪费的次数越多。如果题目还需要「最小质因子」或「积性函数值」,埃筛无法提供——而下面的线性筛一举解决所有问题。


8.2.3 欧拉线性筛(O(n))

目标:让每个合数只被筛一次——被它自己的最小质因子筛掉。

做法:外层枚举 i = 2..n;对每个 i,遍历目前已知素数表 primes,用 i * primes[j] 作为被筛掉的合数。关键的两条终止规则:

  1. 若 i * primes[j] > n,越界,跳出;
  2. 若 i % primes[j] == 0,必须跳出——这是正确性的灵魂(下面解释)。
const int MAXN = 1e7 + 5;
int min_prime[MAXN];      // min_prime[i] = i 的最小质因子(素数则是自身)
vector<int> primes;

// 欧拉线性筛:O(n)
void sieve_linear(int n) {
    for (int i = 2; i <= n; i++) {
        if (min_prime[i] == 0) {        // 0 表示还没被筛过 → i 是素数
            min_prime[i] = i;
            primes.push_back(i);
        }
        for (int p : primes) {
            if ((long long)i * p > n) break;   // 规则 1:越界
            min_prime[i * p] = p;              // p 是 i*p 的最小质因子
            if (i % p == 0) break;             // 规则 2:p 已是 i 的因子
        }
    }
}

为什么规则 2 保证「每个合数只筛一次」?

设合数 m = p × i,其中 p 是 m 的最小质因子。线性筛保证 m 只在「外层数 = m/p、内层素数 = p」这一刻被标记。反证:

  • 若 m 还能写成 m = q × i'(q 为更大的素数,i' = m/q < m/p 对应的外层数),那么在处理外层 i' 时,会先遇到素数 p——而 p | i'(因为 p 是 m 的最小质因子且 p < q,p 必然整除 i')——触发规则 2 提前 break,根本轮不到 q。

一句话总结:外层 i 与内层 p 相乘时,只有当 p ≤ i 的最小质因子,乘积 i×p 的最小质因子才恰好是 p。 一旦 p 超过 i 的最小质因子,i×p 就「轮不到 p 来筛」,继续循环只会造成重复。

逐步追踪:n = 20

外层 i是否素数内层标记(按序)break 原因
2✅ 加入4 (=2×2)i%2==0
3✅ 加入6 (=3×2), 9 (=3×3)i%3==0
4❌8 (=4×2)i%2==0(不能用 3:12 = 4×3 的最小质因子是 2,轮不到 3 筛)
5✅ 加入10 (=5×2), 15 (=5×3)5×5=25 越界
6❌12 (=6×2)i%2==0(12 只被 2 筛一次,不像埃筛还会被 3 筛)
7✅ 加入14 (=7×2)21 越界
8❌16 (=8×2)i%2==0
9❌18 (=9×2)27 越界(27 的确该由 3 筛:外层 9 内层 3 会轮到,只是这里 n=20 越界)
10❌20 (=10×2)i%2==0
11~1911,13,17,19 ✅22+ 均越界越界

最终 primes = {2, 3, 5, 7, 11, 13, 17, 19}——8 个素数,与手算一致。每个合数恰好只出现一次在「内层标记」列里。

筛法对比

方法时间复杂度空间何时用
试除法O(√n) 判单个O(1)少量查询、n 高达 10¹⁸
埃氏筛O(n log log n)O(n) bitset 可压到 n/8只需素数表,代码要短
欧拉线性筛O(n)O(n) + 素数表需要最小质因子 / 积性函数

附赠能力:O(log n) 质因数分解

有了 min_prime[],任意 x ≤ n 的分解只需不断除以最小质因子:

// x 的质因数分解 — O(log x)
vector<pair<int,int>> factorize(int x) {
    vector<pair<int,int>> factors;
    while (x > 1) {
        int p = min_prime[x], cnt = 0;
        while (x % p == 0) { x /= p; cnt++; }
        factors.push_back({p, cnt});
    }
    return factors;
}
// factorize(360) = {(2,3), (3,2), (5,1)},即 360 = 2³ × 3² × 5

8.2.4 欧拉 φ 函数

φ(n) 统计 [1, n] 中与 n 互素(gcd(k, n) = 1)的整数个数:

φ(1) = 1
φ(2) = 1  (只有 1)
φ(6) = 2  (1 和 5)
φ(10) = 4 (1, 3, 7, 9)
φ(p) = p − 1    对任意素数 p(1..p−1 全部与 p 互素)
φ(p^k) = p^k − p^(k−1)

积性公式

若 n = p₁^a₁ × p₂^a₂ × ... × pₖ^aₖ,则:

φ(n) = n × (1 − 1/p₁) × (1 − 1/p₂) × ... × (1 − 1/pₖ)

例: φ(10) = 10 × (1 − 1/2) × (1 − 1/5) = 10 × 1/2 × 4/5 = 4 ✓

单值计算:O(√n)

试除法找质因子,边找边乘 (1 − 1/p)。注意用 result -= result / p 而不是浮点除法:

// 欧拉 φ 单值 — O(√n)
long long euler_phi(long long n) {
    long long result = n;
    for (long long p = 2; p * p <= n; p++) {
        if (n % p == 0) {
            while (n % p == 0) n /= p;   // 除尽 p 的所有幂次
            result -= result / p;        // result *= (1 - 1/p)
        }
    }
    if (n > 1) result -= result / n;     // 剩下一个大质因子
    return result;
}

线性筛同时求 φ(1..n):O(n)

单个求 φ 是 O(√n),求 1~n 全部就是 O(n√n)。但 φ 是积性函数(8.2.5 节),可以搭线性筛的便车,在筛的循环里顺手推出所有 φ 值:

const int MAXN = 1e7 + 5;
long long phi[MAXN];
int min_prime[MAXN];
vector<int> primes;

void sieve_phi(int n) {
    for (int i = 2; i <= n; i++) {
        if (min_prime[i] == 0) {          // i 是素数
            min_prime[i] = i;
            phi[i] = i - 1;               // 素数:φ(p) = p - 1
            primes.push_back(i);
        }
        for (int p : primes) {
            if ((long long)i * p > n) break;
            min_prime[i * p] = p;
            if (i % p == 0) {
                // p | i:phi(p*i) = phi(i) * p(p 的幂次增加一层)
                phi[i * p] = phi[i] * p;
                break;
            } else {
                // gcd(i, p) = 1:积性 → phi(p*i) = phi(i) * (p-1)
                phi[i * p] = phi[i] * (p - 1);
            }
        }
    }
    phi[1] = 1;
}

三种情况对应 φ 的三条性质:

情况公式例
i 是素数φ(i) = i − 1φ(5) = 4
p ∤ i(互素)φ(i·p) = φ(i)·φ(p) = φ(i)·(p−1)φ(6) = φ(2)·φ(3) = 1×2
p | i(p 已是 i 的因子)φ(i·p) = φ(i)·pφ(4) = φ(2)·2 = 2;φ(8) = φ(4)·2 = 4

验证: sieve_phi(10) 后 φ = [1, 1, 2, 2, 4, 4, 2, 6, 4, 6, 4],其中 φ(10) = 4 ✓

φ 的经典应用

  1. 欧拉定理(费马小定理的推广):gcd(a, n) = 1 时,a^φ(n) ≡ 1 (mod n)。当 n 是素数时 φ(n) = n−1 退化为费马小定理——这正是第 8.1 章模逆元的理论基础。
  2. 降幂: 计算 a^b mod n(b 极大)时,若 gcd(a,n)=1,指数可以对 φ(n) 取模。
  3. 计数: 统计 1 ≤ a ≤ b ≤ n 且 gcd(a,b) = 1 的对数 = 1 + Σᵢ₌₂ⁿ φ(i)。
  4. Burnside 项链计数: 答案中出现 Σ φ(d),d 取遍 n 的因子。

8.2.5 积性函数:线性筛的通用引擎

线性筛之所以能「顺手」求 φ,是因为 φ 是积性函数:

定义: 若函数 f 满足 f(ab) = f(a)·f(b) 当 gcd(a, b) = 1,则称 f 为积性函数。 若对所有 a, b 都成立(不要求互素),则称完全积性函数。

常见积性函数速览

函数定义完全积性?备注
恒等函数 id(n) = n—✅
常数函数 1(n) = 1—✅
欧拉 φ[1,n] 中与 n 互素的个数❌φ(4) = 2 ≠ φ(2)·φ(2) = 1
莫比乌斯 μμ(n) = 0(有平方因子);(−1)^k(k 个不同质因子)❌容斥的代数化
约数个数 d(n)n 的因子个数❌n=Πpᵃⁱ 时 d = Π(aᵢ+1)
约数和 σ(n)n 的因子之和❌

为什么积性函数都能用线性筛求?

回顾线性筛标记每个合数 m = i × p 的两个分支:

  • p ∤ i(gcd(i, p) = 1):由积性,f(m) = f(i)·f(p)——两个已知值相乘;
  • p | i:m 与 i 只差一个 p 的幂次,通常存在只依赖 p 的简单递推(如 φ(i·p) = φ(i)·p)。

再加上 f 在素数处的取值往往有封闭形式(φ(p) = p−1、μ(p) = −1),三条信息合起来,线性筛就能以 O(n) 推出整个函数值表。模板结构完全同 8.2.4,只改三个分支的公式。

在 USACO Gold 中最常考的就是 φ;μ 与容斥结合(莫比乌斯反演)多见于 Platinum。


8.2.6 中国剩余定理(CRT)

「一个数除以 3 余 2,除以 5 余 3,除以 7 余 2,问这个数是多少?」——《孙子算经》(约公元 4 世纪)

中国剩余定理(CRT): 若模数 m₁, m₂, ..., mₖ 两两互素,则同余方程组

x ≡ r₁ (mod m₁)
x ≡ r₂ (mod m₂)
...
x ≡ rₖ (mod mₖ)

在模 M = m₁m₂...mₖ 意义下存在唯一解。

合并思路(两方程)

对 x ≡ r₁ (mod m₁)、x ≡ r₂ (mod m₂),gcd(m₁, m₂) = 1:

  1. 写 x = r₁ + m₁·k;
  2. 代入第二式:m₁·k ≡ r₂ − r₁ (mod m₂);
  3. 两边乘 m₁ 的模逆元:k ≡ (r₂ − r₁)·m₁⁻¹ (mod m₂);
  4. 取最小的 k,得 x = r₁ + m₁·k,落在 [0, m₁m₂) 内。
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

// 快速幂 — 见第 8.1 章
ll power(ll a, ll n, ll p) {
    a %= p;
    ll result = 1;
    while (n > 0) {
        if (n & 1) result = result * a % p;
        a = a * a % p;
        n >>= 1;
    }
    return result;
}

// 模逆元(m 为素数时用费马小定理;这里 m1、m2 互素且各自是素数或与 m2 互素即可)
// 通用做法:扩展欧几里得。此处给 extgcd 版本,对任意互素模数都正确。
ll ext_gcd(ll a, ll b, ll &x, ll &y) {
    if (b == 0) { x = 1; y = 0; return a; }
    ll x1, y1;
    ll g = ext_gcd(b, a % b, x1, y1);
    x = y1;
    y = x1 - (a / b) * y1;
    return g;
}

ll mod_inv(ll a, ll m) {          // a 与 m 互素
    ll x, y;
    ext_gcd(a, m, x, y);
    return (x % m + m) % m;
}

// 两方程 CRT(要求 gcd(m1, m2) = 1),返回 [0, m1*m2) 内的唯一解
ll crt(ll r1, ll m1, ll r2, ll m2) {
    ll k = (r2 - r1 % m2 + m2) % m2 * mod_inv(m1 % m2, m2) % m2;
    return r1 + m1 * k;
}

// 多方程:迭代两两合并,返回 {x, lcm};无解返回 {-1, -1}
pair<ll, ll> crt_multi(vector<ll> r, vector<ll> m) {
    ll cur_r = r[0], cur_m = m[0];
    for (int i = 1; i < (int)r.size(); i++) {
        // 广义合并(允许模数不互素)
        ll g = __gcd(cur_m, m[i]);
        if ((r[i] - cur_r) % g != 0) return {-1, -1};   // 无解
        ll lcm = cur_m / g * m[i];
        ll diff = (r[i] - cur_r) / g;
        ll m2g = m[i] / g;
        ll k = diff % m2g * mod_inv(cur_m / g % m2g, m2g) % m2g;
        cur_r = ((cur_r + cur_m % lcm * ((k % lcm + lcm) % lcm)) % lcm + lcm) % lcm;
        // 简化写法:cur_r = (cur_r + cur_m * k) % lcm(cur_m*k 可能需 __int128)
        cur_m = lcm;
    }
    return {cur_r, cur_m};
}

int main() {
    // 孙子算经:x ≡ 2 (mod 3), x ≡ 3 (mod 5), x ≡ 2 (mod 7)
    auto [x, M] = crt_multi({2, 3, 2}, {3, 5, 7});
    cout << "x = " << x << " (mod " << M << ")" << endl;   // x = 23 (mod 105)
    return 0;
}

手算验证: 23 mod 3 = 2 ✓,23 mod 5 = 3 ✓,23 mod 7 = 2 ✓。

USACO CRT 题型模式

「某事件每 A 天发生一次(今天余 r₁ 天)、另一事件每 B 天一次……它们何时首次同时发生?」

就是解 x ≡ r₁ (mod A)、x ≡ r₂ (mod B)……用 crt_multi 迭代合并即可。注意模数不互素时解可能不存在——广义版本已通过 gcd 判断。


⚠️ 常见错误

错误后果修复
线性筛内层先判 i % p == 0 再记录漏筛(i×p 本该由 p 筛掉)先记 min_prime[i*p] = p,再判 i % p == 0 决定 break,顺序不能反
线性筛 i * p 用 int 计算i ≈ 10⁵ 时溢出为负,数组越界写 (long long)i * p > n
埃筛内层从 2*i 起且 j 为 int溢出 + 重复标记浪费从 i*i 起,j 用 long long
phi[1] 忘记初始化phi[1] = 0,但正确值是 1循环从 2 开始,显式补 phi[1] = 1
φ 递推中 p | i 分支写成 phi[i]*(p-1)混淆两个分支(如 φ(4) 会得 2×1=2 碰巧对,但 φ(8) 得 φ(4)×1 = 2 ≠ 4)牢记:p ∤ i 时乘 (p−1),p | i 时乘 p
CRT 中 (r2 - r1) 为负直接取模负数取模结果为负(r2 - r1 % m2 + m2) % m2
CRT 结果 m1 * k 溢出m₁·m₂ > 9.2×10¹⁸ 时爆 long long用 __int128 做中间乘法

📋 章节小结

工具复杂度一句话
试除法O(√n)/次单个数判定,n 可达 10¹⁸
埃氏筛O(n log log n)简单快,但合数被重复筛
欧拉线性筛O(n)每个合数只被最小质因子筛一次
min_prime 分解O(log n)/次筛后免费获得的质因数分解
φ 单值 / 筛求O(√n) / O(n)φ(n) = n·Π(1−1/p),积性
积性函数通筛O(n)三分支:素数处 / 互素处 / 幂次处
CRTO(k log M)同余方程组合并,互素唯一解

🔗 与后续章节的联系

  • 第 8.3 章(组合数学): 容斥 + 素数筛是「整除计数」题的黄金组合;Burnside 引理(用到 φ)出现在项链计数中。
  • 第 8.1 章(数论基础): 欧拉定理 a^φ(n) ≡ 1 是费马小定理的推广,逆元求法由此而来。
  • 第 8.4 章(矩阵快速幂): 与本章的降幂技巧配合处理「超大指数 + 取模」。

🏋️ 练习题

🟢 基础

8.2-P1. 第 10⁷ 个以内的素数个数
用欧拉线性筛统计并输出 1 到 10⁷ 中素数的个数,以及 10⁷ 以内最大的素数。

提示

sieve_linear(10000000) 后输出 primes.size()(应为 620648 与……直接跑)和 primes.back()。注意数组开 1e7 + 5。

8.2-P2. 质因数分解器
给定 T(≤ 10⁵)个整数 aᵢ(≤ 10⁷),对每个输出形如 2^3 * 3^2 * 5 的质因数分解。

提示

先 sieve_linear(1e7) 建好 min_prime,每个查询用 8.2.3 的 factorize,单次 O(log a)。总复杂度 O(N + T log N)。

🟡 中等

8.2-P3. 互素对计数
统计满足 1 ≤ a ≤ b ≤ n 的互素对 (a, b) 个数,n ≤ 10⁶,答案 mod 10⁹+7。

提示

互素对数 = 1 + Σᵢ₌₂ⁿ φ(i)("+1" 是 (1,1) 这一对)。用 sieve_phi(n) 一遍扫出。用 long long 累加,n = 10⁶ 时答案约 2×10¹¹。

🔴 困难

8.2-P4. 三个时钟 (CRT)
某星球有三个周期事件:事件 A 每 4 小时于第 3 小时发生(x ≡ 3 mod 4),事件 B 满足 x ≡ 5 (mod 7),事件 C 满足 x ≡ 1 (mod 9)。求三者首次同时发生的时刻 x(模 252 意义下的最小非负解)。

提示

crt_multi({3, 5, 1}, {4, 7, 9})。手算可以先合并前两个:x ≡ 3 (mod 4) 与 x ≡ 5 (mod 7) → x ≡ 19 (mod 28);再与 x ≡ 1 (mod 9) 合并 → 验证答案对 4、7、9 的余数。