【数论】gcd、exgcd、质数筛法、逆元与同余方程组复习（已废弃）

前言

这是我面对自己稀烂不堪的数论后奋笔疾书写出来的一篇复习笔记。算是对我前几个月学习过的初等数论的一次总结整理，一个下午的时间可能不够全面，往各位大佬海涵。

最大公约数

这里我们使用欧几里得算法（也称为辗转相除法）来计算最大公约数——简称“gcd”。因为是复习笔记，这里就不详细讲解原理了，记住公式即可。

算法流程： $gcd( m , n ) = gcd( n , m \% n)$

那么直到最后变成 $gcd(m , 0)=m$ ,m 最后的取值也就是 m 和 n 的最大公约数。

用于计算 gcd(m,n)的过程：
如果 n=0，返回值 m 则作为结果。通过计算过程结束，否则进入第二步。
第二步：m 除以 n，将余数赋值给 r；
第三步：将 n 的值赋给 m，将 r 的值赋值给 n。返回第一步

我们通过三元运算符?:很容易写出函数 gcd：

1
int gcd(int m, int n)
2
{
3
    return n == 0 ? m : gcd(n, m % n);
4
    // 翻译如下
5
    // if(n == 0) return m;
6
    // return gcd(n, m % n);
7
}

例题 P2441

非常水的一道绿题。

题目链接：P2441 角色属性树 - 洛谷

1
#include <iostream>
2
#include <cstring>
3
#include <cstdio>
4
#include <cmath>
5
#define N 200005
6

7
using namespace std;
8

9
int fa[N]; // father
10
int val[N];
11
int n, k;
12

13
int gcd(int m, int n)
14
{
15
    return n == 0 ? m : gcd(n, m % n);
16
}
17

18
int dfs(int crd, int w)
19
{
20
    int d = -1;
21
    int v = fa[crd];
22
    if(!v) return -1;
23
    if(gcd(w, val[v]) > 1) d = v;
24
    else d = dfs(v, w);
25
    return d;
26
}
27

28
int main()
29
{
30
    scanf("%d%d", &n, &k);
31
    for(int i = 1; i <= n; i++) scanf("%d", &val[i]);
32
    for(int i = 1; i < n; i++)
33
    {
34
        int u, v;
35
        scanf("%d%d", &u, &v);
36
        fa[v] = u;
37
    }
38
    for(int i = 1; i <= k; i++)
39
    {
40
        int u, crd;
41
        scanf("%d%d", &u, &crd);
42
        if(u == 1) printf("%d\n", dfs(crd, val[crd]));
43
        else
44
        {
45
            int a;
46
            scanf("%d", &a);
47
            val[crd] = a;
48
        }
49
    }
50
    return 0;
51
}

扩展欧几里得算法

假设我们现在有一个方程 $am + bn = gcd(m, n)$ ，求 a，b 的解。

计算特解

扩展欧几里得算法可以解决这个问题——“简称 exgcd”。我们在 gcd 函数中添加另外两个参数 a 和 b，代表方程中的 a 和 b。

算法流程是这样的：
先按正常方法算出 gcd(m, n)
回溯时遵从以下公式从回溯时开始，初始化 a = 1; b= 0;
a = b;
b = a - b * (m / n);
最后返回 gcd(m, n)

因为我们计算时 a 与 b 相互引用，加上最后回溯时返回最大公约数，我们定义两个局部变量 c 与 x 当工具人协助我们完成算法。同时为了获取此方程的一个特解，我们把 a 与 b 引用(&)。

代码如下：

1
int exgcd(int m, int n, int& a, int& b) // 回溯时改变a与b的原值
2
{
3
    if(n == 0) { a = 1; b = 0; return m; } // n = 0计算结束，回溯开始
4
    int d = exgcd(n , m % n, a, b); // 工具人d存值
5
    int c = a; a = b; b = c - b * (m / n); // 工具人c
6
    return d;
7
}

计算通解

这时候我们得到的 a 和 b 只不过是一组特解，我们可以通过 d 写出通解：

$\left(a\:+\:\frac{b}{d}k,\:b\:+\:\frac{a}{d}k\right)\:k\epsilon\mathbb{Z}$

扩展情况

假设要求方程 $am + bn = g$ ，求 a，b 的解。

我们可以先按原来的方法，求出 $am + bn = gcd(m, n)$ 的解。

然后把 a，b，gcd(m, n)同时乘上 g / gcd(m, n)就可以了。注意 $gcd(m, n) | g$ 。

质数筛法

很蠢的埃式筛我们就不扯了，它浪费时间的地方就是重复筛——比如 15 被 3 和 5 筛了两次。

我们直接拿出欧拉筛，规定每个数字只被它最小的质因数筛去，复杂度降为 $O(N)$ 。

直接上代码：

1
#define N 100005
2
bool not_prime[N];
3
int prime[N], tot;
4
for(int i = 2; i <= N; i++)
5
{
6
    if(!not_prime[i]) primt[++tot] = i;
7
    for(int j = 1; j <= n && i * prime[j] <= N; j++)
8
    {
9
        not_prime[i * prime[j]] = true;
10
        if(i % prime[j] == 0) break;
11
    }
12
}

逆元

逆元是一个很神奇的东西，你可以理解成广义上的倒数。

什么是倒数？任意一个数乘上一个倒数，结果为 1，相当于除了它自己。那么举个例子，有些题目要求你对答案取模，众所周知乘法加法一边运算一边取模是不会影响正确性的，但是要是来个除，你就不得不把它计算完全之后再取模——有可能爆 long long。

那么逆元就出现了！逆元就是一个数在一个取模环境下的倒数，乘上它就等于除以原来的数字，就与加法乘法一样了。我们这样定义逆元：

假设 a，b 满足以下条件（模 m 意义下的等价类集合）

$a,b\epsilon\mathbb{Z}\_m$

同时满足 $ab\:\equiv\:1\:(\mod\:m)$ ，我们就说 b 是 a 的逆元，或者 a 是 b 的逆元。是不是很简单？

逆元可不只有一个哦！

求解逆元

求解逆元有很多方法，这里我们就只联系上下文，讲一下如何通过 exgcd 求解逆元。

我们简单想一想，假设我们要求 a 的逆元 b，如果 $ab \mod m = 1$ ，是不是说明 $ab - mk = 1$ 呢？而这里我们的a和m是已知的，我去，这不就是 $am + bn = g$ 的形式吗！你说一个减一个加？这有什么关系？不就相当于 $am + (-k)m = 1$ 吗？其实是没有区别的。

那就轻轻松松了，结合我们之前 exgcd 的铺垫，我们直接使用exgcd(a, m, b, k)运行完毕之后b/d就是 a 的一个逆元了。~~当然如果b/d除不尽那就寄啦~~。

但是我说了，逆元并不只有一个，要是题目要求是最小正整数的逆元怎么办？

直接上公式：a的逆元 = (b % m + m) % m;至于为什么，请读者自行证明吧（毕竟是复习笔记）。

卢卡斯定理/Lucas 定理

假设有一道题目让你求 ${C}_{n+m}^{n}\mod p$ 的值，你会怎么求？

也许你会不屑一顾——我逆元在手，不轻轻松松？

但是你忘记了一件事……如果有一部分刚好是 p 的倍数，那不就寄啦？

呵呵，这时候就要请出我们的 Lucas 定理了，专门来解决这种组合数取模的题目。

实现

假设 A、B 是非负整数，p 是质数。AB 写成 p 进制： $A=a[n]a[n-1]...a[0]$ ， $B=b[n]b[n-1]...b[0]$ 。则组合数 $C(A,B)$ 与 $C(a[n],b[n])*C(a[n-1],b[n-1])..C(a[0],b[0]) \mod p$ 同余。即：

$Lucas(n,m,p)=C(n\%p,m\%p)*Lucas(n/p,m/p,p)$

简言之，Lucas 定理是用来求 $C(n,m) \mod p$ ，p 为素数的值。

例题 P3807

题目链接：P3807 【模板】卢卡斯定理/Lucas 定理

1
#include<iostream>
2
#include<algorithm>
3
#include<cstdio>
4
#define N 100005
5

6
using namespace std;
7

8
int t, p;
9
int n, m;
10
long long a[N], b[N];
11

12
long long lucas(int x,int y)
13
{
14
  if(x < y) return 0;
15
  else if(x < p) return b[x] * a[y] * a[x - y] % p;
16
  else return lucas(x / p, y / p) * lucas(x % p,y % p) % p;
17
}
18

19
int main()
20
{
21
  scanf("%d", &t);
22
  while(t--)
23
  {
24
    scanf("%d%d%d", &n, &m, &p);
25
    a[0] = a[1] = b[0] = b[1] = 1;
26
    for(int i = 2; i <= n + m; i++) b[i] = b[i - 1] * i % p; // 处理成p进制
27
    for(int i = 2; i <= n + m; i++) a[i] = (p - p / i) * a[p % i] % p;
28
    for(int i = 2; i <= n + m; i++) a[i] = a[i - 1] * a[i] % p;
29
    printf("%lld\n", lucas(n + m, m));
30
  }
31
  return 0;
32
}

同余方程组

几个 $a\:\equiv\:b\:(\mod\:m)$ 叠在一起就变成了同余方程组。

例如下面的同余方程组：

$\begin{cases}{a\:\equiv\:b_1\:(\mod\:m_1)} \\ {a\:\equiv\:b_2\:(\mod\:m_2)} \\ {a\:\equiv\:b_3\:(\mod\:m_3)}\end{cases}$

我承认这是一个很吓人的东西，但是推导也不难。

虽然可以用中国剩余定理做，但是可惜它只能解决模数互质的情况。

那么我们可以换一种思路，合并同余方程组。

~~可惜过程我忘光了呜呜呜，老师救我数论~~ 有时间再补上推理吧——2023/07/05

直接上代码：

注意暴力计算会超 long long，这里我用了龟速幂等方法成功把其压在了 long long 内。嫌弃麻烦的同学可以直接使用__int128 或者高精度，这里不再赘述。~~我故意保留了一部分注释，让你知道这是我打的~~

1
#include <iostream>
2
#include <cmath>
3
#define ll long long
4

5
using namespace std;
6

7
ll Exgcd(ll a, ll b, ll &x, ll &y)
8
{
9
    if (b == 0)
10
    {
11
        x = 1, y = 0;
12
        return a;
13
    }
14
    ll d = Exgcd(b, a % b, x, y);
15
    ll z = x; x = y; y = z - (a / b) * y;
16

17
    return d;
18
}
19

20
ll qmul(ll a, ll b, ll p)
21
{
22
    a %= p;
23
    b = (b % p + p) % p;
24
    ll res = 0;
25
    while (b)
26
    {
27
        if (b & 1) { res = (res + a) % p; }
28
        a = (a + a) % p;
29
        b >>= 1;
30
    }
31
    return res;
32
}
33

34
int main()
35
{
36
    ll a1, a2;
37
    ll b1, b2;
38
    ll n = 0;
39
    //freopen("1.txt","r",stdin);
40
    scanf("%lld", &n);
41
    scanf("%lld%lld", &a1, &b1);
42
    //printf("1:\nX = %lld (mod %lld)\n", b1, a1);
43
    for(ll i = 2; i <= n; i++)
44
    {
45
        scanf("%lld%lld", &a2, &b2);
46
        //printf("2:\nX = %lld (mod %lld)\n", b2, a2);
47
        ll y1 = 0, y2 = 0;
48
        //a1 * y1 - a2 * y2 = b2 - b1
49
        //y1 y2 dont know
50
        ll d = Exgcd(a1, a2, y1, y2);
51

52
        ll cnt = a2 / d;
53
        //y1 = y1 / d;
54
        //y1 = (y1 % cnt + cnt) % cnt;
55
        //y1 = (y1 = y1 / d * (b2 - b1) % cnt + cnt) % cnt;
56
        //y1 = (qmul(b2 - b1, y1, cnt) + cnt) % cnt;
57
        y1 = (qmul(y1, ((b2 - b1) / d % a2 + a2) % a2, cnt) + cnt) % cnt;
58

59
        b1 = y1 * a1 + b1;
60
        a1 = (a1 * cnt);
61
        b1 = (b1 % a1 + a1) % a1;
62
        //printf("AND:\nX = %lld (mod %lld)\n", b1, a1);
63
    }
64

65
    //printf("X = %lld k + %lld \n", a1, b1);
66

67
    ll ans = (b1 % a1 + a1) % a1;
68
    printf("%lld", ans);
69

70
    return 0;
71
}

【数论】gcd、exgcd、质数筛法、逆元与同余方程组复习（已废弃）

前言

最大公约数

例题 P2441

扩展欧几里得算法

计算特解

计算通解

扩展情况

质数筛法

逆元

求解逆元

卢卡斯定理/Lucas 定理

实现

例题 P3807

同余方程组

觉得这篇文章怎么样？

分享这篇文章

评论区

标签

分类

前言

最大公约数

例题 P2441

扩展欧几里得算法

计算特解

计算通解

扩展情况

质数筛法

逆元

求解逆元

卢卡斯定理/Lucas 定理

实现

例题 P3807

同余方程组

觉得这篇文章怎么样？

分享这篇文章

评论区