C Looooops
Time Limit: 1000MS Memory Limit: 65536K
Total Submissions: 21071 Accepted: 5717

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 < 2k) modulo 2k.

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 < 2k) 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,B,C,k,求
for (variable = A; variable != B; variable += C)statement;
在二進位制k位數下，A每次增加C，增加到B一共進行多少次運算（二進位制每一位改變了就算一次運算）。
題解：
容易得到A + nC = B（mod 2^k）,
也就是 nC = B - A（mod2^k）.
就是求解模線性方程。
我們將模線性方程變形
x*C + y*2^k = B - A(mod2^k)
也就是求這個不定方程的x 最小正整數解
考慮對a*x + b*y = gcd（a,b);
如果我們用輾轉相除求gcd（a，b）是這樣的：
r1 = a%b
r2 = b%r1
r3 = r1%r2
r4 = r2%r3
……
rn = rn-2%rn-1
最終一定有rn = 0。那此時rn-1就是所求的gcd(a,b)。而擴充套件歐幾里得就是在原有的輾轉相除基礎上再倒回去。
可知：
r1 = a%b = a - k1*b = a - (a/b)*b
r2 = b%r1 = b - (b/r1)*r1
r3 = r1%r2 = r1 - (r1/r2)*r2
r4 = r2%r3 = r2 - (r2/r3)*r3
……
rn = rn-2%rn-1 = rn-2 - (rn-2/rn-1)*rn-1
例如：我們假設r3為最終步驟。有r3 = r1 - (r1/r2)*r2，帶入r1，r2
r3 = a - (a/b)b - (（a - (a/b)*b）/（b - (b/r1)*r1）)（b - (b/r1)*r1）
再帶入r1
r3 = a - (a/b)b - (（a - (a/b)*b）/（b - (b/r1)（a - (a/b)b））)（b - (b/（a - (a/b)b）)（a - (a/b)*b））
這就是原來的式子了。用一個例項來更好直觀的表示：
47x+30y=gcd（47，30） = 1

r1 = 17 = 47-30*1
r2 = 13 = 30-17*1
r3 = 4 = 17-13*1
r4 = 1 = 13-4*3
r5 = 0 = 4 - 1*4
倒回去：
r4 = 1=13*1-4*(3)
r4 = 1=13*1-[17-13*1]*3
r4 = 1=13*4-17*3
r4 = 1=[30-17*1]*4-17*3
r4 = 1=30*4-17*7
r4 = 1=30*4-[47-30*1]*7
r4 = 1=30*11+47*(-7)
可知x = -7，y = 11。而gcd（47，30） = 1；
因此擴充套件歐幾里得的演算法寫為

void Exgcd(ll a,ll b,ll& d,ll& x,ll& y){
    if(!b){
        d = a;x = 1;y = 0;
    }else{
        Exgcd(b,a%b,d,y,x);
        y -= x*(a/b);
    }
}

當迴圈到b = 0時，gcd就等於a，並且有1*a + 0*b = a。所以此時x = 1，y = 0
接著每一層和上一層的x和y都要交換，且y -= x*(a/b)。

如此我們求出了一a*x + b*y = gcd（a,b）的特解。但是本題求得是最小正整數解。

我們知道假設一個特解是x0，y0。則通解是x0+k*b/gcd(a,b),y0-k*a/gcd(a,b)。（假設x0*a + y0*b = x1*a + y1*b，（x0-x1）a = （y1-y0）*b，兩邊同除gcd（a，b），（x0-x1）a /d= （y1-y0）*b/d。由於a/d與b/d互質，所以x0-x1一定是b/d整數倍，由此得x1 = x0+k*b/gcd(a,b)）

由通解可知，每一個解的間距為b/gcd(a,b)，則最小整數解一定在（0，b/gcd(a,b)）內，則我們讓x = （x0%b/gcd(a,b)+b/gcd(a,b)）%b/gcd(a,b)，就可以讓它落到（0，b/gcd(a,b)）內。

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <map>
#include <queue>
#include <vector>
#define f(i,a,b) for(i = a;i<=b;i++)
#define fi(i,a,b) for(i = a;i>=b;i--)

using namespace std;
typedef long long ll;

void Exgcd(ll a,ll b,ll& d,ll& x,ll& y){
    if(!b){
        d = a;
        x = 1;
        y = 0;
    }else{
        Exgcd(b,a%b,d,y,x);
        y -= x*(a/b);
    }
}

int main()
{
    ll a,b,c,k;
    while(~scanf("%lld%lld%lld%lld",&a,&b,&c,&k)&&(a||b||c||k)){
        ll M = 1LL<<k;
        ll d,x,y;
        Exgcd(c,M,d,x,y);
        ll C = b-a;
        if(C%d!=0){
            printf("FOREVER\n");
        }else{
            x = (x*(C/d))%M; //乘上倍數
            x = (x%(M/d)+M/d)%(M/d);//最小正整數解。
            printf("%lld\n",x);
        }
    }
    return 0;
}