基础数论笔记

学完基础数论的大部分内容，过来写一个总结，巩固记忆，另外再重新推倒一些公式。

本文章主要涉及的模块：

• 素数
• 最大公因数及最小公倍数
• 快速幂及扩展欧拉定理
• 扩展欧几里得与同余
• 逆元
• 中国剩余定理

素数

定义

因数只有 \(1\) 和这个数本身的数被称作素数。

注意：\(1\) 既不是素数也不是合数，\(2\) 是最小的素数。

两个关于素数的定理

唯一分解定理

对于任意大于 \(1\) 的整数 \(x\)，都可以分解成若干个素数的乘积：

\(x=p_1^{a_1}\times p_2^{a_2}\times p_3^{a_3}\times \cdots \times p_n^{a_n}(a_i \in \Z^+)\)

比如 \(57=3 \times19\)，\(1704=2^3 \times 3\times 71\cdots\)。

费马小定理

若存在正整数 \(a\) 与素数 \(p\)，且 \(\gcd(a,p)=1\)。则有 \(a^{p-1}\equiv 1\pmod p\)。

证明就不给了，有兴趣可以上网找一下。费马小定理在后续会有很多用处。

素数的判别方法

试除法

试判断 \(x\) 是否为素数。

由素数的定义可知，素数的因数只有两个。所以我们可以直接暴力枚举除 \(1\) 以外小于 \(x\) 的数。

• 如果有 \(1<i<x\) 且 \(i \mid x\) 就说明 \(x\) 有超过两个数的因数，所以 \(x\) 不是素数。
• 如果枚举完了，但还是找不到整除 \(x\) 的数，就说明 \(x\) 是素数。

时间复杂度为 \(O(n)\)，能不能优化呢？答案是可以的。

由因数的性质可知，若正整数 \(n\) 有一因数 \(i\) 且 \(i\le\sqrt{n}\)，则 \(n\) 必有另外一因数 \(j\) 且 \(j\ge\sqrt{n}\)。在 \(i=\sqrt{n}\) 时，有 \(i=j\)。

根据这一个性质，我们可以只枚举 \(2\) 到 \(\sqrt{x}\) 的数进行试除，时间复杂度优化到 \(O(\sqrt{n})\)。

代码

inline bool is_prime(int x) {
	for(register int i = 2;i <= sqrt(x);i ++)
		if(x % i == 0)
			return false;
	return true;
}

但当题目要求求出一串数中的素数时，试除法的时间复杂度就变成了 \(O(n\sqrt{n})\)，这在 \(n\) 过大时是不能接受的。所以，我们便有了素数的筛法。

埃氏筛

素数的第一种筛法。那何为筛法呢？顾名思义，就是把不是素数的数筛掉，这样剩下的就只有素数了。那怎么筛呢？

由唯一分解定理可知，任意一个大于 \(1\) 的合数都可以表示成多个素数相乘。此时不难想到，对于每一个素数，我们可以把它的所有的倍数筛掉，这样就可以把题目给定范围内的合数都筛掉，只留下素数。记得要特判 \(1\)，因为它既不是素数也不是合数。这样我们便有了素数的第一种筛法。鉴于这个筛法是古希腊数学家埃拉托斯特尼发明的，所以我们便称这个筛法为埃氏筛。

给大家模拟一下，加深理解：

此时，若我们要筛出 \(12\) 以内的素数。（\(false\) 表示是素数，\(true\) 则反之）

1. 首先，我们要特判 \(1\)：

	\(1\)	\(2\)	\(3\)	\(4\)	\(5\)	\(6\)	\(7\)	\(8\)	\(9\)	\(10\)	\(11\)	\(12\)
\(flag\)	\(\color{red}true\)	\(false\)	\(false\)	\(false\)	\(false\)	\(false\)	\(false\)	\(false\)	\(false\)	\(false\)	\(false\)	\(false\)

2. 接着，我们从 \(2\) 开始一个一个的往后筛：

	\(1\)	\(2\)	\(3\)	\(4\)	\(5\)	\(6\)	\(7\)	\(8\)	\(9\)	\(10\)	\(11\)	\(12\)
\(flag\)	\(true\)	\(\color{orange}false\)	\(false\)	\(\color{red}true\)	\(false\)	\(\color{red}true\)	\(false\)	\(\color{red}true\)	\(false\)	\(\color{red}true\)	\(false\)	\(\color{red}true\)

	\(1\)	\(2\)	\(3\)	\(4\)	\(5\)	\(6\)	\(7\)	\(8\)	\(9\)	\(10\)	\(11\)	\(12\)
\(flag\)	\(true\)	\(false\)	\(\color{orange}false\)	\(true\)	\(false\)	\(\color{red}true\)	\(false\)	\(true\)	$$\color{red}true$$	\(true\)	\(false\)	\(\color{red}true\)

3. 此时，我们发现 \(4\) 已经被筛掉了，所以我们直接跳过它去筛 \(5\)。

	\(1\)	\(2\)	\(3\)	\(4\)	\(5\)	\(6\)	\(7\)	\(8\)	\(9\)	\(10\)	\(11\)	\(12\)
\(flag\)	\(true\)	\(false\)	\(false\)	\(true\)	\(\color{orange}false\)	\(true\)	\(false\)	\(true\)	\(true\)	\(\color{red}true\)	\(false\)	\(true\)

4. 我们发现，在筛到 \(7\) 时，它的倍数已经超出了需要求的范围，所以我们便可以停止了。

重新回看表格，所有的 \(false\) 即为素数：

	\(1\)	\(2\)	\(3\)	\(4\)	\(5\)	\(6\)	\(7\)	\(8\)	\(9\)	\(10\)	\(11\)	\(12\)
\(flag\)	\(true\)	\(\color{red}false\)	\(\color{red}false\)	\(true\)	\(\color{red}false\)	\(true\)	\(\color{red}false\)	\(true\)	\(true\)	\(true\)	\(\color{red}false\)	\(true\)

代码

int n;
bool flag[MAXN];

void prime(int n) {
	flag[1] = true;
	for(register int i = 2;i <= n;i ++)
		if(!flag[i])
			for(register int j = 2;j * i <= n;j ++)
				flag[j * i] = true;
	for(register int i = 1;i <= n;i ++)
		if(!flag[i])
			cout << i << " ";
	return;
}

埃氏筛的时间复杂度是：\(O(n \log \log n)\)。由于证明过程过于繁琐，这里就不赘述了，有兴趣可以去搜一下。一般来说，知道时间复杂度就可以了。

线性筛

埃氏筛的时间复杂度已经很优了，但当题目的数据范围达到 \(10^9\) 时，埃氏筛便处理不了了。此时，我们便需要用到线性筛了。

顾名思义，线性筛的时间复杂度在埃氏筛上的基础上优化到了 \(O(n)\)。具体优化在哪里？

不难发现，许多正整数都有着超过 \(1\) 个的素因数。而在埃氏筛中，正整数会被其所有的素因数都筛一遍，这里浪费了很多时间。但其实，我们只需要用每个正整数最小的素因数去筛便可以了，线性筛就在这个点上优化了。

拿刚刚的例子讲解一下：

在小于 \(12\) 的数中有 \(6=2\times 3,10=2\times5,12=2^2\times3\)。这三个数都被筛了多次。

	\(1\)	\(2\)	\(3\)	\(4\)	\(5\)	\(6\)	\(7\)	\(8\)	\(9\)	\(10\)	\(11\)	\(12\)
\(flag\)	\(true\)	\(false\)	\(false\)	\(true\)	\(false\)	\(\color{red}true\)	\(false\)	\(true\)	\(true\)	\(\color{red}true\)	\(false\)	\(\color{red}true\)

而在线性筛中，\(6\)，\(10\)，\(12\)，都只会被 \(2\) 筛掉 \(1\) 次。

具体的实现方法便是只用小于等于自己最小的素因数的素数来筛掉合数。

算法流程：

• 首先与埃氏筛一样，要特判 \(1\) 的情况。

• 然后，把找到的素数都放进一个 prime 数组里面。假设我们现筛到了 \(i\)，则我们便会用 prime 数组中的素数与 \(i\) 相乘，把相乘所得到的结果都标记为合数。当 \(i\) 除以 prime 中的素数余数为 \(0\) 时停止（先筛完再停止）。这样便做到了用最小的素数筛掉合数。

为什么这是正确的呢？我们来证明一下：

设 prime 数组中的素数为 \(k_1\)，\(k_2\)，\(\cdots\)，\(k_n\)，且 \(k_1\le k_2 \le \cdots \le k_n\)。

• 若 \(i\) 为素数。

  当 \(k_i\nmid i\) 时说明 \(k_i\) 比 \(i\) 小，所以 \(k_i\) 是正整数 \(k_i \times i\) 最小的素因数。

  当 \(k_i \mid i\) 时说明 \(k_i=i\)，此时，\(k_i\) 与 \(i\) 同为 \(k_i \times i\) 的最小素因数。

• 若 \(i\) 不为素数，设其分解素因数形式为：\(p_1^{a_1}\times p_2^{a_2} \times \cdots \times p_n^{a_n}\)。

  则当 \(k_i \nmid i\) 时，同样的，\(k_i\) 必然比 \(p_1\) 要小，所以 \(k_i\) 肯定是正整数 \(k_i \times i\) 的最小素因数。

  当 \(k_i \mid i\) 且是第一次整除时说明 \(k_i=p_1\)，此时，\(k_i\) 与 \(p_1\) 同为 \(k_i \times i\) 的最小素因数。

代码

int prime[MAXN] , cnt;
bool is_prime[MAXN];

inline void Init(int x) {
	is_prime[1] = true;
	for(register int i = 2;i <= x;i ++) {
		if(!is_prime[i])
			prime[++ cnt] = i;
		for(register int j = 1;j <= cnt && i * j <= x;j ++) {
			is_prime[i * prime[j]] = true;
			if(i % prime[j])
				break;
		}
	}
	return;
}

欧拉函数

在做关于素数的题目时经常需要用到关于欧拉函数的知识，并且求法相似，所以把它们俩放在一起讲。

求单个数的欧拉函数

定义

对于任意正整数 \(n\)，欧拉函数 \(\varphi(n)\) 表示在不大于 \(n\) 的正整数中与 \(n\) 互素的数的个数。

举几个例子：

• \(\varphi(1)=1\)，即 \(1\)。
• \(\varphi(6)=2\)，即 \(1\)，\(5\)。
• \(\varphi(12)=4\)，即 \(1\)，\(5\)，\(7\)，\(11\)。

怎样求出一个数的欧拉函数呢？公式呈上：

对于任意正整数 \(n\)，若 \(n=P_1^{k_1} \times P_2^{k_2} \times \cdots \times P_m^{k_m}\)。则有 \(\varphi(n)=n\times (1-\frac{1}{P_1})\times (1-\frac{1}{P_2})\times \cdots \times (1-\frac{1}{P_m})\)。

化简得：\(\varphi(n)=\prod_{i=1}^m (1-\frac{1}{p_i})\)。

有这个公式，我们可以求出单个正整数的欧拉函数。即在分解素因数时计算欧拉函数。

时间复杂度：\(O(\sqrt{n})\)。

代码

注意：在素因数分解完后如果 \(n\) 变成了素数，那 \(ans\) 还要再运算一次。

inline int phi(int n) {
	int ans = n;
	for(register int i = 2;i <= sqrt(n);i ++)
		if(n % i == 0) {
			ans = ans / i * (i - 1);
			while (n % i == 0)
				n /= i;
		}
	if(n > 1)
		ans = ans / n * (n - 1);
	return ans;
}

欧拉函数公式的推导过程

接下来，我们一起推导在 \(O(\sqrt{n})\) 时间内求出欧拉函数的公式。

首先要知道几个性质：

1. 对于任意素数 \(x\)，有 \(\varphi(x)=x-1\)。这个很好理解，即 \(x\) 与所有小于它的数互素。

2. 对于任意素数 \(p\)，有 \(\varphi(p^k)= p^k-p^{k-1}\)。为什么呢？

   证明：因为 \(p\) 是素数，所以 \(p^k\) 的素因数只含有 \(p\)。那么在小于等于 \(p^k\) 的数中与其不互素的数便是 \(p\)，\(2\times p\)，\(3\times p\)，\(\cdots\)，\(p^{k-1}\times p\)，共 \(p^{k-1}\) 个。得证。

3. 欧拉函数是积性函数，但不是完全积性函数。这意味着 \(\varphi(m\times n)=\varphi(m)\times \varphi(n)\) 成立的前提是 \(\gcd(m,n)=1\)。

了解完这几个性质后，我们开始推导：

设 \(N=p_1^{k_1}\times p_2^{k_2}\times p_3^{k_3}\times \cdots \times p_n^{k_n}\)。

因为 \(p_1,p_2,p_3,\cdots,p_n\) 都是素数，所以 \(p_1^{k_1}，p_2^{k_2}，\cdots，p_n^{k_n}\) 两两互素。

所以由性质 3 得：

\[\varphi(N)=\varphi(p_1^{k_1})\times \varphi(p_2^{k_2}) \times \cdots \times \varphi(p_n^{k_n}) \]

再由性质 2 得：

\[\varphi(N)=(p_1^{k_1}-p_1^{k_1-1})\times(p_2^{k_2}-p_2^{k_2-1})\times\cdots\times(p_n^{k_n}-p_n^{k_n-1}) \]

提取公因数得：

\[\varphi(N)=p_1^{k_1}\times p_2^{k_2}\times \cdots \times p_n^{k_2} \times (1-\frac{1}{p_1})\times (1-\frac{1}{p_2})\times\cdots\times (1-\frac{1}{p_n}) \]

最后把前 \(n\) 项合并一下便得到：

\[\varphi(N)=N\times(1-\frac{1}{p_1})\times (1-\frac{1}{p_2})\times\cdots\times (1-\frac{1}{p_n}) \]

推导完毕。

欧拉函数的线性筛法

和素数一样，这样的算法求单个数的欧拉函数时间复杂度很优，但当题目需要求一串数的欧拉函数时便处理不了了。此时的时间复杂度便是 \(O(n\sqrt{n})\)。

那求一串数欧拉函数是否和求素数一样，有跟快甚至线性的方法呢？当然有，由刚刚的几个性质，我们可以在线性筛的基础上改进，在 \(O(n)\) 的时间内求得每个数的欧拉函数。

这里我们需要再引入一个性质：

4. 若 \(i\bmod p=0\)，则有：\(\varphi(p\times i)=\varphi(i)\times p\)。

   证明：

   设 \(b=\gcd(n , i)\)；

   因为 \(n\)，\(i\) 不互素，所以 \(n=k_1b\)，\(i=k_2b\)。

   因为 \(n+i=(k_1+k_2)b\)，所以 \(n+i\) 与 \(i\) 不互素。

   \([1,i]\) 中与 \(i\) 不互素的整数 \(n\) 共 \(i-\varphi(i)\) 个。

   因为 \(n+i\) 与 \(i\) 不互素，所以 \([1+i,i+i]\)，即 \((i,2i]\) 中与 \(i\) 不互素的整数也是 \(i-\varphi(i)\) 个，所以 \([1,i+p]\) 中与 \(i\) 不互素的整数 \(p\times i-p\times \varphi(i)\) 个。

   因为 \(i\bmod p = 0\)，\(i\) 与 \(p\) 没有不同的因数，\([1,i\times p]\) 中与 \(i\times p\) 不互素的整数数量 \(p\times i-p\times \varphi(i)=i\times p-\varphi(i\times p)\)。

   所以：\(\varphi(i\times p)=p\times \varphi(i)\)。

我们一起边回顾线性筛代码，边了解求欧拉函数的方法：

定义phi[MAXN]存储欧拉函数。

首先，在开始筛之前我们先需要特判 \(1\) 的情况。同样，在求欧拉函数时，要先把 \(\varphi(1)\) 赋值为 \(1\)。

 is_prime[1] = true;
 phi[1] = 1;

然后我们开始找素数。在找的过程中，如果遇到了素数 \(x\)，那我们就存储下 \(x\)。同时，由之前讲的性质 \(1\) 得 \(\varphi(x)=x-1\)。所以我们便可以对它的欧拉函数进行赋值。

for(register int i = 2;i <= n;i ++) {
	if(!is_prime[i]) {
		prime[++ cnt] = i;
		phi[i] = i - 1;
	}
	...
}

最后，在用找到的数 \(x\) 来筛素数的过程中。

若用来与 \(x\) 相乘的素数 \(y\) 不是 \(x\) 的最小的素因数，即比 \(x\) 的最小的素因数小，则说明 \(x\) 不包含 \(y\) 这个素因数，所以 \(\gcd(x,y)=1\)。此时，由性质 \(3\) 得到：\(\varphi(x\times y)=\varphi(x)\times\varphi(y)\)。

若用来与 \(x\) 相乘的素数 \(y\) 已经与 \(x\) 不互素了，即此时 \(y\) 是 \(x\) 的最小素因数。则根据性质 \(4\) ，有：\(\varphi(x\times y)=\varphi(x)\times y\)。

for(register int j = 1; j <= cnt , i * prime[j] <= n; j ++) {
	is_prime[i * prime[j]] = true;
	if(i % prime[j] == 0) {
		phi[i * prime[j]] = phi[i] * prime[j];
		break;
	}
	phi[i * prime[j]] = phi[i] * phi[prime[j]];
}

完整代码

int phi[MAXN] , prime[MAXN] , cnt;
bool is_prime[MAXN];

inline void Init_phi(int n) {
	is_prime[1] = true;
	phi[1] = 1;
	for(register int i = 2;i <= n;i ++) {
		if(!is_prime[i]) {
			prime[++ cnt] = i;
			phi[i] = i - 1;
		}
		for(register int j = 1;j <= cnt , i * prime[j] <= n;j ++) {
			is_prime[i * prime[j]] = true;
			if(i % prime[j] == 0) {
				phi[i * prime[j]] = phi[i] * prime[j];
				break;
			}
			phi[i * prime[j]] = phi[i] * phi[prime[j]];
		}
	}
	return;
}

最大公因数与最小公倍数

这一章内容比较少，主要讲的是最大公因数的几种求法。

最大公因数

定义

在两个数的公因数中最大的那个叫做这两个数的最大公因数。

我们一般用 \(\gcd(a,b)\) 来表示 \(a\)，\(b\) 两数的最大公因数。

接下来，我给大家介绍 \(4\) 种求最大公因数的方法。

穷举法

顾名思义，即暴力枚举：从两个数中较小的一个开始从大到小枚举，遇到的第一个可以同时整除这两个数的数就是它们的最大公因数。

至于为什么要从较小的一个开始枚举，因为两个数的公因数只可能小于等于较小的那个数。

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

代码

inline int gcd(int a , int b) {
	int x = min(a , b);
	for(register int i = x;i >= 1;i --)
		if(a % i == 0 && b % i == 0)
			return i;
}

辗转相除法

定义

辗转相除法，又称欧几里得算法，是一种用来求最大公因数的算法。其算法流程主要基于这么一个式子：\(\gcd(a,b)=\gcd(b,a\bmod b)\)。当 \(a\bmod b=0\) 时，便可以得到两数的最大公因数，即 \(\gcd(b,a\bmod b)\) 中的 \(b\)。

由这个式子，我们可以想到一种递归求最大公因数的写法：

代码

int gcd(int a , int b) {
	if(!b)
		return a;
	return gcd(b , a % b);
}

当然，也可以把这个过程改写成一个循环：

inline int gcd(int a , int b) {
	int r = a % b;
	while(r) {
		a = b;
		b = r;
		r = a % b;
	}
	return b;
}

时间复杂度均为：\(O(\log \max(a,b))\)。

辗转相除法的证明

下面给出证明：

设 \(d\) 是 \(a\)，\(b\) 的任意一个公因数且 \(a > b\)。

则必有：\(a=kb+r\)；

移项得：\(r=a-kb\)；

两边同时除以 \(d\) 得：\(\frac{r}{d}=\frac{a}{d}-\frac{kb}{d}\)；

因为 \(d\) 是 \(a\) 与 \(b\) 的公因数，所以 \(d\mid a\) 并且 \(d\mid b\)。由此得到：\(d\mid r\)。

而 \(r=a\bmod b\)。

所以，\(d\) 既是 \(a\)，\(b\) 的公因数，也是 \(b\)，\(a\bmod b\) 的公因数。

现在我们证出 \((a,b)\) 与 \((b,a\bmod b)\) 的公因数相等，而由 \(d\) 的任意性得 \(d\) 可为 \(a\) 与 \(b\) 的最大公因数。故得证。

更相减损法

定义

更相减损法，即把两数不断相减，直到两数相等。其实主要思想和辗转相除法类似。

具体算法流程是：

试求 \(a\)，\(b\) 两数的最大公因数。

1. 若 \(2\mid a\) 且 \(2\mid b\)，则先将两数所含的 \(2\) 给除掉。

2. 接着，我们不断用大数减小数，直到两数相等。此时相等的两数乘上约掉的 \(2\) 即为原本两数的最大公因数。

用公式表达就是：

• \(\gcd(2a,2b)=2\gcd(a,b)\)。
• \(a\ge b\)，\(\gcd(a,b)=\gcd(a,a-b)=\gcd(b,a-b)\)。

我们发现，其实第一步可以省略，直接执行第二步即可。

时间复杂度与辗转相除法一样：\(O(\log \max(a,b))\)。

代码

代码实现比较简单：

inline int gcd(int a , int b) {
	while(a != b) {
		if(a > b)
			a = a - b;
		else
			b = b - a;
	}
	return a;
}

更相减损法的证明

第一个公式十分好证：

因为 \(2a\) 与 \(2b\) 均为 \(2\) 的倍数，所以 \(2\) 必是 \(2a\) 与 \(2b\) 的一个公因数，所以最大公因数中包含 \(2\)。故可以将其提出。

接下来我们来证明第二个公式：

设 \(d=\gcd(a,b)\)。

则 \(a\) 与 \(b\) 必能表示成：\(a=k_1d\)，\(b=k_2d\)。

令 \(a > b\)，所以 \(a-b=d(k1-k1)\)，可以得到 \(d\mid (a-b)\)。

由此可得：\(\gcd(a,b)=\gcd(a,a-b)=\gcd(b,a-b)\)。
故得证。

二进制算法

对于处理数位较少的数的最大公因数，使用辗转相除法就可以了，但当数位达到上百位时便需要用到高精度。此时，辗转相除法就不那么适用了。对于这类问题，我们可以用支持高精度的二进制算法。

算法流程如下：

若 \(a=b\)，则 \(\gcd(a,b)=a\)；

否则，分情况讨论：

1. 若 \(a\) 和 \(b\) 均为偶数，则：\(\gcd(a,b)=2\times\gcd(a\div 2,b\div 2)\)。

2. 若 \(a\) 为偶数，\(b\) 为奇数，则 \(\gcd(a,b)=\gcd(a\div 2,b)\)。

3. 若 \(a\) 与 \(b\) 均为奇数，则 \(\gcd(a,b)=\gcd(a-b,b)\)。

最终答案即是第一个操作中被约掉的 \(2\) 乘操作 \(3\) 的结果。

其实二进制算法和更相减损法十分相似。这里就不给证明了。

代码

这里给出一个完整的高精度最大公因数的代码：

#include<iostream>
#include<iomanip>
#include<cstring>
#include<cmath>
#include<stack>
#include<queue>
#define int long long
using namespace std;
string s1 , s2;
int a[3005] , b[3005] , c[3005];
int len1 , len2;
 
inline void div1(int a[]) {
    int r = 0;
    for(register int i = len1;i >= 1;i --) {
        a[i] = a[i] + r * 10;
        r = a[i] % 2;
        a[i] /= 2;
    }
    while(!a[len1] && len1 > 1)
        len1 --;
    return;
}
 
inline void div2(int b[]) {
    int r = 0;
    for(register int i = len2;i >= 1;i --) {
        b[i] = b[i] + r * 10;
        r = b[i] % 2;
        b[i] /= 2;
    }
    while(!b[len2] && len2 > 1)
        len2 --;
    return;
}
 
inline bool equal(int a[] , int b[]) {
    if(len1 != len2)
        return false;
    for(register int i = 1;i <= len1;i ++)
        if(a[i] != b[i])
            return false;
    return true;
}
 
inline void poww(int a[] , int x) {
    for(register int i = 1;i <= 3000;i ++)
        a[i] *= x;
    for(register int i = 1;i <= 3000;i ++) {
        a[i + 1] += a[i] / 10;
        a[i] %= 10;
    }
    len1 = 3000;
    while(!a[len1] && len1 > 1)
        len1 --;
    for(register int i = len1;i >= 1;i --)
        cout << a[i];
    return;
}
 
inline bool bigger(int a[] , int b[]) {
    if(len2 < len1)
        return false;
    if(len2 > len1)
        return true; 
    for(register int i = len1;i >= 1;i --)
        if(b[i] > a[i])
            return true;
        else if(a[i] > b[i])
            return false;
    return false;
}
 
inline void change(int a[] , int b[]) {
    for(register int i = 1;i <= len1;i ++)
        c[i] = a[i];
    for(register int i = 1;i <= len2;i ++)
        a[i] = b[i];
    for(register int i = 1;i <= len1;i ++)
        b[i] = c[i];
    int k = len1;
    len1 = len2;
    len2 = k;
    for(register int i = 3000;i > len2;i --)
        b[i] = 0;
    for(register int i = 3000;i > len1;i --)
        a[i] = 0;
    return;
}
 
inline void mi(int a[] , int b[]) {
    for(register int i = 1;i <= len1;i ++)
        a[i] = a[i] - b[i];
    for(register int i = 1;i <= len1;i ++)
        if(a[i] < 0) {
            a[i] += 10;
            a[i + 1] --;
        }
    while(!a[len1] && len1 > 1)
        len1 --;
    return;
}
 
inline int ksm(int x) {
    int result = 1 , base = 2;
    while(x > 0) {
        if(x & 1)
            result = result * base;
        x >>= 1;
        base = (base * base);
    }
    return result;
}
//精髓部分
inline void gcd(int a[] , int b[]) {
    int i , j;
    for(i = 0;a[1] & 0;i ++)
        div1(a);
    for(j = 0;b[1] & 0;j ++)
        div2(b);
    if(i > j)
        i = j;
    int k = ksm(i);
    while(true) {
        if(equal(a , b)) {
            poww(a , k);
            exit(0);
        }
        if(bigger(a , b))
            change(a , b);
        mi(a , b);
        while(a[1] & 0)
            div1(a);
    }
}
 
signed main() {
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);
    cin >> s1 >> s2;
    len1 = s1.size();
    len2 = s2.size();
    for(register int i = 0;i < len1;i ++)
        a[len1 - i] = s1[i] - '0';
    for(register int i = 0;i < len2;i ++)
        b[len2 - i] = s2[i] - '0';
    gcd(a , b);
    return 0;
}

c++ 自带函数

在这么多求最大公因数的算法中，应用的比较的多的就只有辗转相除法。但其实 c++ 内部也有自带的求最大公因数的函数，即：__gcd() 函数，用法和手写的一样。若要求 \(a\) 与 \(b\) 的最大公因数，代码应为：cout << __gcd(a,b)。

最小公倍数

定义

两数的公倍数中最小的那个数被称作这两个数的最小公倍数。

我们一般用 \(LCM(a,b)\) 表示两个数的最小公倍数。

最小公倍数的求法

易得：\(LCM(a,b)=a\div \gcd(a,b)\times b\)。

这个公式十分容易得到，这里就不给推导过程了。

个人建议先除以 \(\gcd(a,b)\) 再乘 \(b\)，以防止 \(a\) 和 \(b\) 过大，然后相乘溢出。

代码

int gcd(int a , int b) {
	if(!b)
		return a;
	return gcd(b , a % b);
}

inline int lcm(int a , int b) {
	return a / gcd(a , b) * b;
}

快速幂与扩展欧拉定理

把这两个知识点放一起是因为扩展欧拉定理需要用到快速幂。

快速幂

定义

快速幂，即很快的做幂运算，并且在求幂的同时支持取模操作。但这个算法究竟是怎么实现的，我们来探究一下。

快速幂的原理

假设我们需要求 \(2^8\) 是多少？

最暴力的方法就是一个一个乘：

\[\underbrace{2\times 2\times \cdots \times 2\times 2}_{8 个 2}=256 \]

这样的时间复杂度是：\(O(n)\)。当指数大时，这样的暴力一定 T 到飞起。

当然，聪明的你们（不是我）一定发现了一个简便的方法，那就是让底数不断平方，然后指数不断除以 \(2\)。即：

\[2^{8}=2^{2^4}=4^4=4^{2^2}=16^2=256 \]

不难发现，这样求幂只运算了 \(3\) 次，是 \(\log\) 级别的。但如果指数不是 \(2\) 的倍数怎么办呢？我们再来看一个例子：

求 \(2^{10}\) 是多少？

我们还是先按之前的方法做：

\[2^{10}=2^{2^5}=4^5 \]

现在我们发现指数变成了奇数，不能再除以 \(2\) 了。对于这种情况，我们可以这样处理：把一个底数单独拿出来，使得指数减一而成为偶数。即：

\[4^5=4^4\times 4 \]

经过处理，我们就可以继续重复原来的步骤了：

\[4^4\times 4=16^2\times 4=256\times 4=1024 \]

这两个步骤就是快速幂算法的主要流程。

具体实现：

我们可以定义变量 result，记得一定要初始化为 \(1\)。然后如果指数变成了奇数，则令result乘上底数，否则指数除以二，底数平方。

代码实现

由上面的两个步骤，我们可以写出基本的算法框架：

int ksm(int base , int x) {
	int result = 1;
	while(x) {
		if(x % 2 == 1)
			result *= base;
		x /= 2;
		base *= base;
	}
	return result;
}

但这还不是最快的算法。我们可以将取模运算以及除法运算改成位运算：

inline int ksm(int base , int x) {
	int result = 1;
	while(x) {
		if(x & 1)
			result *= base;
		x >>= 1;
		base *= base;
	}
	return;
}

欧拉定理与扩展欧拉定理

欧拉定理

定义

对于 \(a\)，\(m\) 两数，若 \(\gcd(a,m)=1\) 则有：\(a^{\varphi(m)}\equiv 1\pmod m\)。

其实费马小定理就是欧拉定理的特殊情况，即当 \(m\) 为素数时。

欧拉定理的证明

设 \(A=\left\{x_1,x_2,\cdots,x_{\varphi(m)}\right\}\) 是 \(1-m\) 中与 \(m\) 互素的数的集合。

则这些数模 \(m\) 互不相同，且余数与 \(m\) 互素。

接下来，我们来证明集合 \(B=\left\{ax_1,ax_2,\cdots,ax_{\varphi(m)}\right\}\) 且 \(\gcd(a,m)=1\) 也有这样的性质。

我们先证这些数模 \(m\) 互不相同：

考虑反证法：

若 \(ax_i\equiv ax_j\pmod m\) 且 \(i\neq j\)。

则有 \(ax_i-ax_j\equiv 0\pmod m\)。

即 \(a(x_i-x_j)\equiv 0\pmod m\)。

因为 \(a\) 与 \(m\) 互素且 \(x_i-x_j\nmid m\) 所以这种情况不存在，得证。

接下来，我们证余数与 \(m\) 互素：

因为 \(a\) 与 \(m\) 互素，且 \(x_i\) 与 \(m\) 互素。

所以 \(ax_i\) 与 \(m\) 互素。得证。

所以 \(B\) 集合在模 \(m\) 后便于 \(A\) 集合相等。可得：

\[\prod_{i=1}^{\varphi(m)} ax_i\equiv \prod_{i=1}^{\varphi(m)} x_i\pmod m \]

把 \(a\) 提出得：

\[a^{\varphi(m)}\prod_{i=1}^{\varphi(m)} x_i\equiv\prod_{i=1}^{\varphi(m)} x_i\pmod m \]

同时约掉 \(\prod_{i=1}^{\varphi(m)} x_i\) 得：

\[a^{\varphi(m)}\equiv 1\pmod m \]

得证。

欧拉定理的推论与证明

若 \(a\)，\(m\) 互素，则有 \(a^b\equiv a^{b\bmod \varphi(m)}\pmod m\)。

下面给出证明：

设 \(b=k\times\varphi(m)+t\)，其中 \(t=b\bmod \varphi(m)\)。

把 \(b\) 代入 \(a^b\) 得：

\[\begin{aligned} a^b &= a^{k\times\varphi(m)+t}\\ &= a^{k\times \varphi(m)}\times a^t\\ &=a^{\varphi(m)^k}\times a^t\\ &\equiv 1^k\times a^{b\bmod \varphi(m)}\pmod m\\ &\equiv a^{b\bmod \varphi(m)}\pmod m\\ \end{aligned} \]

得证。

扩展欧拉定理

在大部分情况下，我们都使用扩展欧拉定理，因为它可以处理的情况更多。

定义

对于 \(a^b\bmod m\) 这样一个式子，无论 \(a\) 与 \(m\) 互不互素都有：

\[a^b\bmod m\equiv \begin{cases} b\ge \varphi(m),a^{b\bmod \varphi(m)+\varphi(m)}\pmod m\\ b<\varphi(m),a^b\pmod m \end{cases} \]

• 对于第一种情况，即 \(b\ge\varphi(m)\) 时我们可以用公式化简来求解。

• 对于第二种情况，我们可以直接使用快速幂算出结果。

算法流程十分简单，我们只需要判断 \(b\) 与 \(\varphi(m)\) 的大小，然后按照对应的公式计算即可。

代码

inline int phi(int x) {
	int ans = x;
	for(register int i = 2;i <= sqrt(x);i ++)
		if(x % i == 0) {
			ans = ans / i * (i - 1);
			while(x % i == 0)
				x /= i;
		}
	if(x > 0)
		ans = ans / x * (x - 1);
	return ans;
}

inline int ksm(int base , int x , int m) {
	int result = 1;
	while(x) {
		if(x & 1)
			result = result * base % m;
		x >>= 1;
		base = base * base % m;
	}
	return result;
}

inline int solve(int a , int b , int m) {
	int phi_m = phi(m);
	if(b >= phi_m)
		return ksm(a , b % phi_m + phi_m , m);
	else
		return ksm(a , b , m);
}

证明不放了，以后再 update 吧。

扩展欧几里得与同余

线性同余方程需要用到扩展欧几里得来求解。

同余

简单讲一下同余。

定义

若整数 \(a\)，\(b\) 除以 \(m\) 的余数相等，则称 \(a\) 与 \(b\) 模 \(m\) 同余，记作：

\[a\equiv b\pmod m \]

性质

余数的性质：

• 可加性。
• 可减性。
• 可乘性。
• 可乘方性。

这几个性质光看应该就能理解了，就不过多赘述了。

注意，余数没有可除性！

扩展欧几里得

我们先看一个定理。

裴蜀定理

定义

对于 \(ax+by=c\) 这样一个方程，设 \(d=\gcd(a,b)\)。则只有在 \(d\mid c\) 的情况下，这个方程才有整数解。

感性理解一下，不放证明了。

这个定理可以很好的帮助我们判断线性方程是否有整数解，非常有用。

扩展欧几里得算法

扩展欧几里得的特解

对于 \(ax+by=\gcd(a,b)\) 这样一个方程，我们可以用扩展欧几里得算法来求出一组特解，即：

\[\begin{cases} x=y_1\\ y=x_1-a\div b\times y_1 \end {cases} \]

这是怎么得到得呢？接下来给出推导过程。

扩展欧几里得的推导过程

同样的，我们可以得到另外一个方程：

\[bx_1+(a\bmod b)y_1=\gcd(b,a\bmod b) \]

把 \(a\bmod b\) 化简得：

\[bx_1+(a-a\div b\times b)y_1=\gcd(b,a\bmod b) \]

整理一下等式左边得：

\[\begin{aligned} bx_1+(a-a\div b\times b)y_1 &= bx_1+ay_1-(a\div b\times b)y_1\\ &=ay_1+b(x_1-a\div b\times y_1) \end{aligned} \]

在 \(b\neq 0\) 的情况下，有：\(\gcd(a,b)=\gcd(b,a\bmod b)\)。

将 \(\gcd(b,a\bmod b)\) 替换得：

\[\gcd(a,b)=ay_1+b(x_1-a\div b\times y_1) \]

同时，我们还有：

\[\gcd(a,b)=ax+by \]

对比两式得到：

\[\begin{cases} x=y_1\\ y=x_1-a\div b\times y_1 \end {cases} \]

推导完毕。

扩展欧几里得的算法流程

得到了特解，我们该怎么求这个特解呢？

其实这是一个递归的过程。

不难发现，在不断求解的过程中，等式右边的 \(b\) 终将变为 \(0\)，此时我们就可以得到：

\[\begin{cases} x=1\\ y=0 \end {cases} \]

接着，我们便可以根据上面推导出来的式子不断递归回去而得到特解。

代码

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

扩展欧几里得的通解与最小正整数解

在得到了 \((x_0,y_0)\) 这一组特解后，我们可以得到这个方程的通解。其通解为：

\[\begin{cases} x=x_0+k\frac{b}{\gcd(a,b)} \\ y=y_0-k\frac{a}{\gcd(a,b)} \\ \end{cases} \]

我们还可以得到该方程的最小正整数解，即：

\[x_{min}=(x_0 \% \frac{b}{\gcd(a,b)}+\frac{b}{\gcd(a,b)})\%\frac{b}{\gcd(a,b)} \]

扩展欧几里得的应用

扩展欧几里得在许多地方都有所应用。

这里先讲解用其求解线性方程。

形如 \(ax+by=c\) 这样的方程，都可以用扩展欧几里得来求解。怎么求呢？

首先，由裴蜀定理，在 \(\gcd(a,b) \nmid c\) 时无整数解。

然后，容易想到先求出 \(ax+by=\gcd(a,b)\) 的一组解。

设求得的解为 \(x_0\) 与 \(y_0\)。

有：

\[ax_0+by_0=\gcd(a,b) \]

此时，我们可以在等式两边同时乘上 \(\frac{c}{\gcd(a,b)}\)，得到原方程的一组特解：

\[\begin{cases} x=x_0\times \frac{c}{\gcd(a,b)} \\ y=y_0\times \frac{c}{\gcd(a,b)} \\ \end{cases} \]

若题目要求最小正整数解，则按上面的公式算一遍即可。

扩展欧几里得还能用来求解一些同余方程。

比如若要求 \(ax\equiv c\pmod b\) 的解。

为了求解，我们可以对这个方程进行一些变换：

原方程其实等价于：

\[ax=-by+c \]

为什么 \(y\) 是负的？因为 \(y\) 可以取任意整数，所以正负无所谓了。

对上式进一步简化，得到：

\[ax+by=c \]

此时便可以用扩展欧几里得求解了。

逆元

定义

若整数 \(a\)，\(x\) 满足 \(ax\equiv 1\pmod p\)，则称 \(x\) 为 \(a\bmod p\) 的逆元，记作 \(a^{-1}\)。逆元存在充要条件是 \(\gcd(a,p)=1\)。

本文中讨论的逆元为模意义下的乘法逆元，即 \(x\in \left\{1,2,\cdots p-1\right\}\)。

逆元的求法

本文主要介绍逆元的两种求法。

费马小定理求逆元

在上文中提到的费马小定理就可以用来求解特定情况下数的逆元。

费马小定理的式子是：

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

现在我们对等式右侧进行变化，可得：

\[a\times a^{p-2}\equiv 1\pmod p \]

将得到的式子与逆元的定义式对比：

\[\begin{cases} a\times a^{p-2}\equiv 1\pmod p\\ a\times x\equiv 1\pmod p\\ \end{cases} \]

可得 \(a^{p-2}\) 即为 \(a\bmod p\) 的逆元。

由于计算 \(a^{p-2}\) 涉及到幂运算，所以我们可以用快速幂优化时间，并同时取模。

注意：该算法只能在 \(p\) 为素数时使用。

代码

此处默认 \(p\) 为素数。

int a , p;

inline int ksm(int base , int x , int mod) {
	int result = 1;
	while(x) {
		if(x & 1)
			result = (result * base) % mod;
		x >>= 1;
		base = (base * base) % mod;
	}
	return result;
}

inline int solve(int a , int p) {
	return ksm(a , p - 2 , p);
}

扩展欧几里得算法求逆元

不难发现，逆元的定义式其实就是一个线性同余方程，那自然可以用扩展欧几里得来求解。

具体推导过程和讲解扩展欧几里得算法时的类似，就不再写一次了。这里直接给出最终的等式：

\[ax+py=1 \]

得出等式后，直接求解即可。不过该算法求出来的 \(x\) 可能不在 \(\left\{1,2,\cdots,p-1\right\}\) 中，所以还要做下面这步运算，即：\(x=(x\% \frac{p}{\gcd(a,p)}+\frac{p}{\gcd(a,p)})\% \frac{p}{\gcd(a,p)}\)。

代码

int a , p;

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

inline int solve(int a , int p) {
	int x , y;
	int r = exgcd(a , p , x , y);
	int k = p / r;
	x = (x % k + k) % k;
	return x;
}

逆元的作用

众所周知，模运算符合加法、减法和乘法的性质，但并不满足除法的性质。这也就意味着：\((a\div b)\bmod p \neq (a \bmod p)\div (b\bmod p)\)。但当数据范围过大时，不在运算中取模会使数据溢出，所以取模又是不可避免的。此时，逆元便能在除法的模运算中化除为乘，让代码在取模后仍能得到正确答案。

什么意思呢？即：

\[(a\div b)\equiv (a\times b^{-1})\pmod p \]

其中 \(b^{-1}\) 为 \(b\bmod p\) 的逆元。

这是怎么来的，下面给出推导过程：

\[\because b\times b^{-1}\equiv 1\pmod p \]

\[\therefore (a\div b)\equiv (a\div b\times b\times b^{-1})\pmod p \]

\[\therefore (a\div b)\equiv (a\times b^{-1})\pmod p \]

推导完毕。

中国剩余定理

就详见这里了。

后记

写了一些，但 BSGS 等还没写，以后在 update。