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×(1P11)×(1P21)×...×(1Pn1),这就是欧拉函数的求解公式

比如 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×(121)×(131)=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 1 N N N 中,去掉其全部质因子 P 1 , P 2 , . . . , P n P_1,P_2,...,P_n P1P2,...,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} NP1NP2N...PnN

  2. 加上所有 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+...+Pn1×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 减一次,一共减了两次 ,多减了一次,那么我们需要加一次加回来。

  3. 加上所有 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} S2P1×P2×P3NP1×P2×P4N...Pn2×Pn1×PnN

    即,对 i , j , k ∈ [ 1 , n ] i,j,k \in [1,n] i,j,k[1,n] 的全部组合,进行减

  4. 后续的步骤同理

    若若某个数,同时是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+....+Pn3×Pn2×Pn1×PnN

  5. 对于那些是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} NP1NP2N...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+...+Pn1×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×P3NP1×P2×P4N...Pn2×Pn1×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+....+Pn3×Pn2×Pn1×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×(1P11)×(1P21)×...×(1Pn1)

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

练习题: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×(1P11) 这样的式子,计算时我们采用 N P 1 × ( P 1 − 1 ) \frac{N}{P_1} \times (P_1-1) P1N×(P11),这样能保证得到的结果是整数,而不出现小数。

由于 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×(P11) 就一定是整数,令 S = N P 1 × ( P 1 − 1 ) S = \frac{N}{P_1} \times (P_1-1) S=P1N×(P11),那么在第二轮迭代中,我们需要计算 S × ( 1 − 1 P 2 ) S \times (1-\frac{1}{P_2}) S×(1P21),同样的,我们计算 S P 2 × ( P 2 − 1 ) \frac{S}{P_2} \times (P_2-1) P2S×(P21),那么 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=P1k11×P2k2×...×Pnkn×(P11),所以 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×(Pi1),这样能保证整个计算的过程中,不出现小数,即可保证最终结果的正确性。

筛法

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

然而,有的时候,我们需要求解某一个范围内全部数的欧拉函数(比如求解 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)=i1;而当我们在通过 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×(1p11)×...×(1pn1)

    假设 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×(1p11)×...×(1pk1)

    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×(1p11)×...×(1pk1)=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×(1p11)×...×(1pj1)×...×(1pk1)=(pj1)×ϕ(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)=i1

  • 当这个数是合数

    其一定会被某个 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)=(pj1)×ϕ(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×(121)×(131)=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 ap1modp=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 baa×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 b1modm

其实,通俗的说,就是,对于两个互质的数 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 baa×x ,则称 x x x b b b 在模 m m m 下的乘法逆元,记作 b − 1 m o d    m b^{-1} \mod m b1modm

容易得到: b × b − 1 ≡ 1 b \times b^{-1} \equiv 1 b×b11,即 b × b − 1 m o d    m = 1 b \times b^{-1} \mod m = 1 b×b1modm=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 bp1modp=1,则 b b b 的逆元 b − 1 m o d    p b^{-1} \mod p b1modp就等于 b p − 2 m o d    p b^{p-2} \mod p bp2modp

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 b1modp 就等于 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=1y=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+(abab)x=d,变换一下得 a x + b ( y − ⌊ a b ⌋ x ) ax+b(y-\lfloor\frac{a}{b}\rfloor x) ax+b(ybax)

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

由于我们在递归时会交换 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) (xy)

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(xdb)+b(y+da)=d

通解如下: x = x 0 − b d × k x = x_0-\frac{b}{d} \times k x=x0db×k y = y 0 − a d × k y=y_0-\frac{a}{d} \times k y=y0da×k,其中 k k k 是任意整数

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

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

即,给定 a a a b b b m m m,求解一个 x x x,使得满足 a x ≡ b ax \equiv b axb ( 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 axmy=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) x2(mod3) x ≡ 3 ( m o d    5 ) x \equiv 3(\mod 5) x3(mod5) x ≡ 2 ( m o d    7 ) x \equiv 2(\mod 7) x2(mod7)

该问题的求解步骤如下:

  1. 找出3和5的公倍数中,除以7余1的最小数,即15;找出3和7的公倍数中,除以5余1的最小数,即21;找出5和7的公倍数中,除以3余1的最小数,即70。
  2. 用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
  3. 用 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) xa1(modm1)

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

x ≡ a k ( m o d    m k ) x \equiv a_k (\mod m_k) xak(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} Mi1 ,使得 M i × M i − 1 ≡ 1 ( m o d    m i ) M_i \times M_i^{-1} \equiv 1 (\mod m_i) Mi×Mi11(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×M11+a2×M2×M21+...+ak×Mk×Mk1 ,这个求解公式,就是中国剩余定理

特殊的,满足中国剩余定理的最小解是:
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×M11+a2×M2×M21+...+ak×Mk×Mk1)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 k1a1k2a2=m2m1

这就可以使用扩展欧几里得算法来求解,当然,上式有解的条件是 g c d ( a 1 , a 2 ) ∣ m 2 − m 1 gcd(a_1,a_2) | m_2-m_1 gcd(a1,a2)m2m1,即 m 2 − m 1 m_2-m_1 m2m1 能够被 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 更新:补充了中国剩余定理的练习题题解)

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值