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

Aw顿顿

2021-05-14 21:12:57

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

这个式子的意义是,我们可以缩减 xy 的范围,递归求解了。

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,那么我们简单地称 xa 的逆元,记为 a^{-1}。对于我们目前研究的整数系统,存在最为常见的对称运算是加法与乘法。对于实数 x 来说,它的乘法逆元就是 \frac{1}{x},这也是我们记作 x^{-1} 的原因之一。

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

模意义下的乘法逆元

在模 p 意义下,如果 ax=1,那么 xa 的乘法逆元,记为 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\rfloorj(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;
}

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