扩展欧几里得相关知识

扩展欧几里得是干什么的

扩展欧几里得算法是欧几里得算法的扩展，可以求出方程 $ax+by=\gcd (a,b)$ 的解，并可用于多种算法。

怎么做扩展欧几里得

在此之前

欧几里得算法（又称辗转相除法） 是一种可以高效求两个数的最大公约数的算法。它的原理不是这篇文章的主要内容，这里不赘述。总之代码如下：

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

原理一句话：$\gcd (a,b)=\gcd (b,a\mathrm{mod}b)$

扩展欧几里得

那么对于方程 $ax+by=\gcd (a,b)$ ，我们根据裴蜀定理知道这样的 $x,y$ 必定存在，所以我们将式子变为

$$b\times x1+(a\mathrm{mod}b)\times y1=gcd(a,b)$$
化简为 $a\times y1+b\times (x1-(a÷b)\times y1)=gcd(a,b)$。那么我们就可以得出扩展欧几里得代码：

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

这里的边界条件和欧几里得算法一致。

模版传送门

题面

题目背景

这是一道模板题

题目描述

给定 $n,p$ 求 $1\sim n$ 中所有整数在模 $p$ 意义下的乘法逆元。

这里 $a$ 模 $p$ 的乘法逆元定义为 $ax\equiv 1(\mathrm{mod}p)$ 的解。

输入格式

一行两个正整数 $n,p$。

输出格式

输出 $n$ 行，第 $i$ 行表示 $i$ 在模 $p$ 下的乘法逆元。

输入输出样例 #1

输入 #1

10 13

输出 #1

1
7
9
10
8
11
2
5
3
4

说明/提示

$ 1 \leq n \leq 3 \times 10 ^ 6$，$n < p < 20000528 $。

输入保证 $ p $ 为质数。

思路

对于逆元的求法，即 $ax\equiv 1(\mathrm{mod}p)$，这里 $a$ 为要求逆元的数， $x$ 为 $a$ 的逆元。因为 $p$ 是质数，所以 $x$ 必定存在。我们将原式化为

$$ax+kp\equiv 1(\mathrm{mod}p)$$

这是一个扩展欧几里得方程，其中 $a,p$ 已经确定，所以我们直接求 $x,k$ 即可。最终答案为 $x$ 。

#include<iostream>
using namespace std;
int n,p;
void exgcd(int now,int mod,int &x,int &y)
{
    if(!mod)
    {
        x = now;
        y = 1;
        return;
    }
    exgcd(mod,now % mod,x,y);
    int tmp = x;
    x = y;
    y = tmp - y * (now / mod);
}
int main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);
    cin >> n >> p;
    for(int i = 1;i <= n;i++)
    {
        int x = 0,y = 0;
        exgcd(i,p,x,y);
        cout << (x + p) % p << '\n';
    }
    
}

可惜，这道题的时限是 $500ms$ ，所以这样做会超时。

实际上，求一段区间数的逆元有一种非常简单的递推方法，这里不赘述它的原理。设 $in{v}_i$ 表示 $i$ 的逆元，转移方程为

$$in{v}_i=(p-p÷i)\times in{v}_{p\mathrm{mod}i}\mathrm{mod}p$$

#include<iostream>
using namespace std;
int n;
long long p;
long long inv[3000006];
int main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);
    cin >> n >> p;
    inv[1] = 1;
    cout << 1 << '\n';
    for(int i = 2;i <= n;i++)
    {
        inv[i] = ((p - p / i) * inv[p % i]) % p;
        cout << inv[i] << '\n';
    }
    
}