一:求解不定方程 ax+by=c
先给出两个介绍扩展欧几里德很详细的网站:
http://www.cnblogs.com/frog112111/archive/2012/08/19/2646012.html
http://www.cnblogs.com/ka200812/archive/2011/09/02/2164404.html
扩展欧几里德的递归代码:
//扩展欧几里德的递归代码: int exgcd(int a,int b,int &x,int &y) { if(b==0) { //b==0,也就是说gcd(a,0)==a。 //原式变为ax+by==a --> x==1,y==0。 x=1; y=0; return a; } int d=exgcd(b,a%b,x,y); int t=x; x=y; //更新当前的x y=t-a/b*y; //更新当前的y return d; //返回最大公约数 }
扩展欧几里德的应用:
(一):求解不定方程:
利用扩展欧几里德,可以求出
ax+by=gcd(a,b)的一组解,x1、y1
则通解为:
x = x1+b/gcd*i
y = y1-a/gcd*i
那么对于一般的不定式ax+by=c,
若gcd(a,b)整除c,则方程有解
一组解为:x0=x1*(c/gcd(a,b)),y0=y1*(c/gcd(a,b))。
通解为:x=x0+b/gcd*i,y=y0-a/gcd*i。
否则无解。
若ax+by=c有解,且gcd(a,b)=d,x的特解为x*(c/d)。
令ans=x*(c/d),s=b/d;
则方程的x最小整数解为:(ans%s+s)%s;
(二):求解同余方程
ax+by=c,实际上它等价于同余方程,ax =c mod b,这样同余方程就可以用扩展欧几里德算法求解。
(三):求解模的逆元
另外如果a b互质那么ax+by=gcd(a,b) = 1,转换为同余方程ax = 1 mod b这样x就变成了a的关于mod b的逆元。
求解扩展欧几里德过后,只需x=(x%b+b)%b,求出最小正整数解x,即为a的逆元。
也就是说求逆元也不过是扩展欧几里德算法的一个特殊情况。
两个简单的求解不定方程的例子:
POJ 2115 C Looooops
#include <iostream> #include <stdio.h> using namespace std; long long A,B,C; int k; long long exgcd(long long a,long long b,long long &x,long long &y) { if(b==0) { x=1; y=0; return a; } long long d=exgcd(b,a%b,x,y); long long tmp=x; x=y; y=tmp-(a/b)*y; return d; } int main() { while(scanf("%I64d%I64d%I64d%d",&A,&B,&C,&k)!=EOF) { if(A==0 && B==0 && C==0 && k==0) break; long long a,b,c,x,y,d,s; a=C; b=(long long)1<<k; c=B-A; d=exgcd(a,b,x,y); if(c%d==0) { x=x*c/d; s=b/d; x=(x%s+s)%s; printf("%I64d ",x); } else printf("FOREVER "); } return 0; }
POJ 2142 The Balance
http://www.cnblogs.com/chenxiwenruo/p/3557315.html
二:解方程x^2+y^2=z^2
设不定方程:x^2+y^2=z^2
若正整数三元组(x,y,z)满足上述方程,则称为毕达哥拉斯三元组。
若gcd(x,y,z)=1,则称为本原的毕达哥拉斯三元组。
定理:
正整数x,y,z构成一个本原的毕达哥拉斯三元组且y为偶数,当且仅当存在互素的正整数m,n(m>n),其中m,n的奇偶性不同,
并且满足
x=m^2-n^2,y=2*m*n, z=m^2+n^2
例题:
POJ 1305 Fermat vs. Pythagoras
http://www.cnblogs.com/chenxiwenruo/p/3558297.html
三:佩尔方程x^2-dy^2=1
形如x^2-dy^2=1 (d>1且d不为完全平方数)的不定方程称为佩尔方程。
若d是完全平方数,那么x^2-dy^2=1是无解的。
求解佩尔方程,可以用暴力求解法。
还有一种更高效的解法,连分数法,有兴趣的可以自己看看。
暴力解法代码:
void solve(int d){ y=1; while(1){ x=(long long)sqrt(double(d*y*y+1)); if(x*x-d*y*y==1){ break; } y++; } }
给出两个简单例题:
POJ 1320 Street Numbers
求解两个不相等的正整数n和m (n<M),使得1+2+3+…+n=n+(n+1)+…+m。
思路:实际上求解(2m+1)^2-8n^2=1。令x=2m+1,y=n。可以利用迭代公式求解出前十组的(n,m)。
HDU 3292 No more tricks, Mr Nanguo
求解X^2-NY^2=1的第K大的满足式子的X,如果无解则输出相应的句子。利用矩阵迭代+快速幂即可求出。
#include <iostream> #include <cstdio> #include <string.h> #include <algorithm> #include <math.h> using namespace std; const int mod=8191; int n,k; int vis[30]; long long x,y; void solve(int d){ y=1; while(1){ x=(long long)sqrt(double(d*y*y+1)); if(x*x-d*y*y==1){ break; } y++; } } struct Matrix{ long long m[2][2]; void init(long long x,long long y){ m[0][0]=x;m[0][1]=n*y%mod; m[1][0]=y;m[1][1]=x; } void initE(){ memset(m,0,sizeof(m)); m[0][0]=m[1][1]=1; } }A; Matrix operator*(Matrix a,Matrix b){ Matrix c; for(int i=0;i<2;i++){ for(int j=0;j<2;j++){ c.m[i][j]=0; for(int k=0;k<2;k++){ c.m[i][j]=(c.m[i][j]+a.m[i][k]*b.m[k][j]%mod)%mod; } } } return c; } Matrix quickPow(Matrix a,int b){ Matrix ret; ret.initE(); while(b){ if(b&1) ret=ret*a; a=a*a; b=b>>1; } return ret; } int main() { memset(vis,0,sizeof(vis)); for(int i=1;i*i<30;i++){ vis[i*i]=1; //若是完全平方数,则无解,标记为1 } while(scanf("%d%d",&n,&k)!=EOF){ if(!vis[n]){ solve(n); x=x%mod;y=y%mod; A.init(x,y); Matrix B; B=quickPow(A,k-1); long long xk,yk; xk=(B.m[0][0]*x%mod+B.m[0][1]*y%mod)%mod; printf("%I64d ",xk); } else{ printf("No answers can meet such conditions "); } } return 0; }
四:解高次同余方程A^x=B(mod C)
Baby Step Giant Step算法:复杂度O( sqrt(C) )
请戳: