Poj 2115 C Looooops (模线性方程)

题意：对于循环语句for(i=A ; i!=B ;i +=C)，问在无符号的k位存储系统中循环几次才会结束。若在有限次内结束，则输出循环次数，否则输出死循环。

思路：例如k=16，也就是unsigned int，那么能保存2^16个数据，即最大数为65535，当循环使得i超过65535时，则i会返回0重新开始计数。

于是可得 循环次数 x=[(B-A+2^k)%2^k] /C

即 Cx=(B-A)(mod 2^k)  ，求最小正整数解。

以下参考：http://972169909-qq-com.iteye.com/blog/1104538

解方程：ax == b (mod n);【ax % n == b % n】 
设线性模方程的一个解为x0 
条件①：有d = gcd(a, n) 
条件②：有d = ax1 + ny, 由扩展欧几里得(Egcd)得到x1的值 
条件③：b % d == 0 (有解的条件) 
对条件③进行解释： 
原方程化为：ax + kn = b (设k为某一整数) 
那么如果a与n的最大公约数为d，那么ax + kn 必然可以提取一个d的因子，也就是说b必然有d这个因子，所以如果b%d!=0，说明b没有d这因子，与前面的结论相互矛盾，所以无解 

则x0 = x1*(b/d); 

证明： 
因为：容易求得d = gcd (a, n), 则存在一个x1、y使得d = ax1 + ny①
方程①2边同时模n得：d % n == ax1 % n② 
又因为：b % d == 0, 即b是d的倍数; 
所以(b/d)必为整数; 
所以由②得： b % n == b*d/d%n == d*(b/d) % n == ax1*(b/d) % n == ax % n 
所以很容易可以看出x = x1*(b/d)是方程的一个整数解，得证 

#include <cstdio>
#include <cstring>

#define i64 __int64

i64 Extended_Euclid (i64 a,i64 b,i64 &x,i64 &y)
{//扩展欧几里得算法,求ax+by=gcd(a,b)的一组解(x,y),d=gcd(a,b)
    i64 d;
    if (b==0)
    {
        x=1;y=0;
        return a;
    }
    d=Extended_Euclid(b,a%b,y,x);
    y-=a/b*x;
    return d;
}

bool MLE (i64 a,i64 b,i64 n,i64 &x0) 
//解模线性方程：ax == b (mod n);【ax % n == b % n】
{//无解返回false,有解返回true,解存放在x0
	i64 d,x,y;
	d=Extended_Euclid (a,n,x,y);
	if (b % d)
		return false;
	x0=x * (b/d);
	i64 t=n/d;
	if (t<0) t=-t;//以防万一，有的题目t有可能是负数
	x0=(x0%t + t)%t;
	return true;
}

int main ()
{
	i64 a,b,c,k,x0;
	while (scanf("%I64d%I64d%I64d%I64d",&a,&b,&c,&k),a||b||c||k)
	{
		if (MLE(c,b-a,(i64)1<<k,x0))
			printf("%I64d\n",x0);
		else
			printf("FOREVER\n");
	}
	return 0;
}