POJ 2115 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 < 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

题目大意

对于for(i=A;i!=B;i+=C) 的循环语句，问在 k 位存储系统中循环几次才会结束。比如当 k=3 时，存储的数 num 在 0-7 之间循环，当 i 将大于 7 时，由 0 重新开始。若在有限次内结束，输出次数；否则输出 FOREVER。

解题思路

1、转化题意：设循环次数为 x ，则可转化为求解方程组 (A+C*x)%(2^k)=B，求最小整数解 x；

2、转化线性同余方程：(A+C*x)%(2^k)=B —>(C*x)%(2^k)=B-A

即 C*x≡B-A(mod 2^k)；

3、转化为一元二次方程 C*x+2^k*y=B-A；

4、使用扩展欧几里得算法，求解方程 C*x+2^k*y=gcd(C,2^k)，可求得该方程的解为 x0；

5、因为 ax+by=c 可由方程 ax+by=gcd(a,b) 两边同乘以 (c/gcd(a,b)) 得到，则方程 ax+by=c 的解为 x0*(c/gcd(a,b));

6、（有待理解）利用周期性变化求最小的非负整数解公式：

　　x1=(x1%(b/(gcd(a,b)))+(b/gcd(a,b)))%(b/gcd(a,b))；

　　即若方程 C*x+2^k*y=B-A 的一组整数解为 (x1,y1)，则其任意整数解为 (x1+t*(b/gcd(a,b)),y1-t*(a/gcd(a,b))) ( t 为任意整数)，其中b/gcd(a,b)，a/gcd(a,b) 分别为 x 值，y值的增长周期。

想了两天的扩展欧几里得，回来补问题

有关通解 x=x1+t∗bgcd(a,b)，y=y1−t∗agcd(a,b)的问题

对于方程 a*x1+b*y1=gcd(a,b)

其等价于 (a∗x1+t∗a∗bgcd(a,b))+(b∗y1−t∗a∗bgcd(a,b))=gcd(a,b) （t代表同一个值，且为任意整数）

化简可得： a∗(x1+t∗bgcd(a,b))+b(y1−t∗agcd(a,b))=gcd(a,b)

即通解为 x=x1+t∗bgcd(a,b) y=y1−t∗agcd(a,b)

即bgcd(a,b) agcd(a,b)分别是 x 值，y值的增长周期。

PS：有没有发现我开始学会使用markdown了，傻傻的分式用了好久的 /。致谢友人A O(∩_∩)O

代码实现

#include <iostream>
#include<cstdio>
using namespace std;
#define ll __int64
ll d;
ll ex_gcd(ll a,ll b,ll &x,ll &y)
{
if(!b)
{
d=a,x=1,y=0;
}
else
{
ex_gcd(b,a%b,y,x);
y-=(a/b)*x;
}
return d;
}
int main()
{
ll a,b,c,k;
while(~scanf("%I64d %I64d %I64d %I64d",&a,&b,&c,&k))
{
ll x,y;
if(a==0&&b==0&&c==0&&k==0) break;
ll mod=(ll)1<<k;
ll gcd=ex_gcd(c,mod,x,y);
ll com=b-a;
if(com%gcd!=0)
printf("FOREVER\n");
else
{
x=(x*(com/gcd))%mod;
x=(x%(mod/gcd)+(mod/gcd))%(mod/gcd);
printf("%I64d\n",x);
}
}
return 0;
}