gcd，exgcd与逆元
gcd
最大公因数，也称最大公约数、最大公因子，指两个或多个整数共有约数中最大的一个。\(a\)，\(b\)的最大公约数记为\((a,b)\)，同样的，\(a\)，\(b\)，\(c\)的最大公约数记为\((a,b,c)\)，多个整数的最大公约数也有同样的记号。求最大公约数有多种方法，常见的有质因数分解法、短除法、辗转相除法、更相减损法。与最大公约数相对应的概念是最小公倍数，\(a\)，\(b\)的最小公倍数记为\([a,b]\)。 ----选自百度百科
也就是说，若有两个数\(a,b\)，若存在\(x\)使得\(x\mid b\)且\(x\mid a\)，则称\(x\)为\(a,b\)的公约数，\(a,b\)的所有公约数中，最大的即为最大公约数，记为\((a,b)\)或\(\gcd(a,b)\)。
（补充）相对应地，若存在\(x\)使得\(a\mid x\)且\(b\mid x\)，则称\(x\)为\(a,b\)的公倍数，\(a,b\)的所有公倍数中，最小的即为最小公倍数，记为\([a,b]\)或\(\operatorname{lcm}(a,b)\)。\(\operatorname{lcm}(a,b)=\dfrac{a*b}{\gcd(a,b)}\)。
求\(gcd\)的方法有质因数分解法、短除法、辗转相除法、更相减损法。下面粗略讲解前两个并详细讲第三个，其余自行右转百度百科。
- 质因数分解法
顾名思义，该方法就是将\(a\)和\(b\)都分解质因数，所有公共部分的积即为\(\gcd\)。
比如，若想要求\(\gcd(24,60)\)，首先将\(24\)和\(60\)质因数分解得到\(24=2^3\times 3\)，\(60=2^2\times 3\times 5\)，公共部分即为\(\gcd(24,60)=2^2\times 3=120\)。这种做法保证了得到的数\(x\)是\(a,b\)的公因数，因为\(x\)是从\(a,b\)的因数中选出若干数相乘的；同时确保\(x\)最大，因为剩余部分不存在\(a,b\)的公因数了。 - 短除法
这种做法和上面的做法大体相同，我们将要求\(\gcd\)的数并排写，逐次除以公约数直到商互质，这时所有除过的数的积即为\(\gcd(a,b)\)。如下图可以求出\(\gcd(24,60)\)。
欧几里得算法求gcd
欧几里得算法通常称为辗转相除法。
算法内容
\(\gcd(a,b) = \gcd(b,a\bmod b)\)，故可一直递归直到\(b = 0\)时，此时\(a\)为原来两数的最大公约数。如下面给出求\(\gcd(24,60)\)的过程（其实就是不断使用上述等式）：
证明
令 \(\gcd(a,b) = x\)，则设\(a=px,b=qx\)（显然\(\gcd(p,q)=1\)，即\(p,q\)互质）
再令 \(r = a\bmod b , k=\left\lfloor\dfrac{a}{b}\right\rfloor\)，则\(r=a-kb=px-k\cdot qx=(p-kq)x\)
\(\therefore \gcd(b,a\bmod b)=\gcd(b,r)=\gcd(qx,(p-kq)x)=\gcd(q,p-kq)\cdot x\)
令\(\gcd(q,p-kq)=t\)，又设\(q=p't,p-kq=q't\)（\(\gcd(p',q')=1\)）
则\(p=q't+kq=q't+k\cdot p't=(q'+kp')t\)
\(\because a=px=(q'+kp')tx,b=qx=p'tx\)
\(\therefore \gcd(a,b)=\gcd((q'+kp')tx,p'\cdot tx)=\gcd(q'+kp',p')tx\)
\(\because kp'\)为 \(p'\) 的倍数而 \(q'\) 不是\(\ \ \ \ \ \therefore \gcd(q'+kp',p')=1\)
\(\therefore \gcd(a,b)=\gcd(q'+kp',p')tx=tx\)
又 \(\because \gcd(a,b) = x\ \ \ \ \ \therefore t=\gcd(q,p-kq)=1\)
\(\therefore \gcd(b,a\bmod b)=\gcd(q,p-kq)\cdot x=x=\gcd(a,b)\)，即\(\gcd(a,b)=\gcd(b,a\bmod b)\)
代码
#include<iostream>
#include<cstdio>
using namespace std;
int a,b;
int gcd(int x,int y){
return y==0?x:gcd(y,x%y);
}
int main(){
scanf("%d%d",&a,&b);
printf("%d",gcd(a,b));
return 0;
}
exgcd
\(exgcd\)，即扩展欧几里得算法，也可以叫做扩欧。
扩展欧几里得算法是欧几里得算法的扩展。除了计算\(a,b\)两个整数的最大公约数，此算法还能找到整数\(x,y\)（其中一个很可能是负数）使得\(ax + by = \gcd(a,b)\)。有两个数\(a,b\)，对它们进行辗转相除法，可得它们的最大公约数——这是众所周知的。然后，收集辗转相除法中产生的式子，倒回去，可以得到\(ax+by=\gcd(a,b)\)的整数解。 ----节选自百度百科
让我们看一道例题：
例：P1082 [NOIP2012 提高组] 同余方程
题目意思是给定两个数\(a,b\)，求同余方程\(ax\equiv1\pmod{b}\)的最小正整数解。
考虑到暴力枚举\(x\)一定会超时，我们考虑用\(exgcd\)求解。
把问题转化成\(ax+by=1(y < 0)\)，那么就变成了我们刚刚已知的\(ax + by = \gcd(a,b)\)形式了。
这里多证明一下方程\(ax+by=c\)有解则满足\(c\bmod \gcd(a,b)=0\)，故此题中\(1\bmod \gcd(a,b)\)一定等于\(0\)，即\(\gcd(a,b)=1\)
令\(a=m\gcd(a,b),b=n\gcd(a,b)\)，则一定有\(c=ax+by=(mx+ny)\gcd(a,b)\)，即\(c\bmod \gcd(a,b)=0\)
证明
（下面是\(exgcd\)的推导过程，即方程\(ax+by=\gcd(a,b)\)的最小整数解而非\(ax+by=1\)）
我们已知一组\(a,b\)，想求\(ax + by = \gcd(a,b)\)中\(x\)的最小正整数解，不妨设\(G=\gcd(a,b)\)
我们设两个数\(x_1,y_1\)，令 \(bx_1 + (a\bmod b)y_1 = G\)，则\(ax+by=bx_1 + (a\bmod b)y_1\)
又 \(\because a\bmod b = a - b\cdot \left\lfloor\dfrac{a}{b}\right\rfloor\ \ \ \ \ \ \therefore ax+by=bx_1+(a - b\cdot \left\lfloor\dfrac{a}{b}\right\rfloor)\cdot y_1\)
\(\therefore ax+by=bx_1+ay_1-b\cdot \left\lfloor\dfrac{a}{b}\right\rfloor\cdot y_1=ay_1+b(x_1-\left\lfloor\dfrac{a}{b}\right\rfloor\cdot y_1)\)
所以必然存在一组解： \(\begin{cases}x=y_1\\y=x_1-\left\lfloor\dfrac{a}{b}\right\rfloor\cdot y_1\end{cases}\)
于是，问题变成了如何求\(x_1,y_1\)，此时我们需要求的是\(bx_1+(a\bmod b)y_1=G\)的一组正整数解，那么，我们可以类比刚才的解法，用\(x_2,y_2\)算\(x_1,y_1\)，再用\(x_3,y_3\)算\(x_2,y_2\)……
以此类推，我们发现方程\(x\)和\(y\)的系数从\(a,b\)变成\(a_1=b,b_1=a\bmod b\)再往下变，刚好就是我们用辗转相除法求\(\gcd(a,b)\)时的系数，所以，我们可以知道，到了最后一定会出现\(a_n=G,b_n=0\)，此时，方程为\(Gx_n+0\cdot y_n=G\)，于是必然存在解\(\begin{cases}x_n=1\\y_n=0\end{cases}\)，此时\(y_n\)取\(0\)，是因为我们递归求\(y_i\)，一旦我们的\(y_n\)取\(1\)，\(y_0\)将会非常大。
最后，因为我们求得的\(x\)不一定是最小正整数解，我们将其\(\bmod b\)得到最小正整数解，而且，由于\(x\)可能是负数，所以模完还应该\(+b\)再\(\bmod b\)。
于是，我们可以写出代码：
#include<iostream>
#include<cstdio>
using namespace std;
int a,b;
int x,y;
void exgcd(int a,int b){
if(b==0){
x=1;
y=0;//y最好取0
return;
}
exgcd(b,a%b);
int xx=x;
x=y;
y=xx-a/b*y;//x,y记录该层的答案
return;
}
int main(){
scanf("%d%d",&a,&b);
exgcd(a,b);
x=(x%b+b)%b;
printf("%d",x);
return 0;
}
逆元
定义
如果一个线性同余方程\(ax\equiv1\pmod{b}\)，则\(x\)称为\(a\)在模\(b\)意义下的逆元，记作\(a^{-1}\)。 ----选自\(OI\ Wiki\)
我们在一般的题目中，若使用$$((12\bmod 5)+(3\bmod 5))\bmod5=(12+3)\bmod 5=0$$或$$((12\bmod 5)\times(3\bmod 5))\bmod5=(12\times3)\bmod5=1$$是可行的，但是，若是除法则不能：$$((12\bmod 5)/(3\bmod 5))\bmod5=\dfrac{2}{3}\ne (12/3)\bmod 5=4$$
因此，我们引入逆元来解决这一问题。由定义，假设在模\(b\)意义下\(a\)的逆元为\(a^{-1}\)（即满足\(a*a^{-1}\equiv1\pmod b\)），则\(\dfrac{x}{a}\equiv x*a^{-1}\pmod1\)。证明也很简单，不妨假设\(x=k\times a\)，则\(\dfrac{x}{a}=\dfrac{ka}{a}=k\)，而\(x*a^{-1}=(k*a)*a^{-1}=k*(a*a^{-1})=k\)。
那么，如何求一个数在模\(b\)情况下的逆元呢？
例：【模板】乘法逆元
大意
给定\(n,p\)，求\(1\sim n\)在模\(p\)意义下的逆元。
解法及代码
- 费马小定理
由定理，若\(a,p\)互质，则有\(a^{p-1}\equiv1\pmod p\)，又因为\(a*a^{-1}\equiv1\pmod p\)，故\(a*a^{-1}\equiv a^{p-1}\pmod p\)，所以\(a^{-1}\equiv a^{p-2}\pmod p\)，即\(a^{-1} = a^{p-2}\bmod p\)。
#include<iostream>
#include<cstdio>
#define ll long long
using namespace std;
ll n,p;
ll qp(ll di,ll mi){
if(mi==0){
return 1;
}
ll ans=qp(di,mi/2);
ans=ans*ans%p;
if(mi&1){
ans=ans*di%p;
}
return ans%p;
}
int main(){
scanf("%lld%lld",&n,&p);
for(int i=1;i<=n;i++){
printf("%lld\n",qp(i,p-2));
}
return 0;
}
- 欧拉定理
中的数论定理
由定理，若\(a,p\)互质，\(a^{\varphi(p)}\equiv1\pmod p\)，特殊地，\(p\)是质数时\(\varphi(p)=p-1\)，故式子变成\(a^{p-1}\equiv1\pmod p\)，后面过程同上，代码同上（均\(TLE\)）。 - \(exgcd\)
因为求\(i\)在模\(p\)意义下的逆元相当是解同余方程\(a*a^{-1}\equiv1\pmod p\)的解，直接\(exgcd\)求\(a^{-1}\)即可。
#include<iostream>
#include<cstdio>
using namespace std;
int a,b;
int x,y;
void exgcd(int a,int b){
if(b==0){
x=1;
y=0;
return;
}
exgcd(b,a%b);
int xx=x;
x=y;
y=xx-a/b*y;
return;
}
int main(){
scanf("%d%d",&a,&b);
for(int i=1;i<=a;i++){
exgcd(i,b);
x=(x%b+b)%b;
printf("%d\n",x);
}
return 0;
}
理论可以AC，可是我的代码还是TLE了
4. 递推
设\(k=\left\lfloor\dfrac{p}{a}\right\rfloor,r=p\bmod a\)
因为\(p=k\times a+r\)，则\(k\times a+r\equiv0\pmod p\)
两边同乘\(a^{-1}*r^{-1}\)得到\(k*r^{-1}+a^{-1}\equiv0\pmod p\)，移项得到\(a^{-1}\equiv-k*r^{-1}\pmod p\)
重新代入\(k,r\)得到\(a^{-1}\equiv-\left\lfloor\dfrac{p}{a}\right\rfloor*(p\bmod a)^{-1}\pmod p\)
由于\(p\bmod a<a\)，所以算\(inv_a\)的时候已经算出来了\(inv_{(p\bmod a)}\)（即\((p\bmod a)^{-1}\)），又因为要保证等式右边是正数，所以我们将其加上\(p\)再\(\bmod p\)，这样就能将其变为正数。
综上，递推方程式如下：
初始化如下：
\(tan90^\circ\)即为不存在，但是我们赋值赋\(0\)。
#include<iostream>
#include<cstdio>
#define maxn 3000005
#define ll long long
using namespace std;
int a,b;
ll inv[maxn];
int main(){
scanf("%d%d",&a,&b);
inv[1]=1; printf("%lld\n",inv[1]);
for(ll i=2;i<=a;i++) inv[i]=(ll)b-(ll)b/i*inv[b%i]%b,printf("%lld\n",inv[i]);
return 0;
}
例：【模板】乘法逆元 2
大意
给定 \(n,p,k\) 和长为 \(n\) 的数组 \(a\)，求出 \(\sum\limits_{i=1}^n{\dfrac{k^i}{a_i}}\)。
思路
小学都会的通分。
记 \(t = \prod\limits_{i=1}^n{a_i}\)，\(pre_i\) 和 \(las_i\) 分别为 \(a\) 的前缀积与后缀积。有：
只需要求 \(t\) 的逆元即可。
#include<iostream>
#include<cstdio>
#define maxn 5000005
#define ll long long
using namespace std;
inline ll read(){
char ch; ll res=0; bool neg=1; ch=getchar(); while(ch<'0'||ch>'9'){if(ch=='-') neg=-1;ch=getchar();}
while(ch>='0'&&ch<='9'){res=(res<<3)+(res<<1)+ch-'0';ch=getchar();} return res*neg;
}
ll n,p,k,mul=1,ans=0; ll a[maxn],pre[maxn],las[maxn];
ll inv(int xx){if(xx==1) return 1; return (1LL*p-(p/xx*inv(p%xx)%p)%p)%p;}
int main(){
n=read(); p=read(); k=read();
pre[0]=1LL; for(int i=1;i<=n;i++){a[i]=read();pre[i]=(pre[i-1]*a[i])%p;}
las[n+1]=1LL; for(int i=n;i>=1;i--) las[i]=(las[i+1]*a[i])%p;
for(int i=1,mul=k;i<=n;i++,mul=mul*k%p){ans+=(mul*(pre[i-1]*las[i+1]%p)%p)%p; ans%=p;}
printf("%lld",ans*inv(pre[n])%p);
return 0;
}
完结撒花
