扩展欧几里得算法

扩展欧几里得算法

裴蜀定理（Bézout's lemma）

定义

设 \(a,b\) 是不全为零的整数，对任意整数 \(x,y\)，满足 \(\gcd(a,b)\mid ax+by\)，且存在整数 \(x,y\), 使得 \(ax+by=\gcd(a,b)\).

证明

对于第一点

由于 \(\gcd(a,b)\mid a,\gcd(a,b)\mid b\)

所以 \(\gcd(a,b)\mid ax,\gcd(a,b)\mid by\)，其中 \(x,y\) 均为整数

因此 \(\gcd(a,b)\mid ax+by\)

对于第二点

在欧几里得辗转相除求 \(\gcd\)的最后一步，即 \(b = 0\) 时：

显然有一对整数 \(x = 1,y=0\)，使得 \(a\times 1 + 0\times 0 = \gcd(a,0)\)。

若 \(b > 0\)，则 \(\gcd(a,b) = \gcd(b,a\mod b)\)

假设存在一对整数 \(x,y\),满足 \(b\times x + (a\mod b)\times y = \gcd(b,a\mod b)\)

因为 \(bx + (a\mod b)y = bx + (a-b\lfloor a/b \rfloor)y = ay - b(x-\lfloor a/b \rfloor y)\)

所以令 \(x' = y, y' = x -\lfloor a/b \rfloor y\),就得到了\(ax' + by' = \gcd(a,b)\)

对上述过程用数学归纳法，得证。

扩展欧几里得算法

由上述证明可以推出如何求 \(x,y\)。

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

定义变量 \(d,x_0,y_0\)，调用 d = exgcd(b, a % b, x, y);。注意在上述代码中，\(x_0,y_0\) 是以引用方式传递的。上述程序求出方程 \(ax + by = \gcd(a,b)\) 的一组特解 \(x_0,y_0\)，并返回 \(a,b\) 的最大公约数 \(d\)。

对于更为一般的方程 \(ax + by = c\)，它有解当且仅当 \(d|c\)。我们可以先求出 \(ax + by = d\) 的一组特解 \(x_0,y_0\)，然后令 \(x_0,y_0\) 同时乘上 \(c/d\)，就得到了 \(ax + by = c\) 的一组特解 \((c/d)x_0,(c/d)y_0\)。

事实上，方程 \(ax + by = c\) 的通解可以表示为：

\[x = \frac c d x_0 + k\frac b d,y = \frac c d y_0 - k\frac a d (k\in \Z) \]

其中 \(k\) 取遍整数集合，\(d = \gcd(a,b)\)，\(x_0,y_0\) 是 \(ax + by = \gcd(a,b)\) 的一组特解。

乘法逆元

若整数 \(b,m\) 互质，并且 \(b\mid a\)，则存在一个整数 \(x\)，使得 \(a\div b \equiv a\times x\mod m\)。称 \(x\) 为 \(b\) 的模 \(m\) 乘法逆元，记为 \(b^{-1}\mod m\)

因为

\[a\div b\equiv a\times b^{-1}\equiv a\div b\times b\times b^{-1}\mod m \]

所以

\[b\times b^{-1} \equiv 1\mod m \]

如果 \(m\) 是质数（此时用符号 \(p\) 代替 \(m\)）并且 \(b < p\)，根据费马小定理，\(b^{p -1}\equiv 1\mod p\)，即 \(b\times b^{p - 2}\equiv 1\mod p\)。

因此，当模数 \(p\) 为质数时，\(b^{p - 2}\) 即为 \(b\) 的乘法逆元。

如果只是保证 \(b,m\) 互质，那么乘法逆元可通过求解同余方程 \(b\times x \equiv 1\mod m\) 得到。

这里先按照推出的结论，给出用快速幂求逆元

int qpow(long long a, int b) {
      int ans = 1;
      a = (a % p + p) % p;
      for (; b; b >>= 1) 
      {
            if (b & 1) ans = (a * ans) % p;
            a = (a * a) % p;
      }
      return ans;
}