一.扩展欧几里德(exgcd)
设 a>b。
1,显然当 b=0,gcd(a,b)=a。此时 x=1,y=0;
2,ab!=0 时
设 ax1+by1=gcd(a,b);
bx2+(a mod b)y2=gcd(b,a mod b);
根据欧几里德有 gcd(a,b)=gcd(b,a mod b);
则:ax1+by1=bx2+(a mod b)y2;
即:ax1+by1=bx2+(a-(a/b)*b)y2=ay2+bx2-(a/b)*by2=ay2+b(x2-(a/b)*y2);
根据恒等定理得:x1=y2; y1=x2-(a/b)*y2;
这样我们就得到了求解 x1,y1 的方法:x1,y1 的值基于 x2,y2.
上面的思想是以递归定义的,因为 gcd 不断的递归求解一定会有个时候 b=0,所以递归可以结束。
对于一般的不定式 a*x+b*y==c;有整数解得充分必要条件是(c % gcd(a,b)==0)。
若已经求出 a*x+b*y==gcd(a,b)的解。则原不定式的解是x1=x*(c/gcd(a,b)),y1=y*(c/gcd(a,b))。
那么原不定式的所有解xi=x1+b/gcd(a,b)*t; yi=y1-a/gcd(a,b)*t;其中t为任意常整数。
当gcd(a,b)==1时,形象化理解就是当你的x要+b时,你的y要-a,所以,把x变成非负整数的方法是:x=(x%b+b)%b。
void exgcd(int a,int b,int &x,int &y) { if(b==0){x=1;y=0;} else { exgcd(b,a%b); int t=x;x=y;y=t-a/b*y; } }
二.乘法逆元(模数为p)
求a^(-1)(mod p意义下)。
以下仅讨论gcd(a,p)==1的情况,否则不存在逆元。
在 mod p的意义下我们把x的乘法逆元写作 x^(-1)。
乘法逆元有如下的性质:x*x^(-1)≡1(mod p)。
应用:x/y≡x*y^(-1)(mod p)
1.费马小定理。
欧拉定理:当gcd(a,p)==1时,a^(φ(p))≡1(mod p)
a^φ(p)≡1(mod p)
a*a^(φ(p)-1)≡1(mod p)
∴a^(φ(p)-1)==a^(-1)。
用快速幂即可出解。
通常用在p为质数的情况:当p为质数时,φ(p)=p-1。
a^(p-2)==a^(-1)
2.exgcd。
求ax≡1(mod p)的x。
原方程可化为:ax+py=1。
由于gcd(a,p)==1,所以用exgcd求出x即可,最后把x弄成非负整数。
通常用此计算逆元。
3.线性求逆元(inv[i]表示i的逆元)。
要求p是质数。我们要求inv[1~p-1]。
设当前已知inv[1~i-1]。记a=p/i,b=p%i。则p=ai+b。
ai+b≡0(mod p)
方程两边同乘i^(-1)*b^(-1)
(ai+b)*i^(-1)*b^(-1)≡0(mod p)
a*b^(-1)+i^(-1)≡0(mod p)
i^(-1)≡-a*b^(-1)(mod p)
∴inv[i]=(-p/i*inv[p%i])%p=((p-p/i)*inv[p%i])%p。
初始化:inv[1]=1,inv[0]=1。
#include<cstdio> long long inv[10000000]; int main() { int n,p;scanf("%d%d",&n,&p);inv[0]=inv[1]=1; for(int i=2;i<=n;i++)inv[i]=((long long)(p-p/i)*inv[p%i])%(long long)p; for(int i=1;i<=n;i++)printf("%d ",inv[i]); return 0; }
通过这种方法,我们还可以用递归的方式用O(log p)的时间求某个数的逆元。
即我们要查inv[i],就dfs调用inv[p%i]。
三.BSGS
以下仅讨论p为质数且gcd(a,p)==1的情况
我们要求a^x≡b(mod p)(给出a,b)的最小非负整数的x
我们先考虑暴力枚举x。那么我们枚举的x应该何时停止?
由费马小定理得:a^(p−1)≡1(mod p)
a^x/a^(p−1)^m≡a^x(mod p)(m为非负整数)
a^(x-m*(p-1))≡a^x(mod p)
∵x-m*(p-1)==x%(p-1)
∴a^(x%(p-1))≡a^x(mod p)
所以,我们的x只要从0枚举到p即可。
设x=ki-j。
则a^(ki-j)≡b(mod p)
a^ki≡b*a^j(mod p)
同样的,我们的ki-j只要枚举到p即可。
那么我们令k=√p(上取整),把i从1枚举到k,j从0枚举到k,这样我们的ki-j就只会枚举到p了。
可是这样还是O(k^2)=O(p)的啊。
以下的b*a^j与a^ki均膜p。
我们先枚举j,把所有枚举出来的b*a^j压入hash表,如果有两个j,b*a^j的值是一样的,那么我们取大的一个j,因为这样ki-j最小。
然后再枚举i,对于每个a^ki,我们查表,如果查到一个j满足a^ki==b*a^j,那么我们停止搜索,此时的ki-j即为答案。
时间复杂度(可以使用hash表来代替map,降低时间复杂度,map好写)为O(√plogp)
typedef long long ll; map<ll,ll>mp; ll mod; ll ksm(ll a,ll p) { ll ans=1; while(p) { if(p&1)ans=(ans*a)%mod; a=(a*a)%mod;p>>=1; } return ans; } ll BSGS(ll a,ll b) { if(a%mod==0)return -1; mp.clear(); ll m=(ll)sqrt((double)mod)+1;ll am=ksm(a,m); ll now=1; for(int j=0;j<=m;j++){mp[(b*now)%mod]=j;now=(now*a)%mod;}now=1; for(int i=1;i<=m;i++){now=(now*am)%mod;if(mp[now])return i*m-mp[now];} return -1; }
四.中国剩余定理(CRT)
我们要求一个数x,满足以下这些同余方程(要求p1,p2,...,pn互质)。
x≡a1(mod p1)
x≡a2(mod p2)
...
x≡an(mod pn)
我们设x1是第一个同余方程的一个解,x2是第二个同余方程的一个解...xn是第n个同余方程的一个解
当x2~xn都是p1的倍数,x1,x3~xn都是x2的倍数......x1~xn-1都是xn的倍数时
我们的答案x=x1+x2+...+xn,此时的x可以符合每一个同余方程。
我们不妨假设此时的a1<p1,a2<p2...an<pn(否则可以mod p)
x1 mod p1==a1且是p2*p3*...*pn的公倍数。
x2 mod p2==a2且是p1*p3*...*pn的公倍数。
...
xn mod pn==an且是p1*p2*...*pn-1的公倍数。
我们知道:如果a mod b==1,那么ac mod b==c(c为正整数且c<b)
x1 mod p1==1且是p2*p3*...*pn的公倍数。x1*=a1。
x2 mod p2==1且是p1*p3*...*pn的公倍数。x2*=a2。
...
xn mod pn==1且是p1*p2*...*pn-1的公倍数。xn*=an。
记M=p1*p2*...*pn。
显然,当我们的xi是M/pi的倍数时,我们满足了第2个条件。
我们要找到一个ki,使得M/pi *ki≡1(mod pi),这样就能满足第1个条件。
显然,ki是M/pi的逆元(mod pi意义下)!
于是,xi=M/pi * (M/pi)^(-1) * ai。x=x1+x2+...+xn。
记Mi=M/pi。则x=Σai*Mi*Mi^(-1)。
这只是一个解,其余解可以表示为Σai*Mi*Mi^(-1) +kM(k为整数)(因为我们不管是加上还是减去M,所得的x仍会符合以上的同余方程)。
所以最小非负整数解x'=Σai*Mi*Mi^(-1) %M
扩展中国剩余定理(ExCRT)
考虑合并
x≡a1(mod p1)
x≡a2(mod p2)
x=a1+k1p1=a2+k2p2
k1p1-k2p2=a2-a1
考虑Xp1+Yp2=gcd(p1,p2)
k1=X(a2-a1)/gcd(p1,p2),k2=-Y(a2-a1)/gcd(p1,p2)
因为x=a1+k1p1 >=0 ,所以X为满足Xp1+Yp2=gcd(p1,p2)最小的>=0的数
这样a1+k1p1=a2+k2p2
所以,当x=a1+k1p1=a2+k2p2时可以满足上面两个方程,且x=a1+k1p1 + K*lcm(p1,p2)的x也都可以满足上面两个方程(K为任意自然数)
所以x≡a1+k1p1(mod lcm(p1,p2))的x,就能满足上面两个方程,且不满足该条件的x,不能满足上面两个方程