poj 2115 C Looooops

// 題意: 解同餘方程 Cx≡B-A ( mod 2^k ) ,求出x的最小解

#include<iostream>
#include <math.h>
using namespace std;
long long x_0,y_0,q;    
void extend_eulid(long long a,long long b)            //x_0*a+y_0*b=gcd(a,b)=q
{
    if(b==0)
    {
        x_0=1;y_0=0;q=a;    
    }    
    else
    {
        extend_eulid(b,a%b);
        int temp=x_0;
        x_0=y_0;y_0=temp-a/b*y_0;
    }
}
// 解同餘方程 Cx≡B-A ( mod 2^k ) , 相當于ax≡n (mod b),即 a=C , n=B-A , b=2^k 
int main()
{
    long long A,B,C,k,x;
    while(cin>>A>>B>>C>>k&&k)
    {
        extend_eulid(C,pow(2.0,double(k)));    
        if( (B-A)%q != 0 )
            cout<<"FOREVER\n";
        else
        {
            x = x_0*(B-A)/q;
            long long ff=pow(2.0,double(k))/q;
            long long X=(x%ff+ff)%ff;
            cout<<X<<endl;
        }
    }
    return 0;
}

/*
形如  ax≡n (mod b),可化成 ax + by = n, x有整數解的充分必要條件是 n % （a，b）==0，(a,b)表示gcd(a,b)
方程 a*x+b*y=n；我們可以先用擴充歐幾裡德算法求出一組x_0，y_0。
也就是  a * x_0 + b * y_0 =（a，b）；
然後兩邊同時除以（a，b），再乘以n。
這樣就得到了方程 a * x_0 * n/（a，b）+ b * y_0 * n/（a，b）= n；
方程的一個解: x = x_0 * n/（a，b）, y = y_0 * n/（a，b）
其他的解都為  x_1,2...≡x+(b/(a,b))*i  ,   y_1,2...≡y-(a/(a,b))*i (0 <=i<=(a,b)-1)。 i=0 時即為上面的一組解
總共有gcd(a,b)組解

令ff=b/(a,b),則x的最小非負整數解為 x = x_0 * n/（a，b）,  X = ( x % ff + ff ) % ff ;
*/