数论笔记：Exgcd 与乘法逆元的漫步

· · Personal

这一篇学习笔记真的是伤肝，写完之后要发的时候电脑死机了，结果没保存下来，哭了

我们接下来将要学习模意义下的乘法逆元，具体包含费马小定理（以及欧拉定理，快速幂）求逆元、扩展欧几里得求逆元，线性求连续逆元以及线性求不连续逆元。

为了让自己 AFO 之后还能看懂，写的尽可能细一点吧/kel

Part 1 参考资料与概述

掌握这篇学习笔记，你需要的知识是模意义下的运算以及基本的代数变换。

其他的数学基础与拓展是：

Part 2 Exgcd 扩展欧几里得

我们集中注意力在一个问题上：ax+by=\gcd(a,b)，这个不定方程如何解决？

首先存在 ax_1+by_1=\gcd(a,b)，其次当我们将方程用 (a\bmod b) 和 a 来代入的时候可以发现，同时还存在 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-\bigg(\left\lfloor\dfrac{a}{b}\right\rfloor\times b\bigg)。

于是我们代入原式并化简：

ax_1+by_1=bx_2+\bigg(a-\left\lfloor\dfrac{a}{b}\right\rfloor\times b\bigg)y_2

乘法分配律拆开括号：

ax_1+by_1=bx_2+ay_2-\left\lfloor\dfrac{a}{b}\right\rfloor\times by_2

提取公因式进行因式分解：

ax_1+by_1=ay_2+b(x_2-\left\lfloor\dfrac{a}{b}\right\rfloor y_2)

我们可以将两边简单地对等起来，得到一个具有普遍规律的式子：

x_1=y_2\quad y_1=x_2-\left\lfloor\dfrac{a}{b}\right\rfloor y_2

这个式子的意义是，我们可以缩减 x 和 y 的范围，递归求解了。

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),t=x;
    x=y;y=t-(a/b)*y;
    return d;
}

我们之所以要使用取引用，是因为同时返回两个值不是很方便，不如直接修改变量值。

Part 3 逆元基础

什么是逆元？

假设存在一个自洽的代数系统 S，存在一种运算 \circ。对于其中的某一元素 a，如果存在 a\circ e=a，那么我们称 e 是该系统中对于 \circ 运算的单位元。以具体例子说明，对于整数来说，乘法的单位元就是 1，因为 \forall x\in\mathbb{Z} 都一定存在 x\cdot 1=x。同样，对于整数加法来说单位元就是 0。

而如果存在一个元素 a 与另一元素 x，使得 a\circ x=e，那么 x 称为 a 在该系统中对于该运算的右逆元，记为 a^{-1}_{r}。反之，如果存在 x 满足 x\circ a=e，那么 x 称为 a 在该系统中对于该运算的左逆元，记为 a^{-1}_{l}。

而据我们所知，在整数系统中，存在相当多一部分的运算是左右对称的，即 a\circ x=x\circ a，那么我们简单地称 x 为 a 的逆元，记为 a^{-1}。对于我们目前研究的整数系统，存在最为常见的对称运算是加法与乘法。对于实数 x 来说，它的乘法逆元就是 \frac{1}{x}，这也是我们记作 x^{-1} 的原因之一。

我们之所以要在模意义下求乘法逆元，是因为在取模的条件下，式子 (a\div b)\bmod c 并不等于 (a\bmod c)\div(b\bmod c)，所以我们不得不使用乘法逆元来保证计算的正确性。

模意义下的乘法逆元

在模 p 意义下，如果 ax=1，那么 x 是 a 的乘法逆元，记为 a^{-1}。更加严谨地，存在：

ax\equiv 1\pmod p

为了使得题目更加容易求解，我们可以令 p 是质数，我们先来讨论这种情况。

Part 4 较高复杂度求解

费马小定理与快速幂

\bold{LinearExpectation}：费马小定理与欧拉定理的证明及应用

对于素数 p，当 \gcd(a,p)=1 时，存在：

a^{p-1}\equiv 1\pmod p

将原式替换进去，就可以得到：

a^{p-1}\equiv ax\pmod p

两边同时除以 a，可以得到：

x\equiv a^{p-2}\pmod p

显然，我们只需要求出 a^{p-2}\bmod p 即可。但是这样的一个式子，如果使用循环累乘，达到单个数 O(p) 的复杂度，就过于难看了，无法通过模板题。于是我们使用快速幂，以分治的思想进行解决，可以实现单个数 O(\log p)。对于单点求逆元，这已经是极其优秀的做法了。但是对于模板题所要求的 n 个连续整数逆元，这种做法的时间复杂度 O(n\log p) 则不是很容易被接受。

对于欧拉定理，我们可以得到更为一般的结论，即：

x\equiv a^{\varphi(p)-1}\pmod p

其中欧拉函数可以用链接博客中的方法线性预处理，得到结果后调用即可。

扩展欧几里得 Exgcd

对于题目中给出的式子：

ax\equiv 1\pmod p

很容易转化为一个不定方程：

ax+by=1

那么使用上文说明完毕的扩展欧几里得 Exgcd 很容易解决这个问题。事实上，我们只要在原方程两边同时除以 \gcd(a,b) 即可得到上述的式子。

但是对于每个数，我们的时间复杂度仍是 \log 级别的，虽然后面跟的不再是 p 了，但是这个复杂度，在我们看来仍然偏高了，我们考虑一种线性递推的做法。

Part 5 线性求解

线性求逆元

考虑在 O(n) 的线性时间复杂度中完成在模 p 意义下对于 1\sim n 逆元的求取。

显然，对于边界条件 a=1，存在 a^{-1}\equiv 1\pmod p，其原因是 1\times 1\equiv 1\pmod p 显然恒成立。接着我们考虑进行一般化，对于递归情况 i^{-1}，我们要尝试着拆分 p 使其能被更基本的元素进行表示。

与 i 相关地，我们将其对 i 进行整除，分成商与余数两个部分，即 k=\left\lfloor{\dfrac{p}{i}}\right\rfloor，j(p\bmod i)。这时候我们容易用 k.j 表示出 p，也就是 p=ki+j，对于模 p 意义下的情况，表示为 ki+j\equiv 0\pmod p。

我们容易在 i,j 同时存在时将其转化为逆元的表示方法，具体操作是在两边同时乘上 i^{-1}\times j^{-1}，对于底数相同的部分抵消后就得到 kj^{-1}+i^{-1}\equiv 0\pmod p。然后我们移项，利用主元法将 i^{-1} ；留在左边，并且代入原来设元确认的 j,k 值，即：

i^{-1}\equiv -\left\lfloor{\dfrac{p}{i}}\right\rfloor(p\bmod i)^{-1}\pmod p

显然，这时候我们要求出当前处理数 i（已知）时，唯一需要的值便是 p\bmod i 的逆元。而不难发现该数是一定小于 i 的（即 p\bmod i<i），那么我们可以认为该逆元是已知的，而前方的系数也容易求得。

于是我们得到了通式，即对于 i=1 的情况逆元为 1，而对于 i>1 的情况我们可以得到 i^{-1}\equiv -\left\lfloor{\dfrac{p}{i}}\right\rfloor(p\bmod i)^{-1}\pmod p。

由于一切在模 p 意义下进行，所以存在 -\left\lfloor{\dfrac{p}{i}}\right\rfloor=p-\left\lfloor{\dfrac{p}{i}}\right\rfloor，通过这样我们可以避免负数的问题。

#include<bits/stdc++.h>
#define int long long
using namespace std;
int n,p,inv[20000530];
signed main(){
    cin>>n>>p;
    inv[1]=1;
    for(int i=2;i<=n;++i)
        inv[i]=(p-p/i)*inv[p%i]%p;
    for(int i=1;i<=n;i++)
        cout<<inv[i]<<'\n';
    return 0;
}

然而实现方法并不止这一种，在我们得到 kj^{-1}+i^{-1}\equiv 0\pmod p 时就可以递归求得 i^{-1} 的解，为了防止重复递归，我们可以使用记忆化方法，这样依然能够保证每个逆元只被求取一遍，复杂度仍旧是 O(n)。

至此，我们完成了线性求逆元的过程，即 O(n) 求取 1\sim n 的逆元。

线性求任意逆元

其实要求任意给定的 n 个数的逆元也很简单。

给定 n 个数 a_i 满足 1\le a_i<p，我们求出来 n 个数的前缀积，也就是说对于任意的 i\in[1,n]，存在：

s_i=\prod\limits_{x=1}^{n}a_i

这时候我们可以用 Exgcd 或者快速幂法，在 \log p 时间内求出来 s_n 的乘法逆元，记为 s_{inv_n}。不难发现 (a\times b)^{-1}=a^{-1}\times b^{-1}，所以我们将 s_n 乘以 a_n 的时候，a_n 就会和其逆元抵消掉。也就是说，我们可以通过这种方式求出 s_{inv_{n-1}}。

这样我们很容易依次求出所有前缀积的逆元 s_{inv_{i}}，而这让我们可以通过将前缀积的逆元乘上 i-1 的前缀积以消除所有其他的逆元，这样我们就得到了 a_i 的逆元，也就是说 s_{i-1}\times sv_{i}=a_i^{-1}。

代码实现：

s[0]=1;
for(int i=1;i<=n;++i)s[i]=s[i-1]*a[i]%p;
sv[n]=qpow(s[n],p-2);
for(int i=n;i>=1;--i)sv[i-1]=sv[i]*a[i]%p;
for(int i=1;i<=n;++i)inv[i]=sv[i]*s[i-1]%p;

通过这样我们就在 O(n+\log p) 这一接近线性的时间内完成了任意逆元的求取。

Part 6 例题求解

模板

显然可以用递推法实现。

#include<bits/stdc++.h>
#define int long long
using namespace std;
int n,p,inv[20000530]={0,1};
signed main(){
    cin>>n>>p;
    for(int i=2;i<=n;++i)inv[i]=(p-p/i)*inv[p%i]%p;
    for(int i=1;i<=n;i++)cout<<inv[i]<<endl;
    return 0;
}

同余方程

显然可以用快速幂法进行求解，我们只需要将 p-1 替换成 \varphi(b)-1 就好啦。但是我们注意到 b 的范围很大，所以不可能通过 O(b) 方法进行筛除，我们考虑用单点求 \varphi。同时辅助快速幂进行实现即可。

#include<bits/stdc++.h>
#define int long long
using namespace std;
int a,b;
int phi(int n){
    int ans=1;
    for(int i=2;i*i<=n;i++){
        if(n%i==0){
            n/=i;ans*=(i-1);
            while(n%i==0)n/=i,ans*=i;
        }
    }if(n>1)ans*=(n-1);
    return ans;
}int ksm(int base,int p){
    int ans=1;
    while(p){
        if(p&1)ans=ans*base%b;
        base=base*base%b;
        p>>=1;
    }return ans%b;
}signed main(){
    cin>>a>>b;
    int pb=phi(b);
    cout<<ksm(a,pb-1)%b<<endl;
    return 0;
}

希望对你有用，点个赞吧。