扩展欧几里得和乘法逆元求解

欧几里得算法

为方便后续扩展欧几里得的证明，这里简单回顾一下欧几里得算法的证明

假设欲求整数$a,b$的最大公因子$d=gcd(a,b)$

使用带余除法，则b除a表示为 $a=qb+r,\ (0\leq r< b)$
显然若$r=0$，则$gcd(a,b)=b$
若$r\neq0$，因为 $d|a$ 且 $d|b$ ，所以一定有 $d|r$，接下来用反证法证明$gcd(b,r)=d$
假设$gcd(b,r)=c$且$c>d$，则 $c|b$ 且 $c|(qb+r)=a$，与$gcd(a,b)=d$矛盾，因此$gcd(b,r)=d$

综上即有欧几里得公式$gcd(a,b)=gcd(b,r)=gcd(b,a\mod b)$

扩展欧几里得与ax+by=gcd(a,b)

递归求解

假设有 $ax+by=gcd(a,b)$ 和 $bx’+(a\mod b)y’=gcd(b,a\mod b)$ 成立
则根据欧几里得公式有 $ax+by=bx’+(a\mod b)y’$
将$a\mod b=a-\lfloor\frac{a}{b}\rfloor \times b$带入并整理得 $ax+by=ay’+b*(x’- \lfloor a/b \rfloor y’)$
由此我们得出递归关系

$$\begin{cases} x=y'\quad \\ y=x'-\lfloor a/b \rfloor y'\quad \end{cases}$$

递归边界为$b=0$时$x=1,y=0$，即$a1+b0=gcd(a,b)=a$

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

非递归求解

我们回顾欧几里得算法得推到序列

$$\begin{aligned} &a=q_1b+r_1 && ——gcd(a,b)\\ &b=q_2r_1+r2 && ——gcd(b,r_1)\\ &r_1=q_3r_2+r3 && ——gcd(r_1,r_2)\\ & ...\\ &r_i=q_{i+2}r_{i+1}+r_{i+2} && ——gcd(r_i,r_{i+1})\\ & ...\\ &r_{n-1}=q_{n+1}r_n+0 \end{aligned}$$

其中$r_n$就是$d=gcd(a,b)$

我们令$r_{-1}=a,r_0=b$，则可得

$$\begin{cases} r_i=r_{i-2}-q_ir_{i-1}\\ q_i=\lfloor\frac{r_{i-2}}{r_{i-1}}\rfloor \end{cases}$$

假设$r_i=ax_i+by_i$，代入上式可得

$$\begin{cases} x_i=x_{i-2}-q_ix_{i-1}\\ y_i=y_{i-2}-q_iy_{i-1} \end{cases}$$

显然该式初始值为 $x{-1}=1,y{-1}=0$ 和 $x_0=0,y_0=1$

def exgcd(a, b):
    m, n = 0, 1
    x, y = 1, 0
    while b != 0:
        d = a // b
        m, x = x - d * m, m
        n, y = y - d * n, n
        a, b = b, a % b

    return x, y

ax+by=gcd(a,b)的全部解

上述算法仅仅求出了方程的一组解，但其实答案远远不止一组

设当前解为$x,y$，设下一组解为$x+s1,y+s_2$
那么有$a(x+s{1})+b(y+s{2})=gcd(a,b)=ax+by$，化简得$as{1}=-bs_{2}$，即$\frac{s_1}{s_2}=-\frac{b}{a}$

这样我们就得到了所有解的关系式

$$\begin{cases} x'=x+\frac {b}{gcd(a,b)}\times k\\ y'=y-\frac {a}{gcd(a,b)}\times k \end{cases},\quad k\in Z$$

特别的，原方程的最小非负整数解为$x’=(x\bmod \frac {b}{gcd(a,b)}+\frac {b}{gcd(a,b)})\bmod \frac {b}{gcd(a,b)}$

ax+by=c与ax≡c(mod b)求解

假设有$ax+by=gcd$成立
令等式两边同时乘以$\frac {c}{gcd}$，则有$a\frac {cx}{gcd}+b\frac {cy}{gcd}=c$成立
因此可以先用exgcd求出$ax+by=gcd$的一组解$x,y$，则$x=\frac {cx{1}}{gcd},y=\frac {cy{1}}{gcd}$就是$ax+by=c$的一组解

但是要注意根据裴蜀定理，这样的方程存在解的必要条件为$gcd(a,b)|c$

同样用与上述类似的方法推导出所有解递推式

$$\begin{cases} x'=\frac {cx}{gcd}+\frac {b}{gcd(a,b)}*k \\ y'=\frac {cy}{gcd}-\frac {a}{gcd(a,b)}*k \end{cases},\quad k\in Z$$

原方程的最小非负整数解的关系式为 $x’=(\frac {cx}{gcd}\bmod \frac {b}{gcd(a,b)}+\frac {b}{gcd(a,b)})\bmod \frac {b}{gcd(a,b)}$

对于$ax\equiv c(\mod b)$，其实质上等价于求解$ax+by=c$
但是有许多解在模b意义下是相同的，因此上面所有解的递推式中$k\in[0,gcd(a,b)-1]$是就是所有解

乘法逆元的求解

若整数$a$与$n$互素，则满足$ab\equiv 1(\mod n)$的$b$为$a$在模$n$意义下的乘法逆元

显然乘法逆元的定义可以转化为式$ax+by=1$并用扩展欧几里得求解

​```cpp
int exgcd(int a,int b,int &x,int &y)
{
    if(b==0){x=1;y=0;return a;}
    int gcd=exgcd(b,a%b,x,y);
    int tp=x;
    x=y; y=tp-a/b*y;
    return gcd;
}
int x,y;
exgcd(a,b,x,y);
inv=(x%b+b)%b;
​```

乘法逆元也可以根据费马小定理求解

若p为素数， a为正整数，且a、p互质，则有$a^{p-1}\equiv1(\mod p)$

因此$x=a^{p-2}\mod p$就是$a$在模$p$意义下的乘法逆元，该方法仅在模数为素数时可用

此外还有一种线性求多个乘法逆元的算法

首先显然有$1^{-1}≡1(\mod p)$
设$p=ki+r$ $(r<i$, $1<i<p)$，则有 $ki+r≡0(\mod p)$
两边同时乘上 $i^{-1}+r^{-1}$ 得

$$k*r^{-1}+i^{-1}≡0(\mod p) \\ i^{-1}≡-k*r^{-1}(\mod p) \\ i^{-1}≡-\lfloor p/i \rfloor*(p \mod i)^{-1}(\mod p)$$

如下为其代码，显然算法复杂度为O(n)

inv[1]=1;
for(int i=2;i<=n;++i)
inv[i]=(ll)(p-p/i)*inv[p%i]%p;