gcd算法以及exgcd

1.欧几里得算法

欧几里得是求最大公约数的经典算法，又称辗转相除法。
gcd函数就是用来求(a,b)的最大公约数的。
gcd函数的基本性质：
gcd(a,b)=gcd(b,a)=gcd(-a,b)=gcd(|a|,|b|)
gcd(a,b)=gcd(b,a mod b)

证明：a可以表示成a = kb + r，则r = a mod b
	 假设d是a,b的一个公约数，
	 则有 d | a , d | b ，而r = a - kb，
	 因此 d | r 因此d是(b,a mod b)的公约数
	 假设d 是(b,a mod b)的公约数，
	 则 d | b , d | r ，但是a = kb +r
	 因此d也是(a,b)的公约数
	 因此(a,b)和(b,a mod b)的公约数是一样的，其最大公约数也必然相等，得证
  

所以通过递归，最终可以使得其中一个参数为零，则另一个不为零的数则为最大公约数
板子如下：

typedef long long ll;
ll gcd(ll a,ll b){
	return b==0?a:gcd(b,a%b);
}

当然也可以求最小公倍数
最小公倍数lcm=a * b/gcd(a,b)

2.扩展算法

1.解ax+by=gcd(a,b)

对于不完全为 0 的非负整数 a，b，gcd(a,b)表示 a，b 的最大公约数，必然
存在整数对 x，y ，使得 gcd（a，b）=ax+by。

求解 x，y的方法的理解:
	设 a>b。
1，显然当 b=0，gcd（a，b）=a。此时 x=1，y=0；
2，a>b>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;
说明： 
	a-[a/b]*b即为mod运算。[a/b]代表取小于a/b的最大整数。
也就是ax1+ by1 == ay2+ b(x2- [a / b] *y2);
根据恒等定理得:
		x1=y2;
		y1=x2- [a / b] *y2;
这样我们就得到了求解 x1,y1 的方法：x1，y1 的值基于 x2，y2.
上面的思想是以递归定义的，因为 gcd 不断的递归求解一定会有个时候 b=0，所以递归可以结束。  

板子如下：

int gcd(int a,int b,int &x,int &y){
    if (b==0){
        x=1,y=0;
        return a;
    }
    int q=gcd(b,a%b,y,x);
    y-=a/b*x;
    return q;
}

可以这样思考:
对于a’=b,b’=a%b 而言，我们求得 x, y使得 a’x+b’y=Gcd(a’,b’)
由于b’=a % b=a - a / b * b
(这里的 / 是程序设计语言中的除法)
那么可以得到:
ax+by=gcd(a,b) =>
bx+(a - a / b * b)y = gcd(a’, b’) = gcd(a, b) =>
ay +b(x - a / b * y) = Gcd(a, b)
因此对于a和b而言，他们的相对应的p，q分别是 y和 (x-a/b*y)。
板子如下：

nt gcd(int a,int b,int &x,int &y){
    if (b==0){
        x=1,y=0;
        return a;
    }
    int q=gcd(b,a%b,x,y);
    int t=x;x=y;
    y=t-a/b*y;
    return q;
}

2.解ax+by=c

使用扩展欧几里德算法解决不定方程的办法
对于不定整数方程pa+qb=c，若 c mod Gcd(a, b)=0,则该方程存在整数解，
否则不存在整数解。
通过上面的程序即可求出一组x，y使得gcd（a，b）=ax+by
找其所有解的方法如下：
在找到x * a+y * b = Gcd(a, b)的一组解x0,y0后，
可以得到x * a+y * b = c的一组解:
x1 = x0 * (c/Gcd(a,b)),
y1 = y0 * (c/Gcd(a,b))，

3.解集

假设我们得出ax+by=c一组解x与y，如果x每减去一个t（t为常数），
那么y至少必须加上一个a’* t/b’
(其中a’=a/gcd(a,b)，b’=b/gcd(a,b))，
才能使等式两边成立,即
　　a(x-t)+b(y+a’ * t/b’)=ax+by=c
由于我们需要变换后的x和y仍然是整数，由于t已经是整数，即x-t也一定是整数，
现在只需考虑a’ * t/b’是否为整数，由于gcd(a’,b’)=1，
即a’一定不是b’的倍数，故必须令t=t * b’(即c必须是b’的倍数)，才使得a’ * t/b’为整数，即解集为：
　　x=x-t * b/gcd(a,b)，y=y+t * a/gcd(a,b)，t为任意整数
特别的，当gcd(a,b)=1，解集为：x=x-t * b，y=y+t * a

练习

1.扩展gcd-时间复杂性

题目内容：
计算循环语句的执行频次 for(i=A; i!=B ; i+=C) x+=1;
其中A,B,C,i都是k位无符号整数。

输入描述

A B C k, 其中0<k<32

输出描述

输出执行频次数，如果是无穷，则输出“forever”

输入样例

3 7 2 16

输出样例

2

k位无符号整数则数据范围为0~2^k(具体含义可以搜索无符号整数溢出)
所以设y，题目意思可以理解为
A+x * C - B=y * 2^k => x * C+y * 2^k=B-A
于是可以用exgcd求解
令a=C，b=2^k，c=B-A
首先判断B-A 是否是gcd(a,b)倍数，如果不是则无法求出一组x。
求出一组x0，y0后，则通解为
x=x0 * a/gcd(a,b)*t
y=y0 - b/gcd(a,b)*t
令n=a/gcd(a,b),m= b/gcd(a,b)
则可以得到
①x0 * a+y0 * b=c
②a(x0+n * t)+b(y0-m * t)=c
联立上式
可得
a * t * n - b * t * m = 0
a * n - b * m=0
也就是a * n=b * m，要想n,m最小就是使得a * n，b * m都为a，b的最小公倍数：
a * n = (a * b)/gcd(a,b) => n = b/gcd(a,b);
同理：m = a /gcd(a,b);
则：
x0 = (x % n + n)%n
y0 = (y % m + m)%m

#include<iostream>
using namespace std;
//A + xC - B = y * 2^k
int exgcd(int a,int b,int &x,int &y){
	if(b==0){
		x=1;y=0;
		return a;
	}
		int d=exgcd(b,a%b,x,y);
		int t=x;
		x=y;
		y=t-a/b*y;
	return d;
}
int main(){
	int A,B,C,k;
	cin >> A >> B >> C >> k;
	//特判A=B=C=K=0 
	if(A==0&&A==B&&B==C&&C==k) return !(cout << "forever");
	int b=1<<k,a=C,c=B-A,x=10,y=10;
	int d=exgcd(a,b,x,y);
	if(c%d){
		return !(cout << "forever");
	} 
	x=(x*(c/d))%b;
	x=(x%(b/d)+(b/d))%(b/d);//防止出现负数解
	cout << x;
	return 0;
}