exCRT(扩展中国剩余定理)学习笔记

exCRT(扩展中国剩余定理)学习笔记

背景(?

普通的 CRT 不能解决模数不互质的情况,按照数论算法的命名规则,能解决模数互质的算法就是 exCRT。

然鹅 exCRT 与 CRT 算法本身没有任何关系,个人觉得 exCRT 代码还更简单。

所以 CRT 有什么用?名字里有中国吗?

前置芝士

  • exgcd(扩展欧几里得);
  • 一个简单的同余式转化,$a \equiv b \mod p \Rightarrow a + kp = b$;
  • 脑子

问题

link

解决这样的一个线性同余方程,其中 $a$ 没有特别约束。

n = 2

从简单入手我们可以先考虑 $n = 2$ 的情况:

即 $ x = k_1 b_1 + a_1 = k_2b_2 + a_2 $,得到:

这就是很显的 exgcd 了。用 exgcd 求出 $k_1$,即得 $x_0 = k_1b_1+a_1$,这个方程是否有解即可得到原方程是否有解。

$x_0$ 是一个特解,因为第一个同余方程周期是 $b_1$,第二个是 $b_2$,所以显然解的周期就是 $lcm (b_1, b_2)$,即我们得到了解的表达式:

n > 2

$ n = 2 $ 时的做法相当于合并了两个同余方程,所以我们只要按照这个做法一直合并合并知道只剩最后一个方程为止。

CODE

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
#include <bits/stdc++.h>
using namespace std;

#define int long long

inline int read() {
int x = 0, f = 0; char c = 0;
while (!isdigit(c)) f |= c == '-', c = getchar();
while (isdigit(c)) x = (x << 3) + (x << 1) + (c ^ 48), c = getchar();
return f ? -x : x;
}

#define N 100010

int n, a[N], b[N], P;

int Mul(int x, int k, int r = 0) {
for (; k; k >>= 1, x = (x + x) % P) {
if (k & 1) r = (r + x) % P;
}
return r;
}

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

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

int excrt() {
int A = a[1], B = b[1], k, y;
// A : a1, B : b1, k : k1
for (int i = 2; i <= n; i ++) {
// a[i] : a2, b[i] : b2
int d = gcd(B, b[i]), t = (a[i] - A % b[i] + b[i]) % b[i];
// t : a2 - a1, 在模 b2 意义下
exgcd(B, b[i], k, y);
if (t % d) return -1;
P = b[i] / d;
A += Mul(k, t / d) * B; // a = k1b1 + a1
B = B / d * b[i]; // b = lcm(b1, b2)
A = (A + B) % B;
}
return A;
}

signed main() {
n = read();
for (int i = 1; i <= n; i ++) {
b[i] = read(), a[i] = read();
}
printf("%lld\n", excrt());
return 0;
}

本博客所有文章除特别声明外,均采用 CC BY-SA 4.0 协议 ,转载请注明出处!