2016 September 13

欧几里德算法

求最大公约数

命题：设a % b = r，求证：gcd(a, b) = gcd(b, r)。

由a % b = r，可令a = bq + r。
设u是a和b的公约数，则有u|a且u|b，那么u|(a-bq)，即u|r，可见u也是b与r的公约数。
反过来，设v是b和r的公约数，则有v|b且v|r，那么v|(bq+r)，即v|a，可见v也是a与b的公约数。
综上，a与b的公约数集合跟b与r的公约数集合是相同的，自然最大公约数也相等。

利用上述结论可在对数时间内求出两数的最大公约数。

int gcd(int a, int b) {
    return b ? gcd(b, a % b) : a;
}

拓展：

解二元一次方程

二元一次方程AX + BY = C有整数解的充要条件是gcd(A,B)|C。不失一般性，给定A和B，如何求整数X与Y满足：AX + BY = gcd(A,B)？

以下是通过xgcd过程求解的伪代码：

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

以a=99, b=78为例，手工计算过程如下：

99 = 78 * 1 + 21
78 = 21 * 3 + 15
21 = 15 * 1 +  6
15 =  6 * 2 +  3
 6 =  3 * 2

=> 3 = 15 - 6 * 2 = 15 - (21 - 15) * 2 = 15 * 3 - 21 * 2
     = (78 - 21 * 3) * 3 - 21 * 2
     = 78 * 3 - 21 * 11
     = 78 * 3 - (99 - 78) * 11
     = 99 * (-11) + 78 * 14

有些时候，需要求的是某个变量的最小正整数解，该怎么办呢？不妨设二元一次方程为：AX + BY = C，求最小的正整数解X。

  1. 记d=gcd(A, B)，如果d|C，那么解存在，否则无解。
  2. 等式两边除以d，约去公因子，A'X + B'Y = C'，该方程与原方程同解，此时必有gcd(A', B') = 1。
  3. 根据扩展欧几里德求出方程A'X + B'Y = 1的一组解(x, y)，那么对应原方程的一组解为(X, Y) = (C'x, C'y)。
  4. 最小的正解为X = X mod B，如果此时X为负，用X = X + abs(B)做下调整。

求解模线性方程

考虑求方程：ax = b(mod n)，其中a>0, n>0，求所有满足条件的x的值。

由ax = b(mod n)，得ax = kn + b，也就是ax + kn = b，显然用扩展欧几里德算法求解。

相关定理：

以下是求解的伪代码：

function modular_linear(a, b, n)
    (d, x', y') <- xgcd(a, n)
    if d | b
        x0 <- x'(b/d) mod n
        for i <- 0 to d-1
            print (x0 + i(n/d)) mod n
    else
        print "no solution"

以方程14x = 30(mod 100)为例，a=14，b=30，n=100，调用xgcd后(d, x, y) = (2, -7, 1)，由于2|30，所以有2个解，x0 = -7 * 15 mod 100 = 95，x1 = (95 + 50) mod 100 = 45。

乘法逆元

考虑计算(a / b) mod n，如果计算a/b比较麻烦，由于式子中含有除法，而除法取余不满足分配率，若能将除法转化为乘法就好了，因为乘法取余满足：(x * y) mod z = ((x mod z) * (y mod z)) mod z。下面试图找一个数r，使得(a / b) mod n = (a * r) mod n恒成立，将问题转化。

在上一节求解模线性方程的结论中，当a与n互质时，gcd(a,n)=1，方程ax=b(mod n)有唯一解。特别地，如果b=1，方程变为ax=1(mod n)，满足该条件的x称为a的模n乘法逆元。

可以证明，我们需要寻找的数r，其实就是b的模n乘法逆元。

调用一次扩展欧几里德：(1, x, y) <- xgcd(b, n)，即bx+ny=1，由于bx mod n = (-yn + 1) mod n = 1，因此x是b的模n乘法逆元。由于一般要求乘法逆元非负，所以x为负数时要将它转成正数。

例题：求乘法逆元

def xgcd(a, b):
    if b == 0:
        return (a, 1, 0)
    (d, x, y) = xgcd(b, a % b)
    xx, yy = y, x - a / b * y
    return (d, xx, yy)
m, n = map(int, raw_input().split())
(d, x, y) = xgcd(m, n)
print (x - x / n * n + n) % n

解模线性方程组

设有模线性方程组：x = r[i] (mod m[i])，其中i=1,2,…,n，求满足条件的最小正整数x。

先考虑只有两个的情况：
x = r[1] (mod m[1])
x = r[2] (mod m[2])
=> x = r[1] + k[1]m[1] = r[2] + k[2]m[2]
=> k[1]m[1] - k[2]m[2] = r[2] - r[1]
设d=gcd(k[1], k[2])，如果r[2]-r[1]不能被d整除，则方程组无解。
否则方程两边同时除以d，记A=m[1]/d，B=m[2]/d，C=(r[2]-r[1])/d，则化为
Ak[1] - Bk[2] = C，且gcd(A, B) = 1。
调xgcd(A, B, u0, v0)求出方程Ak[1] - Bk[2] = 1的一组解，
那么(u, v) = (Cu0, Cv0)是方程Ak[1] - Bk[2] = C的一组解。
将u=k[1]代回x = r[1] + k[1]m[1]得到方程组的一个特解，记为x0。
那么方程组的通解可表示为：x = x0 + t * lcm(m[1], m[2])，其是t为任意整数。
该式子等价于：x = x0 (mod lcm(m[1], m[2]))，相当于把最初的两个方程合并了。
一直合并下去，到只剩下一个方程为止，那么此时的x0即为答案。

参考代码：

#include <iostream>
using namespace std;
typedef long long LL;
LL gcd(LL a, LL b) {
    return b ? gcd(b, a % b) : a;
}
void xgcd(LL a, LL b, LL &x, LL &y) {
    if (!b) {
        x = 1; y = 0;
    } else {
        xgcd(b, a % b, y, x);
        y -= a / b * x;
    }
}
int main() {
    LL n, m[1005], r[1005], M, R, d, u, v, A, B, C;
    cin >> n;
    for (int i = 0; i < n; i++)
        cin >> m[i] >> r[i];
    M = m[0]; R = r[0];
    for (int i = 1; i < n; i++) {
        d = gcd(M, m[i]);
        if ((r[i] - R) % d) {
            R = -1;
            break;
        }
        A = M / d; B = m[i] / d; C = (r[i] - R) / d;
        xgcd(A, B, u, v);
        u = u * C % B;
        R = R + u * M;
        M = M * B;
        R = R % M;
        if (R < 0) R += M;
    }
    cout << R << endl;
    return 0;
}
Table of Contents