POJ 2115 C Looooops 题解 (数论)

C Looooops (From http://poj.org/problem?id=2115)

 

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 + Cx = B (mod 2^k) 最小非负整数 x, 即求解模线性方程 Cx = B - A (mod 2^k). 
　　　　　　令 a = C, b = B - A, n = 2^k, 则 方程变为 ax = b (mod n).
　　　　　　记 d = gcd(a, n), 若 d | b 则方程有解， 否则无解. 无解对应题目中的无限循环(FOREVER).
　　　　　　
　　　　　　由数论知识， 存在整数 x', y' 使得 d = ax' + ny'成立, 即 x' 为模线性方程 ax = d (mod n)的一个解.
　　　　　　则方程 ax = b (mod n) 的一个解为 x0 = x'* b/d, 其解的集合 {x0 + i * (n/d) | i = 0, 1, ..., d-1},
　　　　　　其中包含了方程对模 n 的 d 个不同的解.
　　　　　　
　　　　　　由于 ax = b (mod n), 则存在整数 y 使得 ax + ny = b. 有 a'x + n'y = b', 
　　　　　　其中 a' = a/d, n' = n/d, b' = b/d.
　　　　　　可知 方程 a'x = b' (mod n') 与方程 ax = b (mod n) 同解. 
　　　　　　由于 gcd(a', n') = 1, 所以对模 n'方程只有一个解.
　　　　　　则 (x0 mod n' + n') mod n'就是 方程 a'x = b (mod n') 和方程 ax = b (mod n) 的最小非负整数解.
　　　　　　
　　　　　　可以利用扩展欧几里得算法求出 x' 和 d, 进而可求出 x0 和最终答案.
　　　　　　
　　　　　　可参考《算法导论》第二版第三十一章有关数论的算法。

源代码   

#include <stdio.h>

typedef long long int int_t;

int_t EXTENDED_EUCLID(int_t a, int_t b, int_t *px, int_t *py) {
    if (b == 0) { *px = 1, *py = 0; return a; }
    else {
        int_t d, x, y;
        d = EXTENDED_EUCLID(b, a % b, &x, &y);
        *px = y, *py = x - (a/b)*y;
        return d;
    }
}

/**
 *    ax = b (mod n); 0=< a < n, |b| < n 
 */
int_t MODULAR_LINEAR_EQUATION_SOLVER(int_t a, int_t b, int_t n) { 
    int_t x, y, m;
    int_t d = EXTENDED_EUCLID(a, n, &x, &y);
    if (b % d != 0) {
        return -1;
    } else {
        x = x * b/d;
        m = n / d;
        return ((x+m)%m+m)%m; // why?
    }
}

int main() {
    int_t A, B, C, k;
    while (scanf("%lld%lld%lld%lld", &A, &B, &C, &k) != EOF && k != 0) {
        int_t x = MODULAR_LINEAR_EQUATION_SOLVER(C, B-A, 1LL << k);
        if (x == -1) {
            printf("FOREVER\n");
        } else {
            printf("%lld\n", x);
        }
    }
    return 0;
}