POJ 2115 扩展欧几里得

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;
}