POJ 2115 线性模方程 （归根到底还是扩展欧几里得） 转载请标明出处http://www.cnblogs.com/MasNoon/p/5733264.html
题意：给你4个数A,B,C,K,求出（i=A;i!=B;i+=c）这个循环在int k 位里能执行几次。int K位其实就是在2^k这个范围内执行几次，超过2^k，又会重新从0计数。
根据题意列方程：A+X*C=B->x*C=(B-A)%2^K就是一个线性模方程。
那怎么解线性模方程呢？先给出几个定理
定理1：方程ax=b(mod n)对于未知量x有解，当且仅当gcd(a, n)|b（就是gcd(a,b)能整除b,即b%gcd(a,b)=0）
定理2：方程ax=b(mod n)或者对模n有d个不同的解其中d=gcd(a, n),或者无解.
定理3：假设方程ax=b(mod n)有解（即有d|b，其中d=gcd(a, n)），x0是该方程的任意一个解，则该方程对模n恰有d个不同的解，分别为：xi=x0+i(n/d)(i = 1, 2, …, d-1)
再来看方程x*C=(B-A)%2^K，由定理1可知，如果(B-A）%gcd(C,2^K)不等于0，则方程无解，即输出"FOREVER".
如果(B-A）%gcd(C,2^K)等于0，那方程可以写成X*C+n*2^K=B-A,学过扩展欧几里得的孩子都应该觉得很眼熟吧。下面讲讲扩展欧几里得，学过的同学就可以跳着看了。
也是先给出定理
定理：对于不完全为0的非负整数a，b，gcd（a，b）表示a，b的最大公约数d，必然存在整数对x，y，使得gcd（a，b）=d=ax+by
对于gcd(a,b) = d，对(a, b)用欧几里德辗转相除会最终得到(d, 0)。此时对于把a =d, b = 0 代入a*x + b*y = d，显然x = 1，y可以为任意值。 我们可以用a = d, b = 0的情况逆推出来任何gcd(a, b) = d 满足a*x + b*y = d的解。如果x0，y0是b*x + (a%b)*y = d 的解，那么对于a*x + b*y = d的解呢？
b*x + (a%b)*y = d → b*x + (a - [a/b]*b)*y = d → a*y + b*(x - [a/b]*y) = d 所以a*x + b*y = d的解 x1 = y0，y1= x0 - [a/b]*y0;
现在知道（B-A）%gcd(C,2^K)等于0，令L=（B-A）/gcd(C,2^K)
x1*C+y1*2^k=gcd(C,2^k) ①
方程①*L等于方程X*C+n*2^K=B-A
即x=x1*L(其实从这里可以推出线性模方程的定理一，如果（B-A）%gcd(C,2^K)不等于0，L肯定是个小数，而我们知道X是整数，x1也是整数，L*x1肯定不会等于X.所以无解）
又从线性模方程的定理三可知,xi=x0+i(n/d)(i = 1, 2, …, d-1)(n/d就是这里的2^k/gcd(C,2^k))
我们要求的肯定要是最小解，（废话：第一次相等的时候循环就停止了），xi=x0+i(n/d)->x0=xi%(n/d)这样最小解就出来了。
所以答案是x1*L%(2^k/gcd(C,2^k)).
下面贴代码
#include<stdio.h>
#include<cstring>
#include<iostream>
using namespace std;
__int64 ex_Euc( __int64 a, __int64 b, __int64 &x, __int64 &y)//扩展欧几里得，递归
{
if(b==0)
{
x=1,y=0;
return a;
}
__int64 d=ex_Euc(b,a%b,x,y);
__int64 t=x;
x=y;
y=t-a/b*y;
return d;
}
__int64 quickpow( __int64 a, __int64 b)
{
__int64 res=1;
while(b!=0)
{
if(b&1)
{
res*=a;
}
a=a*a;
b>>=1;
}
return res;
}
int main(void)
{
__int64 a,b,c,k;
while ( scanf("%I64d%I64d%I64d%I64d",&a,&b,&c,&k)!=EOF)
{
if(a==0&&b==0&&c==0&&k==0)
{
break;
}
__int64 m=quickpow(2,k);
__int64 x;
__int64 y;
__int64 g=ex_Euc(c,m,x,y);
__int64 n=b-a;
if(n%g!=0)
{
puts("FOREVER");本身就有空行不用再\n
continue;
}
x=(x*(n/g))%m;//(x*n%(m/g))
x=(x%(m/g)+m/g)%(m/g);//避免有负数，mod 运算，(a+b)%b=a%b
printf("%I64d\n",x);
}
return 0;
}
