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;

}