C Looooops

看了半天的同余 扩展欧几里得 练练手

Description

A Compiler Mystery: We are given a C-language style for loop of type

for (variable = A; variable != B; variable += C)
  statement;

I.e., a loop which starts by setting variable to value A and while variable is not equal to B, repeats statement followed by increasing the variable by C. We want to know how many times does the statement get executed for particular values of A, B and C, assuming that all arithmetics is calculated in a k-bit unsigned integer type (with values 0 <= x < 2^k) modulo 2^k.

Input

The input consists of several instances. Each instance is described by a single line with four integers A, B, C, k separated by a single space. The integer k (1 <= k <= 32) is the number of bits of the control variable of the loop and A, B, C (0 <= A, B, C < 2^k) are the parameters of the loop.

The input is finished by a line containing four zeros.

Output

The output consists of several lines corresponding to the instances on the input. The i-th line contains either the number of executions of the statement in the i-th instance (a single integer number) or the word FOREVER if the loop does not terminate.

Sample Input

3 3 2 16
3 7 2 16
7 3 2 16
3 4 2 16
0 0 0 0

Sample Output

0
2
32766
FOREVER

题目的意思就是 a+bx≡c(mod 2^k)典型的线性同余方程 转化为 bx-y*2^k=c-a  当gcd(b,2^k)能整除c-a的时候 就存在解通解好求 跑一边exgcd可以得出x0 然后 基础解 x=(c-a)/gcd(b,2^k)*x0;  通解为 x0+zz*(2^k/gcd(2^k,b); 显然当zz=0的时候最小关键是要求出最小正整数解。对于ax+by=gcd(a,b) 这个方程 我们有通解 x=x0+b/gcd(a,b);那么怎么求x的最小正整数解。呢 首先知道一点 b/gcd(a,b)是解的最小区间这个怎么理解呢

假设c为x的解的最小间距，此时d为y的解的间距，所以x=x0+c*t，y=y0-d*t（x0，y0为一组特解，t为任意整数）

带入方程得：a*x0+a*c*t+b*y0-b*d*t=n,因为a*x0+b*y0=n，所以a*c*t-b*d*t=0，t不等于0时，a*c=b*d

因为a,b,c,d都为正整数，所以用最小的c，d，使得等式成立，ac，bd就应该等于a，b的最小公倍数a*b/gcd(a,b),

所以c=b/gcd(a,b)，d就等于a/gcd(a,b)。

所以，若最后所求解要求x为最小整数，那么x=(x0%（b/gcd(a,b))+b/gcd(a,b))%(b/gcd(a,b))即为x的最小整数解。

x0%（b/gcd(a,b))使解落到区间-b/gcd(a,b)~b/gcd(a,b)，再加上b/gcd(a,b)使解在区间0~2*b/gcd(a,b)，

再模上b/gcd(a,b)，则得到最小整数解 （转自http://www.xuebuyuan.com/1380732.html）

学无止境 挺有意思的这个求最小正整数解的方式

#include <cstdio>
#include <iostream>
#include <string>

using namespace std;
typedef long long ll;
ll exgcd(ll a,ll b,ll &x,ll &y)
{
    if(b==0)
    {
        x=1;
        y=0;
        return a;
    }
    ll temp=exgcd(b,a%b,y,x);
    y-=(a/b)*x;
    return temp;
}
ll get(ll k)
{
    ll temp=1;
    for(int i=1;i<=k;i++)
    {
        temp*=2;
    }
    return temp;
}
int main()
{
    ll a,b,c,k;
    while(cin>>a>>b>>c>>k)
    {
        if(a==0&&b==0&&c==0&&k==0) break;
        ll temp=b-a;
        ll x,y;
        ll zz=get(k);
        //cout<<zz<<endl;
        ll g=exgcd(c,zz,x,y);
        //cout<<g<<endl;
        if(temp%g!=0) cout<<"FOREVER"<<endl;
        else
        {
            ll t=zz/g;//这里有精度问题。。 得多注意
            x=(x*(temp/g)%t+t)%t;// 得到目标的最小解
            cout<<x<<endl;
        }
    }
    return 0;
}