Acwing - 算法基础课 - 笔记（数学知识

文章目录

数学知识（二）
- 欧拉函数
- - 公式法
  - 筛法
  - 欧拉定理
- 快速幂
- 扩展欧几里得算法
- 中国剩余定理

数学知识（二）

这一小节主要讲解的内容是：欧拉函数，快速幂，扩展欧几里得算法，中国剩余定理。

这一节内容偏重于数学推导，做好心理准备。

欧拉函数

公式法

什么是欧拉函数呢？

欧拉函数用 ϕ ( n ) \phi(n) ϕ(n) 来表示，它的含义是， 1 1 1 到 n n n 中与 n n n 互质的数的个数

比如， ϕ ( 6 ) = 2 \phi(6) = 2 ϕ(6)=2，解释：1到6当中，与6互质的数只有1，5，共两个数。

两个数 a a a , b b b 互质的含义是 g c d ( a , b ) = 1 gcd(a,b) = 1 gcd(a,b)=1

欧拉函数如何求解呢？

对于一个数 N N N ，将其写为分解质因数的形式 N = P 1 k 1 × P 2 k 2 × . . . × P n k n N = P_1^{k_1} \times P_2^{k_2} \times ... \times P_n^{k_n} N=P1k1×P2k2×...×Pnkn

则 ϕ ( N ) = N × ( 1 − 1 P 1 ) × ( 1 − 1 P 2 ) × . . . × ( 1 − 1 P n ) \phi(N) = N \times (1-\frac{1}{P_1}) \times (1-\frac{1}{P_2}) \times ... \times (1-\frac{1}{P_n}) ϕ(N)=N×(1−P11)×(1−P21)×...×(1−Pn1)，这就是欧拉函数的求解公式

比如 N = 6 N = 6 N=6， 6有2个质因子2和3，则

ϕ ( 6 ) = 6 × ( 1 − 1 2 ) × ( 1 − 1 3 ) = 2 \phi(6) = 6 \times (1-\frac{1}{2}) \times (1-\frac{1}{3}) = 2 ϕ(6)=6×(1−21)×(1−31)=2

欧拉函数的公式如何证明呢？

欧拉函数的证明需要用到容斥原理，容斥原理的具体内容将在下一节进行讲解，本节直接借用容斥原理的结论。

注：下面的 N P \frac{N}{P} PN 实际指的都是 ⌊ N P ⌋ \lfloor\frac{N}{P}\rfloor ⌊PN⌋，比如 15 7 \frac{15}{7} 715 指的是 ⌊ 15 7 ⌋ = 2 \lfloor\frac{15}{7}\rfloor = 2 ⌊715⌋=2，由于加上取整符号会导致编写latext语法变得非常繁琐，特此申明。

我们求 1 1 1 到 n n n 中，与 n n n 互质的数的个数，仍然采用删数的思想，只要把那些和 n n n 有公约数的数删掉就可以了。那么还是需要先对 n n n 做分解质因数，因为 n n n 的质因数是 n n n 的约数的最小组成部分。

从 1 1 1 到 N N N 中，去掉其全部质因子 P 1 ， P 2 , . . . , P n P_1，P_2,...,P_n P1，P2,...,Pn 的所有倍数。这些数都和 N N N 不互质，因为都有公因子 P i P_i Pi

那么剩余的数的个数就是 N − N P 1 − N P 2 − . . . − N P n N - \frac{N}{P_1} - \frac{N}{P_2} - ... - \frac{N}{P_n} N−P1N−P2N−...−PnN
加上所有 P i × P j P_i \times P_j Pi×Pj 的倍数的个数，其中 i i i 和 j j j 是从 1 1 1 到 n n n 之中任选两个数

假设第一步计算得到的结果是 S 1 S_1 S1，那么第二步要做的就是 S 1 + N P 1 × P 2 + N P 1 × P 3 + . . . N P 1 × P n + N P 2 × P 3 + . . . + N P 2 × P n + . . . + N P n − 1 × P n S_1 + \frac{N}{P_1 \times P_2} + \frac{N}{P_1 \times P_3} + ...\frac{N}{P_1 \times P_n} + \frac{N}{P_2 \times P_3} + ... + \frac{N}{P_2 \times P_n} + ... + \frac{N}{P_{n-1} \times P_n} S1+P1×P2N+P1×P3N+...P1×PnN+P2×P3N+...+P2×PnN+...+Pn−1×PnN

因为在第一步时，会多删掉一些数，比如某个数既是 P 1 P_1 P1 的倍数，也是 P 2 P_2 P2 的倍数（即这个数是 P 1 × P 2 P_1 \times P_2 P1×P2 的倍数），那么这个数会被 P 1 P_1 P1 减一次，被 P 2 P_2 P2 减一次，一共减了两次，多减了一次，那么我们需要加一次加回来。
加上所有 P i × P j × P k P_i \times P_j \times P_k Pi×Pj×Pk 的倍数的个数

若有某个数，同时是3个质因子的倍数，假设一个数同时是 P 1 P_1 P1， P 2 P_2 P2， P 3 P_3 P3 的倍数，那么这个数第一步会被减3次，会在第二步被加3次，实际是没减的。所以还要对这种数减去一次

假设第二步计算得到的结果是 S 2 S_2 S2，那么第三步要做的就是 S 2 − N P 1 × P 2 × P 3 − N P 1 × P 2 × P 4 − . . . − N P n − 2 × P n − 1 × P n S_2 - \frac{N}{P_1 \times P_2 \times P_3} - \frac{N}{P_1 \times P_2 \times P_4} - ... - \frac{N}{P_{n-2} \times P_{n-1} \times P_n} S2−P1×P2×P3N−P1×P2×P4N−...−Pn−2×Pn−1×PnN

即，对 i , j , k ∈ [ 1 , n ] i,j,k \in [1,n] i,j,k∈[1,n] 的全部组合，进行减
后续的步骤同理

若若某个数，同时是4个质因子的倍数，这种数在第一步被减了4次，第二步被加了 C 4 2 = 4 × 3 2 × 1 = 6 C_4^2 = \frac{4 \times 3}{2 \times 1} = 6 C42=2×14×3=6 次，第三步被减了 C 4 3 = 4 C_4^3 = 4 C43=4 次，总共是减了2次，则需要加回一次

所以要做的就是 S 3 + N P 1 × P 2 × P 3 × P 4 + . . . . + N P n − 3 × P n − 2 × P n − 1 × P n S_3 + \frac{N}{P_1 \times P_2 \times P_3 \times P_4} + .... + \frac{N}{P_{n-3} \times P_{n-2} \times P_{n-1} \times P_n} S3+P1×P2×P3×P4N+....+Pn−3×Pn−2×Pn−1×PnN
对于那些是5个质因子的倍数的数，第一次被减了 C 5 1 = 5 C_5^1 = 5 C51=5 次，第二步加了 C 5 2 = 10 C_5^2 = 10 C52=10 次，第三步被减了 C 5 3 = 10 C_5^3 = 10 C53=10 次，第四步被加了 C 5 4 = 5 C_5^4 = 5 C54=5 次，总共是0次操作，所以需要减一次
…

是不是感觉像是无限套娃？！hhhhh

这个过程会在第 n n n 次后停止，而 n n n 是 N N N 的质因数个数

所以， 1 1 1 到 N N N 中所有和 N N N 互质的数的个数，就等于

N − N P 1 − N P 2 − . . . − N P n N - \frac{N}{P_1} - \frac{N}{P_2} - ... - \frac{N}{P_n} N−P1N−P2N−...−PnN

+ N P 1 × P 2 + N P 1 × P 3 + . . . N P 1 × P n + N P 2 × P 3 + . . . + N P 2 × P n + . . . + N P n − 1 × P n + \frac{N}{P_1 \times P_2} + \frac{N}{P_1 \times P_3} + ...\frac{N}{P_1 \times P_n} + \frac{N}{P_2 \times P_3} + ... + \frac{N}{P_2 \times P_n} + ... + \frac{N}{P_{n-1} \times P_n} +P1×P2N+P1×P3N+...P1×PnN+P2×P3N+...+P2×PnN+...+Pn−1×PnN

− N P 1 × P 2 × P 3 − N P 1 × P 2 × P 4 − . . . − N P n − 2 × P n − 1 × P n - \frac{N}{P_1 \times P_2 \times P_3} - \frac{N}{P_1 \times P_2 \times P_4} - ... - \frac{N}{P_{n-2} \times P_{n-1} \times P_n} −P1×P2×P3N−P1×P2×P4N−...−Pn−2×Pn−1×PnN

+ N P 1 × P 2 × P 3 × P 4 + . . . . + N P n − 3 × P n − 2 × P n − 1 × P n + \frac{N}{P_1 \times P_2 \times P_3 \times P_4} + .... + \frac{N}{P_{n-3} \times P_{n-2} \times P_{n-1} \times P_n} +P1×P2×P3×P4N+....+Pn−3×Pn−2×Pn−1×PnN

. . . . .... ....

± N P 1 × P 2 × . . . × P n \pm \frac{N}{P_1 \times P_2 \times ... \times P_n} ±P1×P2×...×PnN

而这一大坨式子，化简后得到的就是

N × ( 1 − 1 P 1 ) × ( 1 − 1 P 2 ) × . . . × ( 1 − 1 P n ) N \times (1-\frac{1}{P_1}) \times (1-\frac{1}{P_2}) \times ... \times (1-\frac{1}{P_n}) N×(1−P11)×(1−P21)×...×(1−Pn1)

这也就是本节开头提到的欧拉函数的公式

练习题：acwing - 873: 欧拉函数

#include <iostream>int euler(int x) {int ans = x;for(int i = 2; i <= x / i; i++) {if(x % i == 0) {while(x % i == 0) x /= i;ans = ans / i * (i - 1);}}if(x > 1) ans = ans / x * (x - 1);return ans;
}int main() {int n, a;scanf("%d" ,&n);while(n--) {scanf("%d", &a);printf("%d\n", euler(a));}
}

注意，欧拉公式的计算过程中，每一项计算完毕后都应该是一个整数，不应该出现小数。

故对于 N × ( 1 − 1 P 1 ) N \times (1-\frac{1}{P_1}) N×(1−P11) 这样的式子，计算时我们采用 N P 1 × ( P 1 − 1 ) \frac{N}{P_1} \times (P_1-1) P1N×(P1−1)，这样能保证得到的结果是整数，而不出现小数。

由于 P 1 P_1 P1 是 N N N 的质因子，故 N P 1 \frac{N}{P_1} P1N 一定是整数，所以 N P 1 × ( P 1 − 1 ) \frac{N}{P_1} \times (P_1-1) P1N×(P1−1) 就一定是整数，令 S = N P 1 × ( P 1 − 1 ) S = \frac{N}{P_1} \times (P_1-1) S=P1N×(P1−1)，那么在第二轮迭代中，我们需要计算 S × ( 1 − 1 P 2 ) S \times (1-\frac{1}{P_2}) S×(1−P21)，同样的，我们计算 S P 2 × ( P 2 − 1 ) \frac{S}{P_2} \times (P_2-1) P2S×(P2−1)，那么 S P 2 \frac{S}{P_2} P2S 一定是整数吗？答案是肯定的。

因为 N = P 1 k 1 × P 2 k 2 × . . . × P n k n N = P_1^{k_1} \times P_2^{k_2} \times ... \times P_n^{k_n} N=P1k1×P2k2×...×Pnkn，容易得到，第一步计算的结果 S = P 1 k 1 − 1 × P 2 k 2 × . . . × P n k n × ( P 1 − 1 ) S=P_1^{k_1-1} \times P_2^{k_2} \times ... \times P_n^{k_n} \times (P_1-1) S=P1k1−1×P2k2×...×Pnkn×(P1−1)，所以 S S S 是一定能被 P 2 P_2 P2 整除的，同理，在第三轮迭代时，第二轮的计算结果仍然能被 P 3 P_3 P3 整除，所以我们就可以在每轮迭代时，计算 S i P i × ( P i − 1 ) \frac{S_i}{P_i} \times (P_i-1) PiSi×(Pi−1)，这样能保证整个计算的过程中，不出现小数，即可保证最终结果的正确性。

筛法

上面的公式法，适用于求解某一个数的欧拉函数，就类似于用试除法判断某一个数是否是质数。

然而，有的时候，我们需要求解某一个范围内全部数的欧拉函数（比如求解 1 1 1 到 N N N 之间所有数的欧拉函数），此时若对每个数依次套用欧拉公式，则整体的时间复杂度为 O ( N × N ) O(N \times \sqrt{N}) O(N×N ) ，因为欧拉函数的计算依赖于分解质因数，而分解质因数的时间复杂度是 O ( N ) O(\sqrt{N}) O(N ) 。这个时间复杂度是不被接受的，太慢了，所以我们需要变换思路。

联想到在求解 1 1 1 到 N N N 之间全部质数的时候，我们采用的是筛法，而不是对每个数依次判断是否是质数。类似的，求解 1 1 1 到 N N N 之间全部数的欧拉函数，我们也可以用类似的思想。

借鉴前面筛选质数时所采用的线性筛法，能够在 O ( N ) O(N) O(N) 的时间复杂度内求解出 1 1 1 到 N N N 每个数的欧拉函数。在本章（数论）后面的学习中，会发现线性筛法在执行过程中不仅仅能求出欧拉函数，还能求出很多其他的内容。

我们先把线性筛法的代码写一遍

#include<iostream>
#include<algorithm>
using namespace std;const int N = 1e6 + 10;bool st[N];int primes[N], ctn;void get_primes_linear(int n) {for(int i = 2; i <= n; i++) {if(!st[i]) primes[ctn++] = i;for(int j = 0; primes[j] <= n / i; j++) {st[primes[j] * i] = true;if(i % primes[j] == 0) break;}}
}

当某一个数 i i i 是质数时，根据欧拉函数的公式，我们容易得出，这个质数的欧拉函数 ϕ ( i ) = i − 1 \phi(i)=i-1 ϕ(i)=i−1；而当我们在通过 p j × i p_j \times i pj×i 来筛数时。分两种情况

当 i m o d p j = 0 i \mod p_j =0 imodpj=0 时

容易推出 ϕ ( p j × i ) = p j × ϕ ( i ) \phi(p_j \times i) = p_j \times \phi(i) ϕ(pj×i)=pj×ϕ(i)

因为一个数 N N N 的欧拉函数 ϕ ( N ) = N × ( 1 − 1 p 1 ) × . . . × ( 1 − 1 p n ) \phi(N) = N \times (1-\frac{1}{p_1}) \times ... \times (1-\frac{1}{p_n}) ϕ(N)=N×(1−p11)×...×(1−pn1)

假设 i i i 一共有 k k k 个质因子 p 1 p_1 p1 ， p 2 p_2 p2，…， p k p_k pk

则 ϕ ( i ) = i × ( 1 − 1 p 1 ) × . . . × ( 1 − 1 p k ) \phi(i) = i \times (1-\frac{1}{p_1}) \times ... \times (1-\frac{1}{p_k}) ϕ(i)=i×(1−p11)×...×(1−pk1)

而 i m o d p j = 0 i \mod p_j =0 imodpj=0 ，说明 p j p_j pj 是 i i i 的一个质因子。则 p j × i p_j \times i pj×i 这个数的质因子也是 p 1 p_1 p1， p 2 p_2 p2，…， p k p_k pk

则 ϕ ( p j × i ) = p j × i × ( 1 − 1 p 1 ) × . . . × ( 1 − 1 p k ) = p j × ϕ ( i ) \phi(p_j \times i) = p_j \times i \times (1-\frac{1}{p_1}) \times ... \times (1-\frac{1}{p_k}) = p_j \times \phi(i) ϕ(pj×i)=pj×i×(1−p11)×...×(1−pk1)=pj×ϕ(i)
当 i m o d p j ≠ 0 i \mod p_j \ne 0 imodpj=0 时

根据类似上面的过程，容易推出

ϕ ( p j × i ) = p j × i × ( 1 − 1 p 1 ) × . . . × ( 1 − 1 p j ) × . . . × ( 1 − 1 p k ) = ( p j − 1 ) × ϕ ( i ) \phi(p_j \times i) = p_j \times i \times (1-\frac{1}{p_1}) \times ... \times (1-\frac{1}{p_j}) \times ... \times (1-\frac{1}{p_k}) = (p_j-1) \times \phi(i) ϕ(pj×i)=pj×i×(1−p11)×...×(1−pj1)×...×(1−pk1)=(pj−1)×ϕ(i)

p j × i p_j \times i pj×i 这个数的质因子，比 i i i 这个数的质因子，多一个 p j p_j pj

由于在线性筛法的执行过程中，对于质数会保留，对于合数会用其最小质因子筛掉。所以线性筛法是会访问到所有数的。而根据上面的推导，在遇到每种情况时，我们都能求出欧拉函数

当这个数是质数： ϕ ( i ) = i − 1 \phi(i)=i-1 ϕ(i)=i−1
当这个数是合数

其一定会被某个 p j × i p_j \times i pj×i 筛掉
- 当 i m o d p j = 0 i \mod p_j = 0 imodpj=0，则 ϕ ( p j × i ) = p j × ϕ ( i ) \phi(p_j \times i) = p_j \times \phi(i) ϕ(pj×i)=pj×ϕ(i)
- 当 i m o d p j ≠ 0 i \mod p_j \ne 0 imodpj=0，则 ϕ ( p j × i ) = ( p j − 1 ) × ϕ ( i ) \phi(p_j \times i) = (p_j-1) \times \phi(i) ϕ(pj×i)=(pj−1)×ϕ(i)
特殊的， ϕ ( 1 ) = 1 \phi(1)=1 ϕ(1)=1（ 1 1 1 既不是质数也不是合数）

于是，我们便可以通过线性筛法来求解出 1 1 1 到 N N N 的所有数的欧拉函数，时间复杂度 O ( N ) O(N) O(N)

练习题2：acwing - 874: 筛法求欧拉函数

#include<iostream>
using namespace std;typedef long long LL;const int N = 1e6 + 10;int primes[N], ctn;bool st[N];LL phi[N];void get_eulers(int n) {phi[1] = 1;for(int i = 2; i <= n; i++) {if(!st[i]) {primes[ctn++] = i;phi[i] = i - 1;}for(int j = 0; primes[j] <= n / i; j++) {st[primes[j] * i] = true;if(i % primes[j] == 0) {phi[primes[j] * i] = primes[j] * phi[i];break;}phi[primes[j] * i] = (primes[j] - 1) * phi[i];}}
}int main() {int n;scanf("%d", &n);get_eulers(n);LL sum = 0;for(int i = 1; i <= n; i++) sum += phi[i];printf("%lld", sum);return 0;
}

欧拉定理

说了这么多，那么欧拉函数有什么实际的作用或者应用场景呢？

有一个定理叫做欧拉定理，它描述这样一个规律：

若 a a a 与 n n n 互质，则 a ϕ ( n ) m o d n = 1 a^{\phi(n)} \mod n = 1 aϕ(n)modn=1

比如 a = 5 a = 5 a=5 ， n = 6 n = 6 n=6 ，由于 ϕ ( 6 ) = 6 × ( 1 − 1 2 ) × ( 1 − 1 3 ) = 2 \phi(6)=6 \times (1-\frac{1}{2}) \times (1-\frac{1}{3})=2 ϕ(6)=6×(1−21)×(1−31)=2

则有 a ϕ ( n ) m o d n = 5 2 m o d 6 = 1 a^{\phi(n)} \mod n = 5^2 \mod 6 = 1 aϕ(n)modn=52mod6=1

那么欧拉定理如何证明呢？

yxc是这样讲的：

首先， 1 1 1 到 n n n 以内与 n n n 互质的数，有 ϕ ( n ) \phi(n) ϕ(n) 个，将其编号如下

a 1 a_1 a1， a 2 a_2 a2，…， a ϕ ( n ) a_{\phi(n)} aϕ(n)

而 a a a 与 n n n 互质，则可以得知，在 m o d n \mod n modn 的语义下。

a × a 1 a \times a_1 a×a1， a × a 2 a \times a_2 a×a2，…， a × a ϕ ( n ) a \times a_{\phi(n)} a×aϕ(n)

这个集合，与上面的与 n n n 互质的数是同一个集合

进而有，两个集合中全部数的乘积，在 m o d n \mod n modn 的语义下是相等的

a 1 × a 2 × . . . × a ϕ ( n ) ≡ a ϕ ( n ) × a 1 × a 2 × . . . × a ϕ ( n ) a_1 \times a_2 \times ... \times a_{\phi(n)} \equiv a^{\phi(n)} \times a_1 \times a_2 \times ... \times a_{\phi(n)} a1×a2×...×aϕ(n)≡aϕ(n)×a1×a2×...×aϕ(n)

将两边的 a 1 × a 2 × . . . × a ϕ ( n ) a_1 \times a_2 \times ... \times a_{\phi(n)} a1×a2×...×aϕ(n) 消掉，则有

a ϕ ( n ) ≡ 1 a^{\phi(n)} \equiv 1 aϕ(n)≡1，也即 a ϕ ( n ) m o d n = 1 a^{\phi(n)} \mod n = 1 aϕ(n)modn=1

（其实我没有完全把推导的原理搞懂，这个有待后续补充，TODO）

特殊的，当 n n n 是一个质数时（此时我们将 n n n 替换成 p p p 来表示），则有

a ϕ ( p ) m o d p = 1 a^{\phi(p)} \mod p = 1 aϕ(p)modp=1

即 a p − 1 m o d p = 1 a^{p-1} \mod p = 1 ap−1modp=1

这个欧拉定理的特例，被称为费马小定理。

快速幂

快速幂，是用来快速求解出 a k m o d p a^k \mod p akmodp 的结果，时间复杂度为 O ( log ⁡ k ) O(\log k) O(logk)，其中 a a a， k k k， p p p 都可以在 1 0 9 10^9 109 内，如果是按照朴素做法的话，则需要 O ( k ) O(k) O(k) 的时间复杂度，如果 k = 1 0 9 k = 10^9 k=109 ，则这个复杂度就比较高了，而如果是 O ( log ⁡ 2 k ) O(\log_2k) O(log2k)，即便 k = 1 0 9 k=10^9 k=109，也大概只需要 30 30 30 次计算即可，非常的快。

那么，快速幂是怎么做到的呢？

快速幂的核心思路是：反复平方法（思想上有点类似逆向二分。二分是每次在当前基础上减一半，快速幂是每次在当前基础上扩大一倍）。

原理说明如下

我们先预处理出来如下的值： a 2 0 m o d p a^{2^0} \mod p a20modp， a 2 1 m o d p a^{2^1} \mod p a21modp， a 2 2 m o d p a^{2^2} \mod p a22modp ，…， a 2 log ⁡ 2 k m o d p a^{2^{\log_2k}} \mod p a2log2kmodp

一共 log ⁡ 2 k \log_2k log2k 个数。随后，我们通过这些数，组合出 a k a^k ak

即，将 a k a^k ak 拆成 a k = a 2 i × a 2 j × . . . = a 2 i + 2 j + . . . a^k=a^{2^i} \times a^{2^j} \times ...=a^{2^i+2^j+...} ak=a2i×a2j×...=a2i+2j+... ，即 k = 2 i + 2 j + . . . k = 2^i+2^j+... k=2i+2j+...

即，将 k k k 拆成 2 0 2^0 20， 2 1 2^1 21， 2 2 2^2 22，…， 2 log ⁡ 2 k 2^{\log_2k} 2log2k 的某些数的和

其实，就只需要把 k k k 转化为二进制即可。

比如 k = ( 54 ) 10 = ( 110110 ) 2 = 2 1 + 2 2 + 2 4 + 2 5 k=(54)_{10}=(110110)_2=2^1+2^2+2^4+2^5 k=(54)10=(110110)2=21+22+24+25

由于一个十进制的数 k k k ，一定可以转化为二进制，即， k k k 一定能拆成若干个 2 i 2^i 2i 的加和，所以 a k a^k ak 一定能按照上面的式子拆开，所以快速幂算法能够生效。

预处理一共计算出 log ⁡ 2 k \log_2k log2k 个数，需要计算 l o g 2 k log_2k log2k 次，将 k k k 拆成二进制表示，并计算结果，需要 l o g 2 k log_2k log2k 次，所以总共的时间复杂度就是 O ( log ⁡ 2 k ) O(\log_2k) O(log2k)。其实编写代码时，可以将上面两步合在一起，实际只需要 log ⁡ 2 k \log_2k log2k 次运算

举个比较简单例子，比如要计算 9 7 9^7 97，对于指数 7 7 7 ，我们可以将其写成 2 2 + 2 1 + 2 0 2^2+2^1+2^0 22+21+20，即二进制的 ( 111 ) 2 (111)_2 (111)2。则计算 9 7 9^7 97 ，相当于计算 9 2 2 + 2 1 + 2 0 9^{2^2+2^1+2^0} 922+21+20，拆一下，得 9 2 2 × 9 2 1 × 9 2 0 9^{2^2} \times 9^{2^1} \times 9^{2^0} 922×921×920，每相邻两项之间都是平方的关系，则我们可以通过反复计算平方的方法，先计算 9 9 9，再计算 9 2 9^2 92，再计算 ( 9 2 ) 2 (9^2)^2 (92)2，再把每一项乘起来即可。

再比如，计算 9 11 9^{11} 911，将 11 11 11写成二进制，为 ( 1011 ) 2 (1011)_2 (1011)2，从最低位开始，第一位是1，需要乘个 9 9 9，第二位也是 1 1 1，则再乘个 9 2 9^2 92，第三位是 0 0 0 ，则不乘，但是我们仍然要计算一下 9 4 9^4 94 （对上一轮的 9 2 9^2 92 做平方），然后第四位，是1，再乘个 9 8 9^8 98 （对上一轮 9 4 9^4 94 做平方）。

是不是有点类似二分的思路呢？将原先的 n n n 次乘法运算，通过反复平方，压缩成了 l o g 2 n log_2n log2n 次

练习题：acwing - 875: 快速幂

#include<iostream>
using namespace std;typedef long long LL;// 快速幂求解 a^k mod p
int qmi(int a, int k, int p) {int res = 1;// 求 k 的二进制表示while(k > 0) {if(k & 1 == 1) res = (LL) res * a % p;k = k >> 1;a = (LL)a * a % p;}return res;
}int main() {int n;scanf("%d", &n);while(n--) {int a, k, p;scanf("%d%d%d", &a, &k, &p);printf("%d\n", qmi(a, k, p));}return 0;
}

练习题：acwing - 876: 快速幂求逆元

乘法逆元的定义：若两个数 b b b 和 m m m 互质，则对任意整数 a a a，如果 b b b 能整除 a a a，则存在一个整数 x x x ，使得 a b ≡ a × x \frac{a}{b} \equiv a \times x ba≡a×x， m o d m \mod m modm

则称 x x x 是 b b b 在模 m m m 下的乘法逆元，记作 b − 1 m o d m b^{-1} \mod m b−1modm

其实，通俗的说，就是，对于两个互质的数 b b b 和 m m m，一定存在一个数 x x x，使得所有能够被 b b b 整除的数 a a a，都有 a b ≡ a × x \frac{a}{b} \equiv a \times x ba≡a×x ，则称 x x x 是 b b b 在模 m m m 下的乘法逆元，记作 b − 1 m o d m b^{-1} \mod m b−1modm

容易得到： b × b − 1 ≡ 1 b \times b^{-1} \equiv 1 b×b−1≡1，即 b × b − 1 m o d m = 1 b \times b^{-1} \mod m = 1 b×b−1modm=1

求一个数 b b b 在模 m m m 下的逆元，即求一个数 x x x，满足 b × x m o d m = 1 b \times x \mod m = 1 b×xmodm=1 即可

若 m m m 是一个质数，将其记为 p p p，则根据上面的费马小定理，有 b p − 1 m o d p = 1 b^{p-1} \mod p= 1 bp−1modp=1，则 b b b 的逆元 b − 1 m o d p b^{-1} \mod p b−1modp就等于 b p − 2 m o d p b^{p-2} \mod p bp−2modp

若 m m m 不是一个质数（但注意 b b b 和 m m m 是互质的），则根据欧拉定理，有 b ϕ ( m ) m o d m = 1 b^{\phi(m)} \mod m = 1 bϕ(m)modm=1，则 b b b 的逆元 b − 1 m o d p b^{-1} \mod p b−1modp 就等于 b ϕ ( m ) − 1 m o d p b^{\phi(m)-1} \mod p bϕ(m)−1modp

这道题考察的就是快速幂+费马小定理

#include<iostream>
using namespace std;typedef long long LL;int qmi(int a, int k, int p) {int res = 1;while(k > 0) {if(k & 1 == 1) res = (LL) res * a % p;k = k >> 1;a = (LL) a * a % p;}return res;
}int main() {int n;scanf("%d", &n);while(n--) {int a, p;scanf("%d%d", &a, &p);if(a % p == 0) printf("impossible\n"); // a 和 p 不互质else printf("%d\n", qmi(a, p - 2, p));}return 0;
}

扩展欧几里得算法

裴蜀定理：对任意的一对正整数 a a a， b b b，一定存在一对非零整数 x x x， y y y，使得 a x + b y = g c d ( a , b ) ax + by = gcd(a,b) ax+by=gcd(a,b)

令 a x + b y = d ax + by = d ax+by=d，则 d d d 一定是 g c d ( a , b ) gcd(a,b) gcd(a,b) 的倍数。这个是显而易见的，令 a a a 和 b b b 的最大公约数是 c c c，即令 g c d ( a , b ) = c gcd(a,b)=c gcd(a,b)=c，则 a a a 一定是 c c c 的倍数， b b b 也一定是 c c c 的倍数，故 a x + b y ax+by ax+by 也一定是 c c c 的倍数。则最小可以凑出的倍数就是 1 1 1。所以裴蜀定理是成立的。

那么给定一对正整数 a a a， b b b，如何求解出一对 x x x， y y y，使得 a x + b y = g c d ( a , b ) ax+by=gcd(a,b) ax+by=gcd(a,b) 成立呢？

求解 x x x， y y y 的过程，就可以采用扩展欧几里得算法。

扩展欧几里得算法，是在欧几里得算法上的扩展。而欧几里得算法，就是前面求解最大公约数时，用到的辗转相除法。

练习题：acwing - 877: 扩展欧几里得算法

#include<iostream>
using namespace std;int gcd(int a, int b, int &x, int &y) {if(b == 0) {x = 1;y = 0;return a;} else {int d = gcd(b, a % b, y, x); // 注意这里要交换 x 和 y 的位置y -= a / b * x;return d;}
}int main() {int n;scanf("%d", &n);while(n--) {int a, b, x, y;scanf("%d%d", &a, &b);gcd(a, b, x ,y);printf("%d %d\n", x, y);}return 0;
}

对代码的解释：

递归到当 b = 0 b=0 b=0 时，找到 g c d ( a , b ) = a gcd(a,b)=a gcd(a,b)=a，所以 a x + b y = a ax+by=a ax+by=a，则此时更新 x = 1 ， y = 0 x=1，y=0 x=1，y=0。

当 b ≠ 0 b \ne 0 b=0 时，先计算 d = g c d ( b , a m o d b ) d = gcd(b, a \mod b) d=gcd(b,amodb)，并且传入 x x x 和 y y y ，注意 x x x 和 y y y 的位置需要交换一下，然后此时有 b y + ( a m o d b ) x = d by + (a \mod b)x = d by+(amodb)x=d，展开得 b y + ( a − ⌊ a b ⌋ b ) x = d by + (a-\lfloor\frac{a}{b}\rfloor b)x = d by+(a−⌊ba⌋b)x=d，变换一下得 a x + b ( y − ⌊ a b ⌋ x ) ax+b(y-\lfloor\frac{a}{b}\rfloor x) ax+b(y−⌊ba⌋x)

则 a a a 的系数 x x x 不更新， b b b 的系数 y y y 需要更新为 y − ⌊ a b ⌋ x y-\lfloor\frac{a}{b}\rfloor x y−⌊ba⌋x

由于我们在递归时会交换 x x x 和 y y y 的位置，则实际两个数都会被交替更新。即我们在执行 gcd(b, a % b, y, x)这一句之后， x x x 和 y y y 都已经满足 b y + ( a m o d b ) x = d by + (a \mod b)x = d by+(amodb)x=d，我们再通过 b b b 与 a m o d b a \mod b amodb 的系数，来更新 a a a 和 b b b 的系数。

注意， x x x 和 y y y 不是唯一的。由于我们的等式是 a x + b y = g c d ( a , b ) ax+by=gcd(a,b) ax+by=gcd(a,b)，其中 a a a， b b b， g c d ( a , b ) gcd(a,b) gcd(a,b) 都是固定的常数，所以这个式子其实可以转化为我们中学时学过的一次线性函数 y = a x + b y=ax+b y=ax+b 的形式，这个函数的图像是一根直线，线上有多个点 ( x , y ) (x,y) (x,y) 都满足这个等式，所以 ( x , y ) (x,y) (x,y) 是有多组的（准确的说，有无穷多组）

我们只要求解出一组 ( x 0 , y 0 ) (x_0,y_0) (x0,y0) 满足上面的等式。就可以求出所有的 ( x ， y ) (x，y) (x，y)

a x + b y = d ax+by=d ax+by=d 可以写成 a ( x − b d ) + b ( y + a d ) = d a(x-\frac{b}{d}) + b(y+\frac{a}{d})=d a(x−db)+b(y+da)=d

通解如下： x = x 0 − b d × k x = x_0-\frac{b}{d} \times k x=x0−db×k， y = y 0 − a d × k y=y_0-\frac{a}{d} \times k y=y0−da×k，其中 k k k 是任意整数

扩展欧几里得算法的应用 —— 求解线性同余方程

练习题：acwing - 878: 线性同余方程

即，给定 a a a， b b b， m m m，求解一个 x x x，使得满足 a x ≡ b ax \equiv b ax≡b ( m o d m ) (\mod m) (modm)

因为上面等式的含义是： a x ax ax 除以 m m m，余数是 b b b ，所以 a x ax ax 一定是 m m m 的某个倍数，加上 b b b，即 a x = m y + b ax = my+b ax=my+b

再变形一下，得 a x − m y = b ax-my=b ax−my=b，令 y ′ = − y y^{'} = -y y′=−y，则有 a x + m y ′ = b ax+my^{'}=b ax+my′=b

这样就转变成了上面使用扩展欧几里得算法能够求解问题，只需要保证 b b b 是 a a a 和 m m m 的最大公约数的倍数即可，否则无解。

为什么要保证 b b b 是 a a a 和 m m m 的最大公约数的倍数呢？解释如下

假设 a a a 和 m m m 的最大公约数是 k k k，那么根据上面的裴属定理，一定有一组 ( x , y ) (x,y) (x,y) ，使得 a x + m y = k ax+my=k ax+my=k 成立，那么对于 k k k 的倍数，比如 t t t 倍，就是 t k tk tk，则一定有 a t x + m t y = t k atx+mty = tk atx+mty=tk。所以，只要保证 b b b 是 a a a 和 m m m 的最大公约数的某个倍数，这个线性同余方程就有解，否则无解。

#include<iostream>
using namespace std;typedef long long LL;int gcd(int a, int b, int &x, int &y) {if(b == 0) {x = 1;y = 0;return a;} else {int d = gcd(b, a % b, y, x); // 注意这里要交换 x 和 y 的位置y -= a / b * x;return d;}
}int main() {int n;scanf("%d", &n);while(n--) {int a, b, m, x, y;scanf("%d%d%d", &a, &b, &m);int d = gcd(a, m, x ,y);if(b % d != 0) printf("impossible\n");else printf("%d\n", (LL)x * b / d % m);}return 0;
}

中国剩余定理

例子：有一个数，它除以3余2，除以5余3，除以7余2，问这个数是多少？

即 x ≡ 2 ( m o d 3 ) x \equiv 2(\mod 3) x≡2(mod3)， x ≡ 3 ( m o d 5 ) x \equiv 3(\mod 5) x≡3(mod5)， x ≡ 2 ( m o d 7 ) x \equiv 2(\mod 7) x≡2(mod7)

该问题的求解步骤如下：

找出3和5的公倍数中，除以7余1的最小数，即15；找出3和7的公倍数中，除以5余1的最小数，即21；找出5和7的公倍数中，除以3余1的最小数，即70。
用15乘以2（除7余2的2），用21乘以3（除5余3的3），用70乘以2（除3余2的2），把3个乘积加起来；得 15 × 2 + 21 × 3 + 70 × 2 = 233 15 \times 2 + 21 \times 3 + 70 \times 2 = 233 15×2+21×3+70×2=233
用 233 除以 3，5，7的最小公倍数 105 ，得到余数 23，答案即为23

这个求解过程，被称为中国剩余定理。更一般的描述如下：

给定一堆两两互质的数 m 1 m_1 m1， m 2 m_2 m2，…， m k m_k mk

存在一个数 x x x，满足

x ≡ a 1 ( m o d m 1 ) x \equiv a_1 (\mod m_1) x≡a1(modm1)

x ≡ a 2 ( m o d m 2 ) x \equiv a_2 (\mod m_2) x≡a2(modm2)

…

x ≡ a k ( m o d m k ) x \equiv a_k (\mod m_k) x≡ak(modmk)

求解这个 x x x，如下：

令 M = m 1 × m 2 × . . . × m k M = m_1 \times m_2 \times ... \times m_k M=m1×m2×...×mk ，令 M i = M m i M_i=\frac{M}{m_i} Mi=miM，即 M i M_i Mi 是除了 m i m_i mi 之外其他所有数的乘积。

而由于对于 j ∈ [ 1 , k ] ， j \in [1,k]， j∈[1,k]， m j m_j mj 之间两两互质，所以 M i M_i Mi 和 m i m_i mi 互质。

而根据上面介绍的欧拉定理和乘法逆元，对于两个互质的数 M i M_i Mi 和 m i m_i mi，一定存在一个乘法逆元 M i − 1 M_i^{-1} Mi−1 ，使得 M i × M i − 1 ≡ 1 ( m o d m i ) M_i \times M_i^{-1} \equiv 1 (\mod m_i) Mi×Mi−1≡1(modmi)

则 x = a 1 × M 1 × M 1 − 1 + a 2 × M 2 × M 2 − 1 + . . . + a k × M k × M k − 1 x = a_1 \times M_1 \times M_1^{-1} + a_2 \times M_2 \times M_2^{-1} + ... + a_k \times M_k \times M_k^{-1} x=a1×M1×M1−1+a2×M2×M2−1+...+ak×Mk×Mk−1 ，这个求解公式，就是中国剩余定理

特殊的，满足中国剩余定理的最小解是：
x = ( a 1 × M 1 × M 1 − 1 + a 2 × M 2 × M 2 − 1 + . . . + a k × M k × M k − 1 ) m o d M x = (a_1 \times M_1 \times M_1^{-1} + a_2 \times M_2 \times M_2^{-1} + ... + a_k \times M_k \times M_k^{-1}) \mod M x=(a1×M1×M1−1+a2×M2×M2−1+...+ak×Mk×Mk−1)modM

上面式子中的乘法逆元，可以用扩展欧几里得算法求解

参考链接：中国剩余定理学习笔记
（强推这篇博文，讲解的十分清晰）

练习题：acwing - 2024: 表达整数的奇怪方式

题解：

这道题，由于对a和m没有任何限制，所以不能完全套用中国剩余定理来求解。

先考虑2个线性同余等式

x m o d a 1 = m 1 x \mod a_1 = m_1 xmoda1=m1， x m o d a 2 = m 2 x \mod a_2 =m_2 xmoda2=m2

可以写成如下的式子

x = k 1 a 1 + m 1 x = k_1a_1 + m_1 x=k1a1+m1， x = k 2 a 2 + m 2 x = k_2a_2 + m_2 x=k2a2+m2

则 k 1 a 1 + m 1 = k 2 a 2 + m 2 k_1a_1 + m_1 = k_2a_2 + m_2 k1a1+m1=k2a2+m2

移项并整理一下，得 k 1 a 1 − k 2 a 2 = m 2 − m 1 k_1a_1-k_2a_2=m_2-m_1 k1a1−k2a2=m2−m1

这就可以使用扩展欧几里得算法来求解，当然，上式有解的条件是 g c d ( a 1 , a 2 ) ∣ m 2 − m 1 gcd(a_1,a_2) | m_2-m_1 gcd(a1,a2)∣m2−m1，即 m 2 − m 1 m_2-m_1 m2−m1 能够被 a 1 a_1 a1 和 a 2 a_2 a2 的最大公约数整除。

用扩展欧几里得算法求出上式的一组解 k 10 k_{10} k10， k 20 k_{20} k20，随后可求出一组通解为

k 1 = k 10 + k a 2 d k_1 = k_{10}+k\frac{a_2}{d} k1=k10+kda2， k 2 = k 20 + k a 1 d k_2=k_{20}+k\frac{a_1}{d} k2=k20+kda1

其中 k k k 为任意整数， d d d 为 g c d ( a 1 , a 2 ) gcd(a_1,a_2) gcd(a1,a2)

则 x = ( k 10 + k a 2 d ) a 1 + m 1 = a 1 k 10 + m 1 + k a 1 a 2 d x = (k_{10}+k\frac{a_2}{d})a_1+m_1=a_1k_{10}+m_1+k\frac{a_1a_2}{d} x=(k10+kda2)a1+m1=a1k10+m1+kda1a2

令 x 0 = a 1 k 10 + m 1 x_0=a_1k_{10}+m_1 x0=a1k10+m1， a = a 1 a 2 d a=\frac{a_1a_2}{d} a=da1a2，则 x = x 0 + k a x=x_0+ka x=x0+ka

这样就将2个不定方程，变成了1个不定方程。

只需要执行n-1次合并，就能将n个不定方程，变成1个不定方程。

最后就是求解一个形如 x = x 0 + k a x=x_0 + ka x=x0+ka 的不定方程，其中 x 0 x_0 x0 和 a a a 是固定值

#include <iostream>
#include <algorithm>
using namespace std;typedef long long LL;// 扩展欧几里得算法
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, y, x);y -= a / b * x;return d;
}int main() {int n;scanf("%d", &n);bool no_ans = false;LL a1, m1;cin >> a1 >> m1;for (int i = 0; i < n - 1; i++) {LL a2, m2;cin >> a2 >> m2;LL k1, k2;LL d = exgcd(a1, a2, k1, k2);if ((m2 - m1) % d != 0) {no_ans = true;break;}k1 *= (m2 - m1) / d;LL t = a2 / d;k1 = (k1 % t + t) % t; // 将k1置为最小的正整数解// 更新a1和m1, 准备下一轮的合并m1 = a1 * k1 + m1;a1 = abs(a1 / d * a2);}if (no_ans) printf("-1");else printf("%lld", (m1 % a1 + a1) % a1);
}

（2021/07/30 更新：补充了中国剩余定理的部分内容，还剩一道练习题的题解）
（2021/08/03 更新：补充了中国剩余定理的练习题题解）