数论复习笔记

· · 算法·理论

质数

定义一个正整数 p 为质数，即它不存在 1 和 p 以外的约数。1 不是质数。

质数的判定

枚举法

可以枚举 [1,\sqrt{n}] 中的整数，判定其是否为 p 的约数。因为 a \mid p 时 \dfrac{p}{a}\mid p，所以无需枚举到 n。复杂度为 O(\sqrt{n})。

这里给出一个 \dfrac{1}{6} 常数的实现，筛掉 2 和 3 的倍数：

inline bool isprime(int x) {
    if (x <= 1) return false;
    if (x == 2 || x == 3) return true;
    if (x % 6 != 1 && x % 6 != 5) return false;
    for (int i = 5; i * i <= x; i += 6) {
        if (x % i == 0 || x % (i + 2) == 0) return false;
    }
    return true;
}

Miller-Rabin 法

先引入二次探测定理，即：若 p 为奇质数，则 x^2 \equiv 1 \pmod p 的解为 x \equiv 1 或 x \equiv p-1。

简证：由 x^2 \equiv 1 \pmod p 得 (x-1)(x+1) \equiv 0 \pmod p，即可得出。

接下来使用费马小定理：若 p 为质数且 \gcd(a,p)=1，则 a^{p-1} \equiv 1 \pmod p。这个结论的讲解可到文章下部。

设 p-1=u \times 2^t，随机一个 a 值并求得 v=a^u \bmod p，然后执行 t 次 v \gets v^2 \bmod p，检查是否满足 a^{p-1} \equiv 1 \pmod p。

在 long long 范围内，只需选取 a \in \{2,325,9375,28178,450775,9780504,1795265022\} 即可保证正确性。复杂度 O(\log n)。

inline long long power(long long a, long long b, long long p) {
    long long res = 1;
    for (a %= p; b; b >>= 1) {
        if (b & 1) res = res * a % p;
        a = a * a % p;
    }
    return res;
}

const long long BASE[] = {2, 325, 9375, 28178, 450775, 9780504, 1795265022LL};

inline bool miller_rabin(long long x) {
    if (x < 3 || !(x & 1)) return x == 2;
    long long u = x - 1, t = 0;
    while (!(u & 1)) u >>= 1, t++;
    for (long long a : BASE) {
        long long v = power(a, u, x);
        if (a % x == 0 || v == 1) continue;
        long long s = 0;
        for (; s < t; s++) {
            if (v == x - 1) break;
            v = v * v % x;
        }
        if (s == t) return false;
    }
    return true;
}

威尔逊定理

若 p 为质数，则

(p-1)! \equiv -1 \pmod p

其逆定理也成立。即若 (p-1)! \equiv -1 \pmod p，则 p 为质数。

算术基本定理

任何一个合数 n 可以唯一分解成有限个质数的乘积。

存在性：用反证法，假设 n 是最小的不能被分解的合数，则存在 n=ab，若 a,b 都可分解，则 n 可以被分解；若 a,b 有不可分解的数，则 a,b 才是最小的数，矛盾。

唯一性：用反证法，假设 n 是最小的存在两种分解的合数，如果 n 存在两种分解 n={p_1}^{a_1} {p_2}^{a_2} \cdots ={q_1}^{b_1} {q_2}^{b_2} \cdots，则 p_1 | {q_1}^{b_1} {q_2}^{b_2} \cdots，也就是 {q_1}^{b_1} {q_2}^{b_2} \cdots 中有一个 {q_i}{b_i} 可以整除 p_1，故 p_1=q_i，同除 p_1，则 {p_2}^{a_2} \cdots 也是存在两种分解的合数，矛盾。

推论

设 n 可以质因数分解为 n=p_1^{c_1} p_2^{c_2} p_3^{c_3} \cdots p_m^{c_m}。

推论1：n 的正约数个数为

(c_1+1)(c_2+1)(c_3+1)\cdots(c_m+1)=\prod_{i=1}^{m} (c_i+1)

推论2：n 的所有正约数和为

(1+p_1+p_1^2+\cdots+p_1^{c_1})\cdots(1+p_m+p_m^2+\cdots+p_m^{c_m})=\prod_{i=1}^{m} [\sum_{j=0}^{c_i} {p_i}^j]

分解质因数

将一个数分解质因数最简单的方法就是试除法，复杂度 O(\sqrt{n})。

int p[N], c[N];
int decompose(int n) {
    int m = 0;
    for (int i = 2; i * i <= n; i++) {
        if (n % i) continue;
        p[++m] = i, c[m] = 0;
        while (n % i == 0) n /= i, c[m]++;
    }
    if (n > 1) p[++m] = n, c[m] = 1;
    return m;
}

质数筛

埃氏筛

小学课本上学过，我们每遍历到一个质数，就把它的倍数划去，最后剩下的未被划去的就是质数。埃氏筛使用 bitset 优化后速度快于线性筛。

bitset<N> isprime;
void solve() {
    isprime.set(), isprime[0] = isprime[1] = false;
    for (int i = 2; i * i <= N; i++) {  // 一个常数优化：筛到 O(sqrt(N)) 即可
        if (!isprime[i]) continue;
        for (int j = 1LL * i * i; j < N; j += i) isprime[j] = false;
    }
}

复杂度为 O(n \log \log n)。但是在实际计算中，bitset 优化后，它比 O(n) 的线性筛更优秀。

线性筛

线性筛虽然跑不过 bitset 埃氏筛，但是它的思想值得学习。

埃氏筛复杂度到不了线性的原因是它会把一个合数划掉两遍。具体地，在标记 i \times prime_j 时，需要确保 i 的最小质因子不小于 p，即 i\bmod prime_j =0 时跳出循环，这样复杂度就降到了 O(n)。

bitset<N> isprime;
int tot, prime[N];
void solve() {
    isprime.set(), isprime[0] = isprime[1] = false;
    for (int i = 2; i < N; i++) {
        if (isprime[i]) prime[++tot] = i;
        for (int j = 1; j <= tot; j++) {
            if (i * prime[j] >= N) break;
            isprime[i * prime[j]] = false;
            if (i & prime[j] == 0) break;
        }
    }
}

区间筛

在一些题目中，我们可能需要求出 [l,r] 中的质数，其中 l,r \le 10^{14} 但 r-l \le 10^7。区间筛可在 O(n \log \log n) 的时间复杂度内解决问题，其中 n=\max(\sqrt{r},r-l)。

观察到 [l,r] 中的合数的最大质因数不会超过 \sqrt{r}，这意味着我们可以用埃氏筛先处理出 [1,\sqrt{r}] 内的质数，再用这些质数去筛掉 [l,r] 中的合数。这也就是埃氏筛那个常数优化的原理。

bitset<N> s, p;
void solve(int a, int b) {  // 筛出[a, b)之间的质数，用 p[i - a] 判断i是否为质数
    s.set(), p.set(), s[0] = s[1] = false;
    if (a <= 1) p[1 - a] = false;
    int z = (int) sqrt(b) + 1;
    for (int i = 2; i < z; i++) {
        for (int j = 2; i * j < z; j++) s[i * j] = false;
    }
    for (int i = 2; i * i <= b; i++) {
        if (!s[i]) continue;
        for (int j = max(i * i, (a + i - 1) / i * i); j < b; j += i) p[j - a] = false;
    }
}

例题

P1147 连续自然数和

由等差数列求和：

\dfrac{(l+r)(r-l+1)}{2}=n

枚举 2n=ij，则

\left\{\begin{matrix} r-l+1=i l+r=j \\ \end{matrix}\right.

解得

\left\{\begin{matrix} l=\dfrac{j-i+1}{2} \\ r=\dfrac{j+i-1}{2} \end{matrix}\right.

显然，i,j 奇偶性应不同，且 i 倒序枚举。复杂度 O(\sqrt{n})。

int n;
void _main() {
    cin >> n; n *= 2;
    for (int i = sqrt(n); i >= 2; i--) {
        if (n % i) continue;
        int j = n / i;
        if ((i & 1) + (j & 1) == 1) cout << (j - i + 1) / 2 << ' ' << (j + i - 1) / 2 << '\n';
    }
}

P1865 A % B Problem

由 m \le 10^6 可以筛一下 m 以内的质数，询问用前缀和处理。远古代码太丑了不想放。

P7960 [NOIP2021] 报数

质数筛的思想延伸。把含有 7 的数先暴力找见，然后把它的倍数筛掉。

const int N = 1e7 + 1000;  // 注意这里要开大点
int t, x, nxt[N];
bitset<N> dis;

inline bool check(int x) {
    for (; x != 0; x /= 10) if (x % 10 == 7) return true;
    return false;
}

void _main() {
    int last = 0;
    dis.reset();
    for (int i = 1; i < N; i++) {
        if (dis[i]) continue;
        if (check(i)) {
            for (int j = 1; j * i < N; j++) dis[j * i] = 1;
            continue;
        }
        nxt[last] = i, last = i;
    }
    for (cin >> t; t--; ) {
        cin >> x;
        cout << (dis[x] ? -1 : nxt[x]) << '\n';
    }
}

P1069 [NOIP 2009 普及组] 细胞分裂

形式化题意：给出 m={m_1}^{m_2} 和 n 个正整数 a_i，求

\min_{i \in [1,n]} \min_{m\mid a_i^k} k

显然这题 m 不能给它算出来，对 m_1 分解质因数，则

m={m_1}^{m_2}=(p_1^{c_1} p_2^{c_2} p_3^{c_3} \cdots p_h^{c_h})^{m_2}=p_1^{c_1m_2} p_2^{c_2m_2} p_3^{c_3m_2} \cdots p_h^{c_hm_2}

这样就得到了 m 的质因数分解。对于每个 a_i，当且仅当 \forall p_j, p_j \mid a_i 时有解，否则 a_i 自乘多少次都不能使 m 成为其约数。

接下来我们求出 p_j 在 a_i 中出现的次数 cnt。则对于 p_j 而言，\min k=\lceil \dfrac{c_j}{cnt} \rceil。然后是一个木桶原理，对所有 \min k 取 \max，为 a_i 的答案。总复杂度 O(n \sqrt{V})，V 是值域。

int n, m1, m2, a[N], tot;

int solve(int x) {
    int res = 0;
    for (int i = 1; i <= tot; i++) {
        if (x % p[i]) return INT_MAX;
        int cnt = 0;
        while (x % p[i] == 0) x /= p[i], cnt++;
        res = max(res, c[i] / cnt + (c[i] % cnt != 0));
    }
    return res;
}

void _main() {
    cin >> n >> m1 >> m2;
    for (int i = 1; i <= n; i++) cin >> a[i];
    tot = decompose(m1);  // 这个函数文章前面有
    for (int i = 1; i <= tot; i++) c[i] *= m2;
    int res = INT_MAX;
    for (int i = 1; i <= n; i++) res = min(res, solve(a[i]));
    cout << (res == INT_MAX ? -1 : res);
}

最大公约数

求解 GCD

STL

C++ STL 中有函数 std::__gcd，复杂度为 O(\log n)，可以直接使用。

分解质因数法

将两数分解质因数：a=p_1^{a_1}p_2^{a_2} \cdots, b=a=p_1^{b_1}p_2^{b_2} \cdots，则 \gcd(a,b)=p_1^{\min(a_1,b_1)}p_2^{\min(a_2,b_2)} \cdots，复杂度为 O(\sqrt{n})。

这证明：求 gcd 和求最小值有共通之处。

辗转相除法

GCD 有如下性质：

\gcd(a,b)=\gcd(b,a \bmod b)

利用此可递归计算 gcd，出口为 \gcd(a,0)=a。

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

更相减损法

GCD 有如下性质：

\gcd(a,b)=\gcd(a-b,b)

如果直接递归计算，复杂度为 O(n)。Stein 算法对此进行了优化：

若 2\mid a 且 2\mid b，则 \gcd(a,b)=2\gcd(\dfrac{a}{2},\dfrac{b}{2})。而若只有 2 \mid a，则 \gcd(a,b)=\gcd(\dfrac{a}{2},b)。这启示我们可以在过程中除掉约数 2，复杂度降到了 O(\log n)。Stein 算法常用于大整数的 GCD。

Stein 算法可以借助二进制内置函数优化，称作 Binary GCD。这份代码的速度比 std::__gcd 更快。可以卡过 P5435。

constexpr int gcd(int a, int b) {
    int az = __builtin_ctz(a), bz = __builtin_ctz(b), z = min(az, bz);
    b >>= bz;
    while (a) {
        a >>= az;
        int d = a - b;
        az = __builtin_ctz(d);
        b = min(a, b), a = abs(d);
    }
    return b << z;
}

ExGCD

ExGCD 用于求解形如 ax+by=\gcd(a,b) 一类的不定方程。

裴蜀定理

对于任意非零整数 a,b，ax+by=c 有整数解当且仅当 c \mid \gcd(a,b)。

不定方程

由 \gcd(a,b)=\gcd(b,a\bmod b) 得：

ax+by=bx'+(a\bmod b)y'\\ =bx'+(a-b\lfloor \dfrac{a}{b} \rfloor)y'\\ =ay'+b(x'-y'\lfloor \dfrac{a}{b} \rfloor)

递归式：x=y',y=x'-y'\lfloor \dfrac{a}{b} \rfloor，这样就可以 O(\log (a+b)) 求解。

void exgcd(int a, int b, int& x, int& y) {
    if (b == 0) return x = 1, y = 0, void();
    exgcd(b, a % b, x, y);
    int t = x; x = y, y = t - a / b * y;
}

乘法逆元

求 x \equiv \dfrac{1}{a} \pmod m 的解。

移项得：ax \equiv 1 \pmod m，即为求解不定方程 ax-bm=1，用 ExGCD 求解即可。

线性方程组

这里的线性方程组是指求一个最小的 x 满足

\left\{\begin{matrix} x \equiv a_1 \pmod {m_1} \\ x \equiv a_2 \pmod {m_2} \\ \cdots \\ x \equiv a_n \pmod {m_n} \end{matrix}\right.

两两考虑。比如我们把前两个方程变成不定方程：x=m_1p+a_1=m_2q+a_2。则 m_1p-m_2q=a_2-a_1，利用 exGCD 解出一组可行解 (p,q)，然后这两个方程组的共同解为 x \equiv (m_1p+a_1) \pmod {\operatorname{lcm}(m_1,m_2)}。然后把 n 个方程都这么合并起来即可。

至于为什么不说 CRT，因为我觉得 exCRT 比 CRT 更容易理解，适用范围更广，代码还好写。

例题

P1072 [NOIP 2009 提高组] Hankson 的趣味题

有一个结论：若 \gcd(a,b)=c，且 a=k_1c,b=k_2c，则 \gcd(k_1,k_2)=1。

证明：考虑反证法，设 K=\gcd(k_1,k_2) \ne 1，则存在不为 1 的正整数 p,q 使得 k_1=pK,k_2=qK，因而 a=pKc,b=qKc，\gcd(a,b)=Kc \ne c，故原结论成立。

由题可得：x=a_1 p，b_1=xq，其中 p,q 为正整数。

由上述结论得：\gcd(\dfrac{a_0}{a_1},\dfrac{x}{a_1})=1。又有 \operatorname{lcm}(x,b_0)=b_1，故 \gcd(\dfrac{b_1}{x},\dfrac{b_1}{b_0})=1。且 x\mid b_1，枚举 x 并判断即可，复杂度 O(\sqrt{b_1} \log V)。

int a0, a1, b0, b1;

void _main() {
    cin >> a0 >> a1 >> b0 >> b1;
    int cnt = 0;
    for (int x = 1; x * x <= b1; x++) {
        if (b1 % x) continue;
        if (x % a1 == 0 && __gcd(x / a1, a0 / a1) == 1 && __gcd(b1 / x, b1 / b0) == 1) cnt++;
        int y = b1 / x;
        if (x == y) continue;
        if (y % a1 == 0 && __gcd(y / a1, a0 / a1) == 1 && __gcd(b1 / y, b1 / b0) == 1) cnt++;
    }
    cout << cnt << '\n';
}

P4549 【模板】裴蜀定理

扩展裴蜀定理。由 \gcd(a,b,c)=\gcd(a,\gcd(b, c)) 和裴蜀定理得

\sum_{i=1}^{n} a_i x_i = \gcd_{i=1}^{n} |a_i|

而题目中记 \sum x_i 为 S，可知 \gcd_{i=1}^{n} |a_i| 是 S 的约数，S 最小时二者相等。也就是求所有数的 \gcd。

int n, res, a;
int main() {   // 两年前的马蜂
    scanf("%d", &n);
    while (n--) scanf("%d", &a), a = abs(a), res = __gcd(res, a);
    printf("%d", res);
    return 0;
}

P1082 [NOIP 2012 提高组] 同余方程

这东西我们在之前的乘法逆元讲过方法了，就是用 exGCD 解方程 ax-by=1。输入保证有解，意味着 \gcd(a,b)=1，所以这就是 exGCD 方程的标准形式。细节上，求出 x 以后要处理负数问题。

template <class T> void exgcd(T a, T b, T& x, T& y) {
    if (b == 0) return x = 1, y = 0, void();
    exgcd(b, a % b, y, x), y -= a / b * x;
}

long long a, b;
void _main() {
    cin >> a >> b;
    long long x, y; exgcd(a, b, x, y);
    cout << (x % b + b) % b;
}

P5656 【模板】二元一次不定方程 (exgcd)

首先由裴蜀定理判无解，若 c 不是 \gcd(a,b) 的倍数则直接无解。

设 \gcd(a,b)=d。考虑用 exGCD 求解方程 ax+by=d，将解记作 x_0,y_0。则

ax_0+by_0=d\\

方程两边同乘 \dfrac{c}{d} 转化为所求：

a\dfrac{cx_0}{d}+b\dfrac{cy_0}{d}=c\\

故原方程的一组解为 x_1=\dfrac{cx_0}{d},y_1=\dfrac{cy_0}{d}。

接下来考虑构造通解形式，设

a(x_1+m)+b(y_0+n)=c

不难发现 m,n 满足条件 am+bn=0。仍然利用裴蜀定理，解得 m=p \times \dfrac{b}{d}, n=-p \times {a}{d}，其中 p 是正整数。于是我们得到原方程的通解

\left\{\begin{matrix} x=x_1+p \times \dfrac{b}{d} \\ y=y_1-p \times \dfrac{a}{d} \end{matrix}\right.

考虑求解数与最值。找 x_{\min} 即为解关于 k 的不等式

x_1+km \ge 1

由于式中的值均为整数，故

k \ge \lceil \dfrac{1-x_1}{m} \rceil

由于 y 随 x 的增大而减小，故 x_{\min} 所对应的 y 就是 y_{\max}。接下来我们用 y_{\max} 推 y_{\min}。推导一下，可以发现

y_{\min}=y_{\max} \bmod n

于是 x_{\max} 加上同样多的 m 即可：

x_{\max}=x_{\min} + \lfloor \dfrac{y_{\max}-1}{n} \times m \rfloor

解数就是 [y_{\min},y_{\max}] 中解的个数，即

\lfloor \dfrac{y_{\max}-1}{n} +1 \rfloor

至此本题在 O(\log V) 内解决，V 为值域。

在实现中，我们可以将 a,b,c 都除去 d，这样可以简化代码。

template <class T> void exgcd(T a, T b, T& x, T& y) {
    if (b == 0) return x = 1, y = 0, void();
    exgcd(b, a % b, y, x), y -= a / b * x;
}

long long a, b, c;

void _main() {
    cin >> a >> b >> c;
    long long d = __gcd(a, b);
    if (c % d) return cout << -1 << '\n', void();
    a /= d, b /= d, c /= d;
    long long x0, y0; exgcd(a, b, x0, y0);
    long long x1 = c * x0, y1 = c * y0;
    int xmin = (x1 > 0 && x1 % b != 0) ? x1 % b : x1 % b + b;
    int ymax = (c - xmin * a) / b;
    int ymin = (y1 > 0 && y1 % a != 0) ? y1 % a : y1 % a + a;
    int xmax = (c - ymin * b) / a;
    if (xmax <= 0) cout << xmin << ' ' << ymin << '\n';
    else cout << (ymax - ymin) / a + 1 << ' ' << xmin << ' ' << ymin << ' ' << xmax << ' ' << ymax << ' ' << '\n';
}

P1516 青蛙的约会

设相遇时两只青蛙跳了 t 次，则

(n-m)t+kL=x-y

即求解 t 的最小非负整数解。在不定方程 ax+by=c 中，我们把 n-m 视作 a，L 视作 b，x-y 视作 c，则应用上一题的方法来求解。

template <class T> void exgcd(T a, T b, T& x, T& y) {
    if (b == 0) return x = 1, y = 0, void();
    exgcd(b, a % b, y, x), y -= a / b * x;
}
long long x, y, m, n, L;

void _main() {
    cin >> x >> y >> m >> n >> L;
    long long a = n - m, b = L, c = x - y;
    if (a < 0) a = -a, c = -c;
    long long d = __gcd(a, b);
    if (c % d) return cout << "Impossible", void();
    a /= d, b /= d, c /= d;
    long long x0, y0; exgcd(a, b, x0, y0);
    long long x1 = c * x0;
    cout << ((x1 > 0 && x1 % b != 0) ? x1 % b : x1 % b + b);
}

P4777 【模板】扩展中国剩余定理（EXCRT）

解线性方程组的方法一般称作 CRT。套用上面的 exGCD 法即可。

const int N = 1e5 + 5;

template <class T> void exgcd(T a, T b, T& x, T& y) {
    if (!b) return x = 1, y = 0, void();
    exgcd(b, a % b, y, x), y -= a / b * x;
}
template <class T> T crt(int n, const T* a, const T* m) {
    T x = 0, y = 0, p = m[1], res = a[1];
    for (int i = 2; i <= n; i++) {
        T a0 = p, b0 = m[i], c = (a[i] - res % b0 + b0) % b0, g = __gcd(a0, b0);
        if (c % g) return -1;
        exgcd(a0, b0, x, y);
        x = (__int128) x * c / g % b0, res += x * p, p = p / __gcd(p, b0) * b0, res = (res % p + p) % p;
    }
    return res;
}

int n;
long long a[N], m[N];

void _main() {
    cin >> n;
    for (int i = 1; i <= n; i++) cin >> m[i] >> a[i];
    cout << crt(n, a, m);
}

P3868 [TJOI2009] 猜数字

对于 (x-a_i) \mid b_i，可以变形为 x-a_i \equiv 0 \pmod {b_i}，然后再移项就是 x \equiv a_i \pmod {b_i}。于是解这个同余方程组即可。

细节上注意 a_i 可能为负，需要对 b_i 取模处理，另外这题会爆 long long，不过上面板子里已经开了 __int128 了。

const int N = 15;
int n;
long long a[N], b[N];

void _main() {
    cin >> n;
    for (int i = 1; i <= n; i++) cin >> a[i];
    for (int i = 1; i <= n; i++) cin >> b[i];
    for (int i = 1; i <= n; i++) a[i] = (a[i] % b[i] + b[i]) % b[i];
    cout << crt(n, a, b);
}

欧拉函数

定义

欧拉函数 \varphi(x) 指 [1,n] 中与 n 互质的数的个数。

若由算术基本定理可将 n 分解为 p_1^{c_1} p_2^{c_2} p_3^{c_3} \cdots p_m^{c_m}，则

\varphi(x)=n (1-\dfrac{1}{p_1}) (1-\dfrac{1}{p_2}) (1-\dfrac{1}{p_3}) \cdots (1-\dfrac{1}{p_m}) =n\times \prod_{i=1}^{m} (1-\dfrac{1}{p_i})

由此有一个 O(\sqrt{n}) 计算欧拉函数的方法：

inline int phi(int n) {
    if (n == 1) return 1;
    int res = n;
    for (int 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;
}

性质

性质 1：若 p 是质数，则 \varphi(p)=p-1。

性质 2：若 a,b 互质，则 \varphi(ab)=\varphi(a) \times \varphi(b)。这表明：欧拉函数是积性函数，因而可以用线性筛筛出。

性质 3：\sum_{d \mid n} \varphi(d)=n。

性质 4：对于 n>1，[1,n] 中与 n 互质的数的和为 \dfrac{1}{2}n\times \varphi(n)。

性质 5：若 n=p^k，其中 p 是质数，则 \varphi(n)=p^k-p^{k-1}。

筛法

使用性质 2 并结合线性筛法，可以在 O(n) 时间内预处理 [1,n] 的欧拉函数值：

int len, phi[N], prime[N];
bitset<N> isprime;

inline void eular(int n) {
    phi[1] = 1, isprime.set(), isprime[0] = isprime[1] = false;
    for (int i = 2; i <= n; i++) {
        if (isprime[i]) prime[++len] = i, phi[i] = i - 1;
        for (int j = 1; j <= len && i * prime[j] <= n; j++) {
            isprime[i * prime[j]] = false;
            if (i % prime[j] == 0) {
                phi[i * prime[j]] = phi[i] * prime[j];
                break;
            }
            phi[i * prime[j]] = phi[i] * phi[prime[j]];
        }
    }
}

例题

SP4141 ETF - Euler Totient Function

线性筛欧拉函数板子。直接上代码：

const int N = 1e6 + 5;
void _main() {
    eular(N - 1);
    int t, x; cin >> t;
    while (t--) cin >> x, cout << phi[x] << '\n';
}

UVA10179 Irreducable Basic Fractions

题意翻译：求 \dfrac{0}{n},\dfrac{1}{n},\dfrac{2}{n},\cdots \dfrac{n-1}{n} 中多少个分数为最简分数。认为 \dfrac{0}{n} 最简。

我们发现最简分数 \dfrac{x}{n} 满足 \gcd(x,n)=1，于是所求为 \varphi(n)。用分解质因数法求即可，代码不放了。

P2158 [SDOI2008] 仪仗队

由样例图可以发现能看见的位置关于对角线对称，只需要考虑一个三角形即可。以左下角为原点 (0,0) 建立坐标系，则一个点被阻挡的条件就是到原点连线的斜率相同。设两点 (x_1,y_1),(x_2,y_2)，则 \dfrac{y_1}{x_1}=\dfrac{y_2}{x_2}。当且仅当这个分数已经为最简形式时它不会被挡住，这就和上题类似了。

所求即为

2\sum_{i=1}^{n-1} \varphi(i)+1

这是因为 (2,2) 满足条件而我们无法统计进去。注意特判 n=1。

用线性筛筛出欧拉函数即可，复杂度 O(n)，代码：

int n;
void _main() {
    cin >> n;
    if (n == 1) return cout << 0, 0;
    int res = 1;
    eular(40000);
    for (int i = 1; i < n; i++) res += 2 * phi[i];
    cout << res;
}

P2398 GCD SUM

这题做法很多，这里写一个欧拉函数做法。

可以发现 \gcd(i,j) 的值只有 n 种，不妨考虑 gcd(i,j)=k 的数目。对于一对互质的 i,j，有 \gcd(ik,jk)=k，所以可以得到结果为 k 的数目为

2\sum_{i=1}^{\lfloor \frac{n}{k} \rfloor} \varphi(i)-1

和上面那题是一样的，因为 \gcd(i,j)=\gcd(j,i) 有对称性，而 (1,1) 会重复计算，要减一。

然后我们用线性筛处理出 \varphi(i) 的前缀和，O(n) 计算即可。

const int N = 1e5 + 5;
int n;
long long pre[N];

void _main() {
    cin >> n;
    eular(n);
    for (int i = 1; i <= n; i++) pre[i] = pre[i - 1] + phi[i];
    long long res = 0;
    for (int i = 1; i <= n; i++) res += 1LL * i * (2 * pre[n / i] - 1);
    cout << res;
}

双倍经验：P1390。这个题没有对称性，把系数 2 去掉即可。

常用定理

费马小定理

若 p 为质数，且 \gcd(a,p)=1，则 a^{p-1} \equiv 1 \pmod p。

变形式：a^p\equiv a \pmod p 和 a^{p-2}\equiv \dfrac{1}{a} \pmod p。其中第二个式子表明，当 p 为质数时可以直接用快速幂求逆元。

证明：取一个不为 p 的倍数的数 a，构造序列 A=\{1,2,3, \cdots p-1\}，则：

\prod_{i=1}^{p-1} A_i \equiv \prod_{i=1}^{p-1} (A_i\times a) \pmod p

考虑每一个 A_i 都不是 p 的约数易证。则令 m=(p-1)!，则

m \equiv \prod_{i=1}^{p-1} (A_i\times a) \pmod p

又有 a^{p-1} \times f \equiv f \pmod p，故

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

欧拉定理

欧拉定理是费马小定理的扩展形式。若 \gcd(a,n)=1，则 a^{\varphi(n)} \equiv 1 \pmod n。证明与费马小定理类似，构造一个与 n 互质的序列即可。

推论 1：若 \gcd(a,n)=1，则 a^b \equiv a^{b \bmod \varphi(n)} \pmod n。

推论 2：若 \gcd(a,n)=1，则满足 a^x \equiv 1 \pmod n 的最小正整数是 \varphi(n) 的约数。

扩展欧拉定理

扩展欧拉定理在 OI 中常用于对幂取模的情况。

a^b \equiv \left\{\begin{matrix} a^{b \bmod \varphi(n)} & \gcd(a,n)=1 \\ a^b & \gcd(a,n) \ne 1,b < \varphi(n) \\ a^{(b \bmod \varphi(n))+\varphi(n)} & \gcd(a,n) \ne 1, b \ge \varphi(n) \end{matrix}\right. \pmod n

例题

P4139 上帝与集合的正确用法

题里给的那个东西实际上是 2^{2^{2^{2^ \cdots}}} \bmod p。对指数塔不断递归，用扩展欧拉定理降幂即可。\varphi(n) 提前筛出来。

由于指数塔是无限层，因此 b \ge \varphi(n)，直接认为是第三种情况就行了。

long long p;
long long power(long long a, long long b, long long p) {
    long long res = 1; for (a %= p; b; a = a * a % p, b >>= 1) {
        if (b & 1) res = res * a % p;
    } return res;
}
long long f(long long p) {return p == 1 ? 0 : power(2, f(phi[p]) + phi[p], p);}

void _main() {   // 这里已经筛过phi[n]了
    cin >> p; cout << f(p) << '\n';
}

P10496 The Luckiest Number

设答案为 x 个 8 连在一起的整数 n，则

\begin{aligned} n&=8+8\times 10+8\times 10^2 + \cdots + 8 \times 10^x\\ &=\sum_{i=0}^{x} 8 \times 10^i \\ &=\dfrac{8(10^x-1)}{9} \end{aligned}

设 d=\gcd(L,8)，则由题知 L \mid \dfrac{8(10^x-1)}{9}，故 9L \mid 8(10^x-1)，即

\dfrac{9L}{d} \mid (10^x-1)

因此，10^ x \equiv 1 \pmod {\dfrac{9L}{d}}。

由欧拉定理的推论 2，我们求出 \varphi(\dfrac{9L}{d}) 并枚举其约数检查即可，复杂度 O(\sqrt{L} \log L)。

long long power(long long a, long long b, long long p) {
    long long res = 1; for (a %= p; b; a = (__int128) a * a % p, b >>= 1) {
        if (b & 1) res = (__int128) res * a % p;
    } return res;
}
vector<long long> factors(long long n) {
    vector<long long> res;
    for (long long i = 1; i * i <= n; i++) {
        if (n % i) continue;
        res.emplace_back(i);
        if (n / i != i) res.emplace_back(n / i);
    }
    sort(res.begin(), res.end());
    return res;
}

long long n;

void _main() {
    for (int kase = 1; ; kase++) {
        cin >> n;
        if (n == 0) break;
        cout << "Case " << kase << ": ";
        n = 9 * n / __gcd(n, 8LL);
        if (__gcd(n, 10LL) != 1) {cout << "0\n"; continue;}
        vector<long long> f = factors(phi(n));
        for (long long x : f) {
            if (power(10, x, n) == 1) {cout << x << '\n'; break;}
        }
    }
}

数论分块

数论分块用于快速计算形如

\sum_{i=1}^{n} f(i) g(\lfloor \dfrac{k}{i} \rfloor)

的式子，其中 f(i) 可处理出前缀和或可以快速计算 f(x)-f(y)。如果这是 O(1) 的，则数论分块可在 O(\sqrt{n}) 的时间内得出结果。

原理

以函数 y=\lfloor \dfrac{100}{x} \rfloor 为例，图像如下：

我们注意到对于每个固定的 y，x 的取值总是一个固定区间。事实上，\lfloor \dfrac{k}{i} \rfloor 不变时，i \le \lfloor \frac{k}{\lfloor \frac{k}{i} \rfloor} \rfloor。

证明：设 m=\lfloor \dfrac{k}{i} \rfloor，则 m \le \dfrac{k}{i}，所以 \lfloor \dfrac{k}{m} \rfloor \ge \lfloor \frac{k}{\frac{k}{i}} \rfloor=i。

因此，i_{\max}=\lfloor \dfrac{k}{m} \rfloor=\lfloor \frac{k}{\lfloor \frac{k}{i} \rfloor} \rfloor。

板子

template <class F_t, class G_t>
long long sqrt_decomposition(long long n, long long k, F_t f, G_t g) {
    long long res = 0;
    for (long long l = 1, r = 0; l <= n; l = r + 1) {
        r = (k / l) ? min(n, (k / (k / l))) : n;
        res += f(r, l - 1) * g(k / l);
    }
    return res;
}

这个板子第一个参数传入 n，第二个参数传入 k，第三个参数传入一个函数计算 f(x)-f(y)，第四个参数传入 g(x) 的函数。

例题

UVA11526 H(n)

【模板】数论分块。

int n;
void _main() {
    cin >> n;
    cout << sqrt_decomposition(n, n,
        [](long long x, long long y) {return x - y;}, // f(x) = 1
        [](long long x) {return x;} // g(x) = x
    ) << '\n';
}

P2261 [CQOI2007] 余数求和

首先有一个结论是 a \bmod b=a-b \times \lfloor \dfrac{a}{b} \rfloor，然后推一波式子：

\begin{aligned} G(n, k) &= \sum_{i = 1}^n k \bmod i \\ &= \sum_{i = 1}^n k -i \times \lfloor \dfrac{k}{b} \rfloor \\ &= nk-\sum_{i = 1}^n i \times \lfloor \dfrac{k}{b} \rfloor \end{aligned}

后面这个东西可以数论分块来做。

long long n, k;

void _main() {
    cin >> n >> k;
    cout << n * k - sqrt_decomposition(n, k, 
        [](long long x, long long y) {return (x - y) * (x + y + 1) / 2;}, // f(x) = x
        [](long long x) {return x;} // g(x) = x
    );
}

P2260 [清华集训 2012] 模积和

和上面的题很像。

推柿子：

\begin{aligned} ans &=\sum_{i=1}^{n} \sum_{j=1}^{m} (n \bmod i) \times (m \bmod j), i \neq j \\ &=\sum_{i=1}^{n} \sum_{j=1}^{m} (n \bmod i) \times (m \bmod j)-\sum_{i=1}^{\min(n,m)} (n \bmod i) \times (m \bmod i)\\ &=\sum_{i=1}^n(n-\lfloor\frac{n}{i}\rfloor\times i)\times\sum_{j=1}^m(m-\lfloor\frac{m}{j}\rfloor\times j)-\sum_{i=1}^{\min(n,m)}(n-\lfloor\frac{n}{i}\rfloor\times i)\times(m-\lfloor\frac{m}{i}\rfloor\times i)\\ &=(n^2-\sum_{i=1}^ni\times\lfloor\frac{n}{i}\rfloor)\times(m^2-\sum_{i=1}^mi\times\lfloor\frac{m}{i}\rfloor)-\sum_{i=1}^{\min(n,m)}(nm-mi\times\lfloor\frac{n}{i}\rfloor-ni\times\lfloor\frac{m}{i}\rfloor+i^2\times\lfloor\frac{n}{i}\rfloor\times\lfloor\frac{m}{i}\rfloor)\\ &=(n^2-\sum_{i=1}^ni\times\lfloor\frac{n}{i}\rfloor)\times(m^2-\sum_{i=1}^mi\times\lfloor\frac{m}{i}\rfloor)-nm\times \min(n,m)+\sum_{i=1}^{\min(n,m)}mi\times\lfloor\frac{n}{i}\rfloor+\sum_{i=1}^{\min(n,m)} ni\times\lfloor\frac{m}{i}\rfloor -\sum_{i=1}^{\min(n,m)} i^2\times\lfloor\frac{n}{i}\rfloor\times\lfloor\frac{m}{i}\rfloor \end{aligned}

最后这个 i^2\times\lfloor\frac{n}{i}\rfloor\times\lfloor\frac{m}{i}\rfloor 用数论分块的时候不太一样。首先我们需要保证块内的 \lfloor\frac{n}{i}\rfloor\times\lfloor\frac{m}{i}\rfloor 相同，将右端点改为 \min(\lfloor \frac{n}{\lfloor \frac{n}{i} \rfloor} \rfloor,\lfloor \frac{m}{\lfloor \frac{m}{i} \rfloor} \rfloor)。然后用到一个结论 \sum_{i=1}^{n} i^2=\frac{n(n+1)(2n+1)}{6}。

另外本题一堆取模，笔者用了自己写的 modint。

long long n, m;
modint pre(modint n) {return n * (n + 1) * (n * 2 + 1) / modint(6);}
modint k1() {  // 这里没法套板子，自己写一个
    modint res = 0;
    for (long long l = 1, r = 0; l <= min(n, m); l = r + 1) {
        r = min(n / (n / l), m / (m / l));
        res += modint(n / l) * modint(m / l) * (pre(r) - pre(l - 1));
    }
    return res;
}

void _main() {
    cin >> n >> m; long long k = min(n, m);
    modint n0 = n, m0 = m;
    modint n1 = n0 * n0 - sqrt_decomposition(n, n,
        [](modint x, modint y) {return (x - y) * (x + y + 1) / 2;}, 
        [](long long x) {return x;} 
    );
    modint m1 = m0 * m0 - sqrt_decomposition(m, m,
        [](modint x, modint y) {return (x - y) * (x + y + 1) / 2;}, 
        [](long long x) {return x;}
    );
    modint n2 = sqrt_decomposition(k, n, 
        [](modint x, modint y) {return (x - y) * (x + y + 1) / 2;}, 
        [&](long long x) {return m0 * x;}
    );
    modint m2 = sqrt_decomposition(k, m, 
        [](modint x, modint y) {return (x - y) * (x + y + 1) / 2;}, 
        [&](long long x) {return n0 * x;}
    );
    cout << n1 * m1 - n0 * m0 * k + n2 + m2 - k1();
}

P2424 约数和

首先求 f(x) 的前缀和 g(x)。有一个结论是 i 的约数在 n 以内有 \lfloor \dfrac{n}{i} \rfloor 个，所以

g(x) = \sum_{i=1}^{n} i \times \lfloor \dfrac{n}{i} \rfloor

然后套板子去求即可。

long long g(long long n) {
    if (n <= 1) return n;
    return sqrt_decomposition(n, n, 
        [](long long x, long long y) {return (x - y) * (x + y + 1) / 2;},
        [](long long x) {return x;}
    );
}

long long l, r;
void _main() {
    cin >> l >> r;
    cout << g(r) - g(l - 1);
}