SGU 106 The equation 【扩展欧几里得算法】

欧几里得算法：（原文链接）

1. 对于等式ax+by=c，a，b，c皆为整数，c如果是gcd(a, b)的倍数，则方程有解，否则方程无解。（定理1）
2. 因为等式ax+by=gcd(a, b)必定有解（定理1），所以可以解出来，解法如下：
   　　因为gcd(a, b) = gcd(b, a % b)，所以有bx1+(a%b)y1=gcd(a,b)，注意！此时x1并不等于x，y1也不等于y！这个过程可以循环（就想求最大公约数一样），将b视为新的a'，a%b视为新的b'，可以得到b'x2+(a'%b')y2=gcd(a,b)一直到b''''(n个')=0时，等式就变成：a''..'xn+0y=gcd(a,b)，此时根据辗转相除法gcd(a, b)就应该等于a''...'，所以xn = 1, yn = 任意数，取0就好。然后根据ax+by=bx1+(a%b)y1 => ax+by=bx1+(a-a/b*b)y1 => ax+by=ay1+b(x1+a/b*y1) 根据恒等定理得x=y1，y=(x1+a/b*y1)，接着……递归吧。
   
3. 对于等式ax+by=c，abc皆为整数且c是gcd(a, b)的倍数，且(x1, y1)是方程ax+by=gcd(a, b)一组整数解，则(x1*(c/gcd(a, b)), y1*(c/gcd(a, b)))是方程ax+by=c的一组解
4. 对于等式ax+by=c，abc皆为整数且c是gcd(a, b)的倍数，且(x1, y1)是方程一组整数解，那么方程的所有整数解为：x=x1+k*(b/gcd(a,b))，y=y1-k*(a/gcd(a,b))，k是任意整数，对于一组x, y，k相同。
   　　因为(x1, y1)是原方程的一组解，所以有ax1+by1=c，则易得a(x1+k*b)+b(y1-k*a)=c，但是x每增加b，y就会减少a，是否有更小的范围，设d是a，b的约数，但不是最大公约数（若最大公约数是1则d=1）易得a(x1+k*(b/d))+b(y1-k*(a/d))=c，可以知道增加量变小了，现在要求什么时候最小，当然是分母最大的时候变化值最小，也就是当x1+b/gcd(a,b)，就是x1的下一个整数解！

代码实现：返回值为gcd(a,b)

int ex_gcd(int a, int b, int &x, int &y) {
    if (b == 0) { x = 1, y = 0; return a; }
    int r = ex_gcd(b, a%b, x, y);
    int t = x;
    x = y; y = t - a/b *y;
    return r;
}

SGU 106 

裸的扩展欧几里得

/*
    扩展欧几里得算法
    计算给定区间内的不定方程ax+by+c=0的所有解的个数
    注意点：负数整除的进位问题，因此wa到爆了……
*/
#include <cstdio>
#include <algorithm>
#include <iostream>
using namespace std;
typedef long long LL;
LL ex_gcd(LL a, LL b, LL &x, LL &y) {
    if (b == 0) { x = 1, y = 0; return a; }
    LL r = ex_gcd(b, a%b, x, y);
    LL t = x;
    x = y; y = t - a/b *y;
    return r;
}
int main() {
    LL a, b, c, x1, x2, y1, y2;
    cin >> a>>b>>c>>x1>>x2>>y1>>y2;
    c = -c;
    if (a == 0 && b == 0) {
        if (c == 0) cout << (x2-x1+1)*(y2-y1+1) << endl;
        else cout << 0 << endl;
        return 0;
    }
    if (a == 0) {
        if (c % b == 0 && y1 <= c/b && c/b <= y2) cout << x2-x1+1 << endl;
        else cout << 0 << endl;
        return 0;
    }
    if (b == 0) {
        if (c % a == 0 && x1 <= c/a && c/a <= x2) cout << y2-y1+1 << endl;
        else cout << 0 << endl;
        return 0;
    }

    LL x, y, gcd;
    gcd = ex_gcd(a, b, x, y);
    if (c % gcd != 0) {
        cout << 0 << endl;
        return 0;
    }

    c /= gcd;
    x *= c; x1 -= x; x2 -= x;
    y *= c; y1 -= y; y2 -= y;
    c = -y2; y2 = -y1; y1 = c;
    b /= gcd;
    a /= gcd;

    if (b < 0) { c = -x2; x2 = -x1; x1 = c; b = -b; }
    if (a < 0) { c = -y2; y2 = -y1; y1 = c; a = -a; }
    LL t1, t2, r1, r2;
    t1 = (x1<0 || x1%b == 0) ? x1/b : x1/b+1;
    t2 = (x2>0 || x2%b == 0) ? x2/b : x2/b-1;
    r1 = (y1<0 || y1%a == 0) ? y1/a : y1/a+1;
    r2 = (y2>0 || y2%a == 0) ? y2/a : y2/a-1;

    cout << min(t2, r2)-max(t1, r1) + 1 << endl;

    return 0;
}