ACM数论——欧几里得与拓展欧几里得
欧几里得算法:
欧几里德算法又称辗转相除法,用于计算两个整数a,b的最大公约数。
基本算法:设a=qb+r,其中a,b,q,r都是整数,则gcd(a,b)=gcd(b,r),即gcd(a,b)=gcd(b,a%b)。
int gcd(int a,int b) { return b ? gcd(b,a%b) : a; }
扩展欧几里德算法:
基本算法:对于不完全为 0 的非负整数 a,b,gcd(a,b)表示 a,b 的最大公约数,必然存在整数对 x,y ,使得 gcd(a,b)=ax+by。
证明:设 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;
根据恒等定理得:x1=y2; y1=x2-(a/b)*y2;
这样我们就得到了求解 x1,y1 的方法:x1,y1 的值基于 x2,y2.
上面的思想是以递归定义的,因为 gcd 不断的递归求解一定会有个时候 b=0,所以递归可以结束。
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; }
扩展欧几里德算法的应用主要有以下三方面:
(1)求解不定方程;
(2)求解模线性方程(线性同余方程);
(3)求解模的逆元;
用扩展欧几里得算法解不定方程ax+by=c:
bool linear_equation(int a,int b,int c,int &x,int &y) { int d=exgcd(a,b,x,y); if(c%d) return false; int k=c/d; x*=k; y*=k; //求得的只是其中一组解 return true; }
求出解之间的间隔,那么就可以求出模的线性方程的解集:
bool modular_linear_equation(int a,int b,int n) { int x,y,x0,i; int d=exgcd(a,n,x,y); if(b%d) return false; x0=x*(b/d)%n; //特解 for(i=1;i<d;i++) printf("%d\n",(x0+i*(n/d))%n); return true; }
用扩展欧几里得求解逆元是一种常用的方法
你是否经常遇到过类似的问题 ,(A/B)%Mod 。此时,要先计算B%Mod的逆元p, 其实他是用逆元解决的典型题目。但是在使用逆元时候你需满足一下两个条件才能保证得到正确的结果。
- gcd(B,Mod) == 1,显然素数肯定是有逆元的。
- 这里B需要是A的因子
下面就给出扩展欧几里得的典型式子:a*x + b*y = 1 。求得x即为a%b的逆元; y即为b%a的逆元。
另一种方法是:p = b^(Mod-2) % Mod,因为b^(Mod-1)%Mod = 1(这里需要Mod为素数),因为这种方法不常用,因此这里不再详细介绍。
下面就给出求解逆元的模版(代码非原创)
扩展欧几里德求逆元模板:
#include<iostream> #define __int64 long long using namespace std; //举例 3x+4y=1 ax+by=1 //得到一组解x0=-1,y0=1 通解为x=-1+4k,y=1-3k inline __int64 extend_gcd(__int64 a,__int64 b,__int64 &x,__int64 &y)//ax+by=1返回a,b的gcd,同时求的一组满足题目的最小正整数解 { __int64 ans,t; if(b==0){x=1;y=0;return a;} ans=extend_gcd(b,a%b,x,y);t=x;x=y;y=t-(a/b)*y; return ans; } //(a/b)%mod=c 逆元为p,(p*b)%mod=1 //(a/b)*(p*b)%mod=c*1%mod=c // (p*b)%mod=1 等价于 p*b-(p*b)/mod*mod=1其中要求p,b已知 等价于 ax+by=1 //其中x=p(x就是逆元),y=p/mod,a=b,b=b*mod 那么调用extend_gcd(b,b*mod,x,y)即可求(a/b)%mod的逆元等价于a*p%mod int main() { __int64 a,b,x,y,c,gcd,mod,p;//ax+by=c while(cin>>a>>b>>c) { gcd=extend_gcd(a,b,x,y); cout<<x<<" "<<y<<endl; if(c%gcd){cout<<"无解!"<<endl;continue;} cout<<"x="<<x*c/gcd<<" y="<<y*c/gcd<<endl; } return 0; }
void extend_Euclid(int a, int b) { if(b==0) { x = 1; y = 0; return; } extend_Euclid(b, a%b); int t = x; x = y; y = t - a/b*y; } int main() { //b%mod的逆元 int b,mod; while(cin>>b>>mod){ // x=0;y=0; extend_Euclid(b,mod); cout<<(x%mod+mod)%mod<<endl; } return 0; }