进阶数论
TonyYin
2020-02-11 09:50:33
# 进阶数论
# By YYWD
## www.bjeaedu.com/
部分资料来自www.jisuanke.com
## 扩展欧几里得
### 概述
扩展欧几里得算法是用来在已知$a,b$的情况下求解一组$x,y$，使满足$ax+by=gcd(a,\,b)=d$.（$gcd$为最大公因数，方程一定有解）
### 特解
要计算$gcd(a,\,b)$，并求出满足条件的$x,y$.
将下一个状态 $b$ 和 $(a\%b)$ 表示为：$b\cdot x_1+(a\%b)\cdot y_1=d$
则：
$$
\because b\cdot x_1+(a\%b)\cdot y_1=d,
$$
又
$$\because a\%b=a-\left\lfloor\frac{a}{b}\right\rfloor\times b
$$
$$
\begin{aligned}
\therefore gcd(a,\,b)=d & =b\cdot x_1+(a-\left\lfloor\frac{a}{b}\right\rfloor\times b)\cdot y_1\\
&=b\cdot x_1+a\cdot y_1-\left\lfloor\frac{a}{b}\right\rfloor\times b\times y_1\\
&=a\times y1+b\times(x_1-\left\lfloor\frac{a}{b}\right\rfloor\times y_1)
\end{aligned}
$$
所以满足条件的$x=y_1,\,y=x_1-\left\lfloor\frac{a}{b}\right\rfloor\times y_1$
~~~cpp
int exgcd(int a, int b, int &x, int &y) {//函数返回值为gcd(a,b)
if(b == 0) {
x = 1;//gcd(a,b) = a = 1a+0b
y = 0;
return a;
}
int r = exgcd(b, a % b, x, y);//先递归
int t = x;//再计算x,y
x = y;
y = t - a / b * y;
return r;
}
~~~
### 通解
假设特解为$(x_0,y_0)$，满足$ax_0+by_0=gcd(a,\,b)=d$
则有
$$
\begin{aligned}
a(x_0+k\frac{b}{d})+b(y_0-k\frac{a}{d})=d
\end{aligned}
$$
(k是整数)
所以方程的通解为：
$$
\begin{aligned}
x=x_0+k\frac{b}{d} \\ y=y_0-k\frac{a}{d}
\end{aligned}
$$
对于其它二元一次不定方程$ax+by=c$（$c$为任意正整数）只有当$d\mid c$时才有整数解，通解形式：
$$
\begin{aligned}
x=\frac{c}{d}x_0+k\frac{b}{d}\\y=\frac{c}{d}y_0-k\frac{a}{d}
\end{aligned}
$$
## 同余方程（模方程）
### 定义：
指形如$ax\equiv b\pmod c$的一种方程。其中$a,b,c$是参数，$a,c$是正整数，$b$是小于$c$的非负整数，$x$是未知数。
### 扩展欧几里得解同余方程
根据同余定义，$ax\equiv b\pmod c$ 可化为 $xa+kc=b$. 其中$a,b,c$已知，$x,k$未知，所以能用扩展o欧几里得解同余方程。
### 通解
如果$b\nmid gcd(a,c)$，方程无整数解.
否则同扩展欧几里得酸法，求得通解：
$$
x=\frac{b}{d}x_0+k\frac{c}{d}\\k=\frac{b}{d}p_0-k\frac{a}{d}
$$
其中k为整数
### 最小非负整数解
求$x$的最小非负整数解
因为 $x=\frac{b}{d}x_0+k\frac{c}{d}$，其中$x_0$是方程$xa+kc=gcd(a,\,c)$的特解， $\frac{b}{d}x_0$ 是原同余方程的特解，
若$x>0$，则有$x\equiv \frac{b}{d}x_0 \pmod{ \frac{c}{d}} $。最小非负整数解为$\frac{b}{d}x_0\,\%\,\frac{c}{d}$.
若$x<0$，同理，最小非负整数解为$(\frac{b}{d}x_0\bmod \frac{c}{d}+\frac{c}{d} )\bmod\frac{c}{d}$
综上，最小非负整数解=$(\frac{b}{d}x_0\bmod\frac{c}{d}+\frac{c}{d})\bmod \frac{c}{d}$
~~~cpp
int exgcd(int a, int b, int &x, int &y) {
if(b == 0) {
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;
}
int main() {
int a,b,c,x,p;
cin>>a>>b>>c;
int d=exgcd(a,c,x,p);
if(b%d!=0){
cout<<"No Solution"<<endl;
}
else{
cout<<(b/d*x%(c/d)+c/d)%(c/d)<<endl;
}
return 0;
}
~~~
## 逆元
### 乘法逆元
给定两个整数$a$和$p$，假设存在$x$使得$ax\equiv 1\pmod p)$.称 $x$ 为 $a$ 关于 $p$ 的**乘法逆元**。$a$关于$p$的逆元存在的充分必要条件是$gcd(a,p)=1$
### 求一个数的逆元
将上面式子变形，变成$ax+kp=1$，用扩展欧几里得能够求出一个$x$.
如果$p$是质数，有另一种方法求$a$关于$p$的逆元：
$$
a^{p-1}\bmod p=(a^{p-2}\times a)\bmod p=(x\times a)\bmod p=1
$$
所以$a$关于$p$的逆元是 $a^{b-2}\bmod b$ 的值，用快速幂优化。
~~~cpp
int exgcd(int a, int b, int &x, int &y) {
if(b == 0) {
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;
}
int pow_mod(int x, int n, int mod) {
int res = 1, temp = x;
for (; n; n /= 2) {
if (n & 1) {
res = res * temp % mod;
}
temp = temp * temp % mod;
}
return res;
}
bool is_prime(int x) {
if (x == 1) {
return false;
}
for (int i = 2; i * i <= x; i++) {
if (x % i == 0) {
return false;
}
}
return true;
}
int main() {
int a,b;
cin>>a>>b;
int x,y;
int d=exgcd(a,b,x,y);
if(d!=1){
cout<<"No Solution"<<endl;
}
else{
cout<<(x+b)%b<<endl; //方法一
}
if(is_prime(b)){
cout<<pow_mod(a,b-2,b)<<endl;//方法2（p为质数时才能用）
}
return 0;
}
~~~
### 线性预处理逆元
有些情况，需要用到区间内的逆元。对于$\leq p-1$的正整数$n$，可以用$\mathcal{O}(n)$复杂度求解$1-n$的逆元。（用 $inv_i$表示 $i$ 关于 $p$ 的逆元）
方法：
$$
inv_i=(p-\left\lfloor\frac{p}{i}\right\rfloor)inv_{p\%i}\%p
$$
证明：要证$i\cdot inv_i\%p=1$
计蒜客方法：
$$
\begin{aligned}
i\cdot inv_i\%p&=i\times (p-\frac{p}{i})inv_{p\%i}\%p\\
&=i\times(p-\frac{p-p\%i}{i})inv_{p\%i}\%p\\
&=p\%i\times inv_{p\%i}\%p\\
&=1
\end{aligned}
$$
证明方法：
设$p=ki+r$，
则有：
$$
\begin{aligned}
ki+r&\equiv 0\pmod p\quad\quad\quad (2)\\
(1)\times inv_i\times inv_r&=(2),\\
k\cdot inv_r+inv_i&\equiv0\pmod p\quad\quad\quad(2)\\
\because p&=ki+r\\
\therefore k=\left\lfloor\frac{p}{i}\right\rfloor&,r=p\%i\\
\therefore \left\lfloor\frac{p}{i}\right\rfloor inv_{p\%i}+inv_i&\equiv0\pmod p\\
\therefore -\left\lfloor\frac{p}{i}\right\rfloor inv_{p\%i}&\equiv inv_i\pmod p\\
\therefore (p-\left\lfloor\frac{p}{i}\right\rfloor) inv_{p\%i}&\equiv inv_i\pmod p
\end{aligned}
$$
实现：
~~~cpp
// p 必须为质数，p / i 为整除。
inv[1] = 1;
for (int i = 2; i <= n; ++i) {
inv[i] = (p - p / i) * inv[p % i] % p;
}
~~~
### 组合数取模
设$invs_i=inv_1\times inv_2\times \cdots \times inv_i\%p,fact_i=1\times 2\times 3\times \cdots\times i\%p$
$C_n^m\%p=\frac{n!}{m!(n-m)!}\%p=fact_n\cdot invs_m\cdot invs_{n-m}\%p$
实现：
~~~cpp
fac[0]=1;
for(int i=1;i<=n;i++){
fac[i]=fac[i-1]*i%p;
}
inv[1]=1;
for(int i=2;i<=n;i++){
inv[i]=(p-p/i)*inv[p%i]%p;
}
inv[0]=1;
for(int i=1;i<=n;i++){
inv[i]=inv[i]*inv[i-1]%p;
}
cout<<fac[n]*inv[m]*inv[n-m]%p<<endl;
~~~
## 中国剩余定理
### 内容
中国剩余定理一元同余线性方程组：
$$
\begin{aligned}&x\equiv a_1\pmod {m_1}\\&x\equiv a_2\pmod {m_2}\\&\cdots\\&x\equiv a_n\pmod {m_n}\\\end{aligned}
$$
假设整数$m_1,m_2,\cdots,m_n$**两两互质**，则方程组有解，通解如下：
设$M=m_1\times m_2\times \cdots\times m_n=\prod m_i$，并设$M_i=M/m_i$， $t_i=M_i^{-1}$.
$t_i$表示$M_i$模$m_i$意义下的逆元，即$M_it_i\equiv1\pmod m_i$.
通解为：
$$
x=a_1t_1M_1+a_2t_2M_2+\cdots+a_nt_nM_n+jM=\sum_{i-1}^{n}a_it_iM_i+kM
$$
### 实现
程序输出的是最小非负整数解
~~~cpp
#include <iostream>
using namespace std;
int exgcd(int a, int b, int &x, int &y) {
if (b == 0) {
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;
}
int a[101], m[101];
int CRT(int n){
int M=1;
int ans=0;
for(int i=1;i<=n;i++){
M*=m[i];
}
for(int i=1;i<=n;i++){
int x,y;
int Mi=M/m[i];
exgcd(Mi,m[i],x,y);
ans=(ans+Mi*x*a[i])%M;
}
if(ans<0){//最小非负整数解就是在模 M 意义下进行求解，
//如果答案是负数就加 M 变成非负数，最后输出答案
ans+=M;
}
return ans;
}
int main() {
int n;
cin>>n;
for(int i=1;i<=n;i++){
cin>>m[i]>>a[i];
}
cout<<CRT(n)<<endl;
return 0;
}
~~~
## 一般同余线性方程组
### 概述
由若干个同余方程组成的方程组叫做同余方程组$\ldots\ldots$<font size=2>废话</font>
与中国剩余定理的区别：模数不一定互质
### 算法
考虑两个同余方程：
$$
\begin{aligned}x\equiv a_1(mod\;m_1)\\x\equiv a_2(mod\;m_2)\end{aligned}
$$
可以写为：
$$
\begin{aligned}x+k_1m_1=a_1\quad\quad\quad(1)\\ x+k_2m_2=a_2\quad\quad\quad(2)\end{aligned}
$$
$(1)-(2)$得到：
$$
m_1k_1-m_2k_2=a_1-a_2
$$
其中$m1,m2,a1,a2$已知，发现方程形如$ax+by=c$.
所以方程有解条件是$gcd(m1,m2)\mid (a_1-a_2)$
可以用扩展欧几里得算法求得$k_1$的最小非负整数解$K$，通解为：
$$
k_1=K+\frac{m2}{gcd(m_1,m_2)}\times C
$$
其中$C$为整数.
设$d=gcd(m_1,m_2)$，对方程两边除以$d$，得到：
$$
\frac{m_1k_1}{d}=\frac{a_1-a_2}{d}+\frac{m_2}{d}\times k_2
$$
写成同余方程的形式：
$$
\frac{k_1m_1}{d}\equiv \frac{a_1-a_2}{d}(mod\;\frac{m_2}{d})
$$
同余方程两边同乘$d$：
$$
k_1m_1\equiv a_1-a_2(mod\;\frac{m_2}{d})
$$
将$k_1$的通解$k_1=K+\frac{m2}{d}\times C$代回方程$x+k_1m_1=a_1$，得到：
$$
x=a_1-Km_1-\frac{m_1m_2C}{d}
$$
写成同余方程的形式：
$$
x\equiv a_1-Km_1(mod\;\frac{m_1m_2}{d})
$$
这样把两个同余方程合并成一个。同理可以将$n$个方程组成的同余方程组化为一个同余方程，最后用扩展欧几里得算法求解$x$。
### 实现
~~~cpp
#include <iostream>
using namespace std;
int exgcd(int a, int b, int &x, int &y) {//ax+by=d;
if (b == 0) {
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;
}
int a[101], m[101];
int exCRT(int a[],int m[],int n){
int M=m[1],A=a[1],x,y,d;
//a是余数数组，m是模数数组，n是方程数量
//M记录合并后方程的模数，A是合并后方程的余数x是M特解,y是mi特解
for(int i=2;i<=n;i++){
d=exgcd(M,m[i],x,y);
if((A-a[i])%d!=0){
return -1;
}
x=((A-a[i])/d*x%(m[i]/d)+m[i]/d)%(m[i]/d);//最小非负整数解
A-=x*M;//新的余数就是原来的余数减去x倍原来的模数再对新的模数取余
M=M/d*m[i];//新的模数是原来的模数/d*m[i];
A=(A%M+M)%M;//防止A为负数
}
return A;
}
int main() {
int n;
cin>>n;
for(int i=1;i<=n;i++){
cin>>m[i]>>a[i];
}
cout<<exCRT(a,m,n)<<endl;
return 0;
}
~~~
## 练习题
| **算法** | **题号** 点击可跳转 |
| ------------------ | ------------------------------------------------------------ |
| 扩展中国剩余定理 | [$\color{#9d3dcf}{P4777}$](https://www.luogu.com.cn/problem/P4777) |
| 中国剩余定理 | [$\color{#53c41a}{P1495}$](https://www.luogu.com.cn/problem/P1495) |
| 矩阵快速幂 | [$\color{#ffc116}{P3390}$](https://www.luogu.com.cn/problem/P3390) |
| 区间预处理乘法逆元 | [$\color{#ffc116}{P3811}$](https://www.luogu.com.cn/problem/P3811) |
| 扩展欧几里得算法 | [$\color{#53c41a}{P1082}$](https://www.luogu.com.cn/problem/P1082) |
| 矩阵快速幂 | [$\color{#53c41a}{P1962}$](https://www.luogu.com.cn/problem/P1962) |
| 扩展欧几里得算法 | [$\color{#3498db}{P1516} $](https://www.luogu.com.cn/problem/P1516) |
| 中国剩余定理 | [$\color{#3498db}{P1495}$](https://www.luogu.com.cn/problem/P1495) |
### 子段乘积
[$\textcolor{red}{Click\;\;\; here}$](https://ac.nowcoder.com/acm/contest/3005/C)
给出一个长度为 n 的数列 $a_1,a_2,\ldots,a_n$ 求其长度为 $k$ 的连续子段的乘积对 $998244353$ 取模余数的最大值。