2013-team4/code/linear-mod-equation

从 Trac 迁移的文章

这是从旧校内 Wiki 迁移的文章，可能存在一些样式问题，您可以向 memset0 反馈。

原文章内容如下：

{{{
typedef long long LL;

const int N = 10009;
LL A[N], B[N], M[N];

LL gcd(LL a, LL b){
    while(b){
        LL c = a % b; a = b; b = c;
    }
    return a;
}
// 扩展 gcd 求 ax + by = gcd(a,b)
// 解的范围 |x| <= b, |y| <= a
// 若有一组解 x0, y0, 则 x0+k*dx, y0-k*dy 也是一组解, 与方程右边无关
// dx = b/gcd, dy = a/gcd
LL extgcd(LL a, LL b, LL &x, LL &y){
    if(b == 0){
        x = 1; y = 0; return a;
    }
    else{
        LL g = extgcd(b, a % b, y, x);
        y -= x * (a / b);
        return g;
    }
}
// 求 a 在模 m 下的逆元
// 不存在返回 -1
LL inv(LL a, LL m){
    LL x, y;
    if(extgcd(a, m, x, y) != 1) return -1;
    return (x + m) % m;
}
// 求解 A[i]x = B[i] (mod M[i])
// 方程编号 0 到 n-1
// 返回 x 的最小值 同时 x+m, x+2m, x+3m... 都是解
// 即 X % m = x
bool gao(int n, LL A[], LL B[], LL M[], LL &x, LL &m){
    x = 0; m = 1;
    for(int i = 0; i < n; i++){
        LL a = A[i] * m, b = B[i] - A[i] * x, d = gcd(a, M[i]);
        if(b % d) return false;
        LL t = b / d * inv(a / d, M[i] / d) % (M[i] / d);
        x += m * t;
        m *= M[i] / d;
    }
    return true;
}
}}}

typedef long long LL;
const int N = 10009;
LL A[N], B[N], M[N];
LL gcd(LL a, LL b){
    while(b){
        LL c = a % b; a = b; b = c;
    }
    return a;
}
// 扩展 gcd 求 ax + by = gcd(a,b)
// 解的范围 |x| <= b, |y| <= a
// 若有一组解 x0, y0, 则 x0+k*dx, y0-k*dy 也是一组解, 与方程右边无关
// dx = b/gcd, dy = a/gcd
LL extgcd(LL a, LL b, LL &x, LL &y){
    if(b == 0){
        x = 1; y = 0; return a;
    }
    else{
        LL g = extgcd(b, a % b, y, x);
        y -= x * (a / b);
        return g;
    }
}
// 求 a 在模 m 下的逆元
// 不存在返回 -1
LL inv(LL a, LL m){
    LL x, y;
    if(extgcd(a, m, x, y) != 1) return -1;
    return (x + m) % m;
}
// 求解 A[i]x = B[i] (mod M[i])
// 方程编号 0 到 n-1
// 返回 x 的最小值 同时 x+m, x+2m, x+3m... 都是解
// 即 X % m = x
bool gao(int n, LL A[], LL B[], LL M[], LL &x, LL &m){
    x = 0; m = 1;
    for(int i = 0; i < n; i++){
        LL a = A[i] * m, b = B[i] - A[i] * x, d = gcd(a, M[i]);
        if(b % d) return false;
        LL t = b / d * inv(a / d, M[i] / d) % (M[i] / d);
        x += m * t;
        m *= M[i] / d;
    }
    return true;
}