一元线性同余方程（扩展欧几里得算法，C Looooops）

---

前言

贝祖定理：
若a,b是整数,且gcd(a,b)=d，那么对于任意的整数x,y,ax+by都一定是d的倍数，特别地，一定存在整数x,y，使ax+by=d成立。

但如何求出满足条件的x、y就是我们接下来的任务。

需要注意的是这个方程是丢番图方程。
丢番图方程（Diophantine Equation）：有一个或者几个变量的整系数方程，它们的求解仅仅在整数范围内进行。最后这个限制使得丢番图方程求解与实数范围方程求解有根本的不同。丢番图方程又名不定方程、整系数多项式方程，是变量仅容许是整数的多项式等式。

作为不定方程，变量个数多于方程个数，将会存在多个解，我们只需要找到其中的特解，通过特解来表示通解。

若我们找到了一组特解x，y可使得方程成立,
那么通解可以表示为：
$X=x+\frac{b}{d}t$
$Y=y-\frac{a}{d}t$
$\frac{b}{d}$与$\frac{a}{d}$为互质关系，从而通过整数 t 的不同取值可以得到所有可能的值。
两式后面的余项在方程中可以进行抵消，从而不影响等式的成立。

---

一、推导

1.反递推起点

在欧几里得算法中，我们通过辗转相除可以得到：
当前的a = gcd，b = 0；
此时由于b 等于0，我们可以解出在方程 $ax+by=d$中，
x = 1；
y的值随意，不妨令其为 0；

2.反递推

设当前层的a、b为A、B，根据辗转相除可知下一层的a、b为B、A%B。

若我们已经求出来在下一层中满足ax+by=d的x、y，并以X、Y来表示。
那么此时满足方程：$B\ast X+（A\%B）\ast Y=d$
若$A=K\ast B+R$，有$R=A\%B=A-\frac{A}{B}\ast B$
$d=B\ast X+（A-\frac{A}{B}\ast B）\ast Y$
可化简得：

$d=A\ast Y+B\ast （X-\frac{A}{B}\ast Y）$

当前层的$x=Y，y=X-\frac{A}{B}\ast Y$
那么我们就可以根据这个来反推出在最初的式子中的x、y；

---

二、代码实现

在原来欧几里得算法的基础上进行修改，在求最大公约数的同时求出满足条件的x、y。

int exgcd(int a, int b, int &x, int &y){//最终返回的是gcd(a,b)
	if(!b) {
		x=1;y=0;
		return a; //到达递归边界返回
	}
	int r=exgcd(b,a%b,x,y);
	int temp=y; //把x y变成上一层的
	y=x-a/b*y;
	x=temp;
	return r; //得到a b的最大公因数
}
void exgcd(int a, int b, int &d, int &x, int &y){
	if(!b) d=a, x=1, y=0;
	else{
		exgcd(b,a%b,d,y,x);
		y-=x*(a/b);
	}
}

---

三、例题·C Looooops

C Looooops

问题描述

给定A,B,C,K
以$2^K$为模，问能否使得循环
$for(i=A;i!=B;i+=C)$结束：
若无法循环出结果，则输出$FOREVER$；
若可以结束，则输出循环次数。

多组输入，当$A==B==C==K==0$时结束程序。

---

分析

若循环可以结束，设循环次数为$x$，则存在式子：
$A+xC\equiv B(mod{2}^K)\to xC\equiv B-A(mod{2}^K)$
可以用拓展欧几里得求解：
$xC+y{2}^K=gcd(C,2^k)$

若$(B-A)\%gcd(C,2^K)!=0$，说明方程无解，输出"FOREVER";
当$(B-A)\%gcd(C,2^K)==0$，说明方程有解，再根据这个解求出满足原方程的最小正整数解。

注意：给定k后，要把k换算为$2^k$，在默认使用$longlong$的前提下，如果是使用位运算，要注意写为$k=(ll)1<<k;$
如果使用函数$pow()$,则不必注意。

代码

#include<iostream>
#include<cmath>
using namespace std;
typedef long long ll;

ll exgcd(ll a, ll b, ll &x, ll &y){//最终返回的是gcd(a,b)
	if(!b) {
		x=1;y=0;
		return a; //到达递归边界返回
	}
	ll r=exgcd(b,a%b,x,y);
	ll temp=y; //把x y变成上一层的
	y=x-a/b*y;
	x=temp;
	return r; //得到a b的最大公因数
}

int main(){
	ll a,b,c,k;
	ll x,y,d,e;
	while(1){
		cin>>a>>b>>c>>k;
		if(a==b&&b==c&&c==k&&k==0)break;
		k=pow(2,k);
		e=b-a;
		d=exgcd(c,k,x,y);
		if(e%d)
			cout<<"FOREVER"<<endl;
		else {
			ll s=k/d;
			ll ans=(x*e/d%s+s)%s;
			cout<<ans<<endl;
		}		
	}
}