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

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

EveSunMaple Lv3
  2023-07-04 00:07 2023-07-04 00:07 创建   2024-02-23 12:02:20 2024-02-23 12:02:20 更新      2.6k 字  12 分钟  

前言

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

最大公约数

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

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

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

用于计算gcd(m,n)的过程：

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

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

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

例题P2441

非常水的一道绿题。

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
#include <iostream>
#include <cstring>
#include <cstdio>
#include <cmath>
#define N 200005

using namespace std;

int fa[N]; // father
int val[N];
int n, k;

int gcd(int m, int n)
{
    return n == 0 ? m : gcd(n, m % n);
}

int dfs(int crd, int w)
{
    int d = -1;
    int v = fa[crd];
    if(!v) return -1;
    if(gcd(w, val[v]) > 1) d = v;
    else d = dfs(v, w);
    return d;
}

int main()
{
    scanf("%d%d", &n, &k);
    for(int i = 1; i <= n; i++) scanf("%d", &val[i]);
    for(int i = 1; i < n; i++)
    {
        int u, v;
        scanf("%d%d", &u, &v);
        fa[v] = u;
    }
    for(int i = 1; i <= k; i++)
    {
        int u, crd;
        scanf("%d%d", &u, &crd);
        if(u == 1) printf("%d\n", dfs(crd, val[crd]));
        else 
        {
            int a;
            scanf("%d", &a);
            val[crd] = a;
        }
    }
    return 0;
}

扩展欧几里得算法

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

计算特解

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

算法流程是这样的：

  1. 先按正常方法算出gcd(m, n)
  2. 回溯时遵从以下公式
    从回溯时开始，初始化a = 1; b= 0;
  3. a = b;
  4. b = a - b * (m / n);
  5. 最后返回gcd(m, n)

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

代码如下：

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

计算通解

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

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

扩展情况

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

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

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

质数筛法

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

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

直接上代码：

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

逆元

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

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

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

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

a,bϵZma,b\epsilon\mathbb{Z}_m

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

逆元可不只有一个哦！

求解逆元

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

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

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

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

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

卢卡斯定理/Lucas 定理

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

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

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

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

实现

假设A、B是非负整数，p是质数。AB写成p进制：A=a[n]a[n−1]...a[0]A=a[n]a[n-1]...a[0]，B=b[n]b[n−1]...b[0]B=b[n]b[n-1]...b[0]。则组合数C(A,B)C(A,B)与C(a[n],b[n])∗C(a[n−1],b[n−1])..C(a[0],b[0])mod  pC(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(n,m,p)=C(n\%p,m\%p)*Lucas(n/p,m/p,p)

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

例题P3807

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
#include<iostream>
#include<algorithm>
#include<cstdio>
#define N 100005

using namespace std;

int t, p;
int n, m;
long long a[N], b[N];

long long lucas(int x,int y)
{
    if(x < y) return 0;
    else if(x < p) return b[x] * a[y] * a[x - y] % p;
    else return lucas(x / p, y / p) * lucas(x % p,y % p) % p;
}

int main()
{
    scanf("%d", &t);
    while(t--)
    {
        scanf("%d%d%d", &n, &m, &p);
        a[0] = a[1] = b[0] = b[1] = 1;
        for(int i = 2; i <= n + m; i++) b[i] = b[i - 1] * i % p; // 处理成p进制
        for(int i = 2; i <= n + m; i++) a[i] = (p - p / i) * a[p % i] % p;
        for(int i = 2; i <= n + m; i++) a[i] = a[i - 1] * a[i] % p;
        printf("%lld\n", lucas(n + m, m));
    }
    return 0;
}

同余方程组

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

例如下面的同余方程组：

{a ≡ b1 (mod   m1)a ≡ b2 (mod   m2)a ≡ b3 (mod   m3)\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
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
#include <iostream>
#include <cmath>
#define ll long long

using namespace std;

ll Exgcd(ll a, ll b, ll &x, ll &y)
{
    if (b == 0) 
    {
        x = 1, y = 0;
        return a;
    }
    ll d = Exgcd(b, a % b, x, y);
    ll z = x; x = y; y = z - (a / b) * y;

    return d;
}

ll qmul(ll a, ll b, ll p)
{
    a %= p;
    b = (b % p + p) % p;
    ll res = 0;
	while (b)
	{
        if (b & 1) { res = (res + a) % p; }
        a = (a + a) % p;
        b >>= 1;
    }
    return res;
}

int main()
{
    ll a1, a2;
    ll b1, b2;
    ll n = 0;
    //freopen("1.txt","r",stdin);
    scanf("%lld", &n);
    scanf("%lld%lld", &a1, &b1);
    //printf("1:\nX = %lld (mod %lld)\n", b1, a1);
    for(ll i = 2; i <= n; i++)
    {
        scanf("%lld%lld", &a2, &b2);
        //printf("2:\nX = %lld (mod %lld)\n", b2, a2);
        ll y1 = 0, y2 = 0;
        //a1 * y1 - a2 * y2 = b2 - b1
        //y1 y2 dont know
        ll d = Exgcd(a1, a2, y1, y2);
        
        ll cnt = a2 / d;
        //y1 = y1 / d;
        //y1 = (y1 % cnt + cnt) % cnt;
        //y1 = (y1 = y1 / d * (b2 - b1) % cnt + cnt) % cnt;
        //y1 = (qmul(b2 - b1, y1, cnt) + cnt) % cnt;
        y1 = (qmul(y1, ((b2 - b1) / d % a2 + a2) % a2, cnt) + cnt) % cnt;
        
        b1 = y1 * a1 + b1;
        a1 = (a1 * cnt);
        b1 = (b1 % a1 + a1) % a1;
        //printf("AND:\nX = %lld (mod %lld)\n", b1, a1);
    }

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

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

    return 0;
}
  • 标题: 【数论】gcd、exgcd、质数筛法、逆元与同余方程组复习（已废弃）
  • 作者: EveSunMaple
  • 创建于 : 2023-07-04 00:07:00
  • 更新于 : 2024-02-23 12:02:20
  • 链接: https://old.saroprock.com/post/c22a7f0c.html
  • 版权声明: 本文章采用 CC BY-NC-SA 4.0 进行许可。
推荐阅读
【数论】费马小定理 【数论】费马小定理 【题解】7月11到7月17典例复习 【题解】7月11到7月17典例复习 【线性代数】矩阵乘法与线性DP优化 【线性代数】矩阵乘法与线性DP优化
推荐阅读
【数论】费马小定理 【数论】费马小定理 【题解】7月11到7月17典例复习 【题解】7月11到7月17典例复习
 评论
此页目录
【数论】gcd、exgcd、质数筛法、逆元与同余方程组复习（已废弃）
  1. 1. 前言
  2. 2. 最大公约数
    1. 2.1. 例题P2441
  3. 3. 扩展欧几里得算法
    1. 3.1. 计算特解
    2. 3.2. 计算通解
    3. 3.3. 扩展情况
  4. 4. 质数筛法
  5. 5. 逆元
    1. 5.1. 求解逆元
    2. 5.2. 卢卡斯定理/Lucas 定理
  6. 6. 同余方程组