POJ 2115 C Looooops （模线性方程）

坑点：

long long int mod = 1LL << k，否则会炸int。

题意：

给你A B C K四个变量，问你对于这样一个for循环在循环多少次之后跳出
for(int i = A ; i != B; i += c){
nothing;
}
注意这里对 i 的bit做了限制，只有k个bit，默认为无符号，也就是说 ${ i \pmod {2^{k}}}$

思路：

设次数为x，得到方程：

$${A + C * x = B \pmod {2^{k}}}$$
稍微变个形： $${ C * x = B - A\pmod {2^{k}}}$$
于是乎形如 ${ a * x = b \pmod n}$ （1） 的方程就新鲜出炉啦，所以啊，我们只需要把在方程（1）的x的解得到，本题的解也就得到了。

关于模线性方程的通解：

形如这样的方程：

$${ a * x = b \pmod n}$$
我们可以把他转化为： $${a*x + n*(-y) = b \tag{1} }$$
首先判断是否有解，利用Bezout定理：若设a,b是整数，则存在整数x,y，使得 ${ax+by=gcd(a,b)}$
又因为 ${d = gcd(a,n)}$，所以当 ${ b \mod d = 0 }$ 时，即b是d的整数倍时，有解。

接着，根据 ${d = gcd(a,n)}$ 我们得到这个式子：

$${a*x' + n*y' = d \tag{2}}$$
是不是感觉跟（1）有些像啊，所以咱们从这里入手，把（2）的左右同乘 b/d得到下式：
$${a * (x'*b/d) + (n*b/d)*y' = b \tag{3}}$$
把（1）和（3）式一对比，哎呦我去，想要知道（1）式的x，我们只需知道（3）式的x’就够了。

dakala（所以）。。。
接着求x’就好了，，，

而怎么求x’呢？答案：扩展欧几里得算法。
根据式子（2）运用${extend-gcd(a,n,x,y)}$ 算出x’。所以 ${x =(x'*b/d) }$ 也就得到了。
扩展欧几里得算法：

lli x = 0,y = 0;
lli e_gcd(lli a,lli b,lli &x,lli &y){
    if(b == 0){
        x = 1;
        y = 0;
        return a;
    }
    lli ans = e_gcd(b,a%b,x,y);
    lli temp = x;
    x = y;
    y = temp - a/b*y;
    return ans;
}

但是这里又有个问题，解出来的x有可能是负的，而我们想要的是正的，所以还要进行一步操作：
${x = (x \mod (n/d) + n/d) \mod (n/d)}$
有的童鞋可能要问了，这里为什么周期是n/d呢？请看（2）式，在这个式子里，如果d不等于1的话，左右两边是可以约分的，也就是说，周期并没有原来的n那么大，所以是n/d。

AC代码：

#include <iostream>
#include <cstdio>
#include <cmath>

typedef long long int lli;
using namespace std;

lli x = 0,y = 0;
lli e_gcd(lli a,lli b,lli &x,lli &y){
    if(b == 0){
        x = 1;
        y = 0;
        return a;
    }
    lli ans = e_gcd(b,a%b,x,y);
    lli temp = x;
    x = y;
    y = temp - a/b*y;
    return ans;
}

// ax = b (mod n)
lli mod_line(lli a,lli b,lli n){
    lli d = e_gcd(a,n,x,y);
    if(b % d == 0){
        lli minx0 = (x * b / d) % n;
        return (minx0 % (n/d) + n/d)% (n/d);
    }
    else
        return -1;
}

int main()
{
    lli A,B,C,k;
    while(scanf("%I64d%I64d%I64d%I64d",&A,&B,&C,&k),k){
        lli status = mod_line(C,B-A,1LL<<k);
        if(status == -1){
            puts("FOREVER");
        }
        else{
            printf("%I64d\n",status);
        }
    }
    return 0;
}