数据结构与算法 - 逆元

发表于 2021-01-24 更新于 2021-11-07 分类于数据结构与算法

对于加减乘运算都会满足：

　　\[ (a + b) \bmod p = a \bmod p + b \bmod p \]

　　\[ (a - b) \bmod p = a \bmod p - b \bmod p \]

　　\[ (a \times b) \bmod p = a \bmod p \times (b \bmod p) \]

但是对于除法而言 $ (a/b) p (a p) / (b p) $

这个可以通过举例验证一下或者可以通过 $b \bmod p$ 可以为 0 但是分母不能为 0 得到对于除法而言不满足这个运算的结论。

但是如果非要求 $(a/b) \bmod p$ 要怎么办呢？

逆元的定义

为了解决模意义下的除法问题，我们引入了逆元 $\text{inv}(a)$。 $\text{inv}(a)$ 其实可以看做模 $p$ 意义下的 $\displaystyle\frac{1}{a}$ ，那么在模 $p$ 意义下， $\displaystyle\frac{a}{b}$ 就可以变形为 $a\cdot \text{inv}(b))$

定义：如果一个线性同余方程 $ax \equiv 1 \pmod b$ ，则 $x$ 称为 $a \bmod b$ 的逆元，记作 $\text{inv}(a)$ 。

逆元的求解

这里介绍三种计算逆元的方法：拓展欧几里得，费马小定理，线性递推。

拓展欧几里得

欧几里得算法

如果我们已知两个数 $a$ 和 $b$ ，如何求出二者的最大公约数呢？

不妨设 $a > b$

我们发现如果 $b$ 是 $a$ 的约数，那么 $b$ 就是二者的最大公约数。下面讨论不能整除的情况，即 $a = b \times q + r$ ，其中 $r < b$ 。

我们通过证明可以得到 $\gcd(a,b)=\gcd(b,a \bmod b)$ ，过程如下：

设 $a=bk+c$ ，显然有 $c=a \bmod b$ 。设 $d|a\ \ \ d|b$ ，则 $c=a-bk$ $\frac{c}{d}=\frac{a}{d}-\frac{b}{d}k$ 由右边的式子可知 $\frac{c}{d}$ 为整数，即 $d|c$ 所以对于 $a,b$ 的公约数，它也会是 $a \bmod b$ 的公约数。

反过来也需要证明

设 $d|b\ \ \ d|(a \bmod b)$ ，我们还是可以像之前一样得到以下式子 $\frac{a\bmod b}{d}=\frac{a}{d}-\frac{b}{d}k$ $\frac{a\bmod b}{d}+\frac{b}{d}k=\frac{a}{d}$ 因为左边式子显然为整数，所以 $\frac{a}{d}$ 也为整数，即 $d|a$ ，所以 $b,a\bmod b$ 的公约数也是 $a,b$ 的公约数。

既然两式公约数都是相同的，那么最大公约数也会相同。

所以得到式子 $\gcd(a,b)=\gcd(b,a\bmod b)$

既然得到了 $\gcd(a, b) = \gcd(b, r)$ ，这里两个数的大小是不会增大的，那么我们也就得到了关于两个数的最大公约数的一个递归求法。

int gcd(int a, int b) {
    if (b == 0) return a;
    return gcd(b, a % b);
}

递归至 b==0 （即上一步的 a%b==0 ) 的情况再返回值即可。

拓展欧几里得算法

扩展欧几里德定理（Extended Euclidean algorithm, EXGCD），常用于求 $ax+by=\gcd(a,b)$ 的一组可行解。

证明：

设

$ax_1+by_1=\gcd(a,b)$

$bx_2+(a\bmod b)y_2=\gcd(b,a\bmod b)$

由欧几里得定理可知： $\gcd(a,b)=\gcd(b,a\bmod b)$

所以 $ax_1+by_1=bx_2+(a\bmod b)y_2$

又因为 $a\bmod b=a-(\lfloor\frac{a}{b}\rfloor\times b)$

所以 $ax_1+by_1=bx_2+(a-(\lfloor\frac{a}{b}\rfloor\times b))y_2$

$ax_1+by_1=ay_2+bx_2-\lfloor\frac{a}{b}\rfloor\times by_2=ay_2+b(x_2-\lfloor\frac{a}{b}\rfloor y_2)$

因为 $a=a,b=b$ ，所以可令 $x_1=y_2,y_1=x_2-\lfloor\frac{a}{b}\rfloor y_2$

将 $x_2,y_2$ 不断代入递归求解直至 GCD（最大公约数，下同）为 0 递归 x=1,y=0 回去求解。

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

函数返回的值为 GCD，在这个过程中计算 $x,y$ 即可。

线性同余方程

形如 $ax \equiv c \pmod b$ 的方程被称为线性同余方程 (Congruence Equation)。

根据以下两个定理，我们可以求出同余方程 $ax \equiv c \pmod b$ 的解。

定理 1 ：方程 $ax+by=c$ 与方程 $ax \equiv c \pmod b$ 是等价的，有整数解的充要条件为 $\gcd(a,b) \mid c$

根据定理 1，方程 $ax+by=c$，我们可以先用扩展欧几里得算法求出一组 $x_0,y_0$ ，也就是 $ax_0+by_0=\gcd(a,b)$，然后两边同时除以 $\gcd(a,b)$，再乘 $c$。然后就得到了方程 $a\dfrac{c}{\gcd(a,b)}x_0+b\dfrac{c}{\gcd(a,b)}y_0=c$，然后我们就找到了方程的一个解。

定理 2 ：若 $\gcd(a,b)=1$，且 $x_0$、$y_0$ 为方程 $ax+by=c$ 的一组解，则该方程的任意解可表示为：$x=x_0+bt$，$y=y_0-at$ , 且对任意整数 $t$ 都成立。

根据定理 2，可以求出方程的所有解。但在实际问题中，我们往往被要求求出一个最小整数解，也就是一个特解 $x=(x \bmod t+t) \bmod t$，其中 $t=\dfrac{b}{\gcd(a,b)}$。

代码如下：

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

扩展欧几里得法和求解线性同余方程是一个原理，在这里不展开解释。

费马小定理

费马小定理是数论里的重要定理，叙述如下：

若 $p$ 是质数，且 $\gcd(a,p)=1$ ，则有 $a^{p-1}\equiv1\pmod{p}$

从逆元的定义推导，可得：

\[ a\cdot\text{inv}(a)\equiv1\equiv a^{p-1}\pmod{p} \]

于是有：

\[\text{inv}(a)\equiv a^{p-2}\pmod{p}\]

于是对 $a^{p-2}$ 算一下快速幂就好了。注意这个方法只对 $p$ 是质数的情形有效。

代码如下：

// 快速幂算法 (a ^ n) mod p
inline int qpow(long long a, long long n, long long p) {
    long long ans = 1;
    a = (a % p + p) % p;
    for (; n; n>>= 1) {
        if (n & 1) ans = (a * ans) % p;
        a = (a * a) % p;
    }
    return ans;
}
// a mod b
inline long long inv(long long a, long long p) {
    return qpow(a, p - 2, p);
}

快速幂算法

快速幂算法解读（以 3^10 为例）：

3^10 = (3*3*3*3*3*3*3*3*3*3)

尽量想办法把指数变小来，这里的指数为 10

3^10 = (3*3)*(3*3)*(3*3)*(3*3)*(3*3)
     = (3*3)^5
     = 9^5

此时指数由 10 缩减一半变成了 5，而底数变成了原来的平方
求 3^10 原本需要执行 10 次循环操作，求 9^5 却只需要执行 5 次循环操作
但是 3^10 却等于 9^5, 我们用一次（底数做平方操作）的操作减少了原本一半的循环量，特别是在幂特别大的时候效果非常好
例如 2^10000=4^5000, 底数只是做了一个小小的平方操作，而指数就从 10000 变成了 5000，减少了 5000 次的循环操作。

现在我们的问题是如何把指数 5 变成原来的一半，5 是一个奇数，5 的一半是 2.5
但是我们知道，指数不能为小数，因此我们不能这么简单粗暴的直接执行 5/2
然而，这里还有另一种方法能表示 9^5

9^5 =（9^4）*（9^1）

此时我们抽出了一个底数的一次方，这里即为 9^1，这个 9^1 我们先单独移出来
剩下的 9^4 又能够在执行缩指数操作了，把指数缩小一半，底数执行平方操作

9^5 =（81^2）*(9^1)

把指数缩小一半，底数执行平方操作

9^5 =（6561^1）*(9^1)

此时，我们发现指数又变成了一个奇数 1
按照上面对指数为奇数的操作方法，应该抽出了一个底数的一次方，这里即为 6561^1
这个 6561^1 我们先单独移出来，但是此时指数却变成了 0，也就意味着我们无法再进行缩指数操作了

9^5 =（6561^0）*(9^1)*(6561^1)
    = 1*(9^1)*(6561^1)
    = (9^1)*(6561^1)
    = 9*6561
    = 59049

我们能够发现，最后的结果是 9*6561
而 9 是怎么产生的？是不是当指数为奇数 5 时，此时底数为 9
那 6561 又是怎么产生的呢？是不是当指数为奇数 1 时，此时的底数为 6561

所以我们能发现一个规律：最后求出的幂结果实际上就是在变化过程中所有当指数为奇数时底数的乘积