/***求a,b最大公约数***/
long long gcd(long long a, long long b) {
if(b == 0)
return a;
else
return gcd(b, a % b);
}
证法一
证法二
/***扩展的欧几里德算法a*x + b*y = Gcd(a,b)的一组整数解,结果存在x,y中***/
void extend_gcd(long long a, long long b, long long& x, long long &y) {
if(b == 0) {
x = 1;
y = 0;
return;
}
extend_gcd(b, a % b, x, y);
long long tmp = x;
x = y;
y = tmp - a / b * y;
}
a * (x + b) + b * (y - a) = a * x + b * y = (a, b)。
然而在应用上,往往并不是如此简单,很多时候会求解不定方程a * x + b * y = n。这个时候还是应用上面的算法:
- 求(a,b), 设c = (a,b),如果! c|n,则不存在整数解。因为将上式左右两边都除以c,可以知道,左边为整数,右边为非整数,故矛盾。
- 将左右两边同时除以c,设得到新的方程为a' * x + b' * y = n',应用上述算法求a' * x + b' * y = 1的解(由第一步知道(a',b') = 1)。设结果为x', y'。
- x = x' * n' , y = y' * n'是方程a * x + b * y = n。这个比较好理解,将a' * x + b' * y = 1两边同时扩大n'倍就行了。
- x = x' * n' + t * b, y = y' * n' - t * a(t为整数)是原方程a * x + b * y = n的所有解。
/*
扩展欧几里得定理
扩展欧几里得定理:对于两个不全为0的整数a、b,必存在一组解x,y,
使得ax+by==gcd(a,b);
拓展欧几里得实现
下面面的程序中,x和y用全局变量保存
int gcd(int a,int b)
{
int t,d;
if(b==0)
{
x=1;
y=0;
//此时b==0,也就是说gcd(a,0)==a。原式变为ax+by==a=gcd(a,b)--> x==1,y==0
return a; //返回a,b最大公约数的值
}
d=gcd(b,a%b); //欧几里得求最大公约数应用
t=x;
x=y;
y=t-(a/b)*y; //不明处2
return d; //返回a,b最大公约数的值
}
//不明处2 解释 ax+by==gcd(a,b)(1)
说明规则,x,y表示第一次递归时的值,x1,y1表示第二次递归时的值。
那么gcd(a,b)==gcd(b,a%b),同时都代入式1,
有ax+by==b*x1+(a%b)*y1。将右边变形一下
b*x1+(a%b)*y1==b*x1+(a-(a/b)*b)*y1==a*y1+b*(x1-(a/b)*y1),
最终得到ax+by==a*y1+b*(x1-(a/b)*y1)
也就是说:
上一深度的x等于下一深度的y1,
上一深度的y等于下一深度的x1-(a/b)*y1。
*需要注意,上面推导时用的除法都是整型除法
因此,得到了不定式ax+by==gcd(a,b)的一组解,x、y。
那么对于一般的不定式ax+by==c,它的解应该是什么呢。
很简单,x1=x*(c/gcd(a,b)),y1=y*(c/gcd(a,b))
再深入一点,就解出这么一组解其实一般来说是解决不了什么问题的,
我们现在要得到所有的解,那么这所有的解究竟是什么呢?
假设d=gcd(a,b). 那么x=x0+b/d*t; y=y0-a/d*t;其中t为任意常整数。
*/
#include<stdio.h>
#include<stdlib.h>
/*
求关于 x 的同余方程 ax ≡ 1 (mod b)的最小正整数解。
即求ax=mb+r 1=nb+r
变形ax+(n-m)b=1,此方程即拓展欧几里德的应用ax+by=gcd(a,b),(n-m相当于y)
事实上ax ≡1(mod b) 有解的必要条件是gcd(a,b)|1,即gcd(a,b)=1;
使用拓展欧几里得可知 ax+by=1(x,y是整数)
*/
int Extragcd(int a,int b,int *x,int *y)
{
int d,t;
if(b==0) //递归调用终止条件,当根据欧几里得辗转相除法则,余数为0停止
{
*x=1;
*y=0;
return a;
}
else
{
d=Extragcd(b,a%b,x,y);
t=*x; //根据下一个x1,y1的值,倒推前一个x,y的值
*x=*y;
*y=t-a/b*(*y);
return d;
}
}
int main()
{
int a,b,x,y;
int ans;
freopen("mod.in","r",stdin);
freopen("mod.out","w",stdout);
scanf("%d%d",&a,&b);
ans=Extragcd(a,b,&x,&y);
if(ans!=1) return 0;
//根据若x是方程的一个解,则方程的所有解为x+k*b k为整数
x=x%b; //保证最小的正整数解x ,且x属于{0,1,2,3...b-1}
while(x<=0)
x+=b ;
printf("%d\n",x);
return 0;
}
ax+by=c问题
问题:ax+by=c,已知a、b、c,求解使该等式成立的一组x,y。其中a、b、c、x、y均为整数
a,b的最大公约数为gcd(a,b)。如果c不是gcd(a,b)的倍数,则该等式无解,因为等式左边除以gcd(a,b)是整数,
而等式右边除以gcd(a,b)后为小数。(根据解方程的时候,在等式的左右两边同时除以非0的整数,等式依然成立)
因此,只有当c是gcd(a,b)的倍数的时候,该等式有解。这样,可以通过计算使ax1+by1=gcd(a,b)成立的x1、y1,
然后有x=(c/gcd(a,b))*x1,y=(c/gcd(a,b))*y1,得到x,y。
问题就被转换为求使ax+by=gcd(a,b)成立的一组x,y。这可以用扩展欧几里德算法求解。如下:
根据欧几里得的辗转相除法,最终如果b为零,则gcd(a,b)=a,那么x=1,y=0为一组解。(一组重要的特解)
如果b不为零,根据欧几里德定理,可以设
ax1+by1=gcd(a,b)=gcd(b,a%b)=bx2+(a%b)y2=bx2+(a-(a/b)*b)y2
进一步化简可以得到 -> a*y2+b*(x2-(a/b)*y2)
化简后有x1=y2,y1=x2-(a/b)y2。因此x1,y1依赖于x2,y2,同理依次类推递归调用求出x3,y3,x4,y4……,
类似于辗转相除,直到b=0时,求出xn,yn,便可以由后往前倒推出x1,y1的值。
// 扩展欧几里德算法,解gcd(a, b) = ax + by
// 结果存储在x,y中,用户调用时保证a、b、c都是整数
// 返回a,b的最大公约数
int extended_euclid(int a, int b, int &x, int &y)
{
if(b == 0) // gcd(a, b) == a
{
x = 1;
y = 0;
return a;
}
int n = extended_euclid(b, a%b, x, y);
int tmp = x;
x = y;
y = tmp - static_cast<int>(a/b)*y;
return n;
}
// 等式ax+by=c,已知a、b、c,求x和y。
// 解该线性方程等同于解同余式ax = c(mod b)
// 返回值表示是否有解,true有解,false无解
bool linear_equation(int a, int b, int c, int &x, int &y)
{
int n = extended_euclid(a, b, x, y);
if(c%n)
return false;
int k = c/n;
x *= k;
y *= k;
return true;
}
linear_equation函数也可以用来解同余式ax=c(mod b)。
由ax=c(mod b),可以得到ax = mb+r;c = nb+r。化简可以得到ax+(n-m)b=c。调用linear_equation可以求出x。