扩展 CRT
突然发现博客里中国剩余定理的文章少之又少。。。
然后昨天考场上忘了 CRT… 然后就自己 YY 了一个扩展 CRT 发现是对的。。。
(换句话说就是我想了半天想起来了)
于是印象加深了很多,就写篇文章记录一下。。。
首先扩展 CRT 是用来解决模数不互质情况下的模意义一元线性方程组的
思想就是把方程两两合并。
比如:
$x \equiv a _ 1 \text{ mod } p _ 1$
$x \equiv a _ 2 \text{ mod } p _ 2$
其中 $x$是未知数,$p _ 1,p _ 2$不互质
我们按套路转化一下式子:
$x = p _ 1 k _ 1 + a _ 1$
代入第二个式子:
$p _ 1 k _ 1 + a _ 1 \equiv a _ 2 \text{ mod } p _ 2$
移项:
$p _ 1 k _ 1 \equiv a _ 2 – a _ 1 \text{ mod } p _ 2$
OK,是不是很眼熟?
用扩展欧几里得呀!
设 $d = gcd(p _ 1, p _ 2)$
由扩展欧几里得定理可知方程有解当且仅当 $d \mid a _ 2 – a _ 1$,否则无解
然后等号两边同除一个 $d$:
$\frac {p _ 1} {d} k _ 1 \equiv \frac{a _ 2 – a _ 1} {d} \text{ mod } \frac {p _ 2} {d}$
$k _ 1 \equiv \frac {a _ 2 – a _ 1} {d} \times (\frac {p _ 1} {d}) ^ {-1} \text { mod } \frac {p _ 2} {d}$
其中 $(\frac {p _ 1} {d}) ^ {-1}$表示 $\frac {p _ 1} {d}$在模 $\frac {p _ 2} {d}$意义下的逆元,很显然这个逆元就是前面欧几里得求出的 $k _ 1$
设 $k = \frac {a _ 2 – a _ 1} {d} \times (\frac {p _ 1} {d}) ^ {-1}$
则 $k _ 1 = k + \frac {p _ 2} {d} \times i$
回代:
$x = p _ 1 k _ 1 + a _ 1$
$=p _ 1(k + \frac {p _ 2} {d} \times i) + a _ 1$
$=\frac {p _ 1\times p _ 2} {d} \times i + p _ 1\times k + a _ 1$
$\equiv p _ 1 \times k + a _ 1 \text { mod } \frac {p _ 1\times p _ 2} {d}$
另外为了方便,我们其实不需要单独计算 $k$的值,直接用 $k _ 1$代替即可。
因为:
$k \equiv k _ 1 \text { mod } \frac {p2} {d}$
所以:
$p _ 1 \times k \equiv p _ 1 \times k _ 1 \text { mod } \frac {p _ 1 \times p _ 2} {d}$
因此用 $k _ 1 和 k$是等价的
所以合并的式子也可以写成:
$x \equiv x _ 0 \text { mod } lcm(p _ 1, p _ 2)$
$x _ 0$是这组方程的解。
例题
【模板】扩展中国剩余定理(EXCRT)—— LUOGU – 4777
代码:
#include <bits/stdc++.h>
#define NS (100005)
using namespace std;
typedef long long LL;
template <typename _Tp> inline void IN(_Tp& dig)
{
char c; bool flag = 0; dig = 0;
while (c = getchar(), !isdigit(c)) if (c == '-') flag = 1;
while (isdigit(c)) dig = dig * 10 + c - '0', c = getchar();
if (flag) dig = -dig;
}
int n;
LL A[NS], B[NS], M, X;
LL mul(LL a, LL b, LL M)
{
a = a % M, b = b % M;
return ((a * b - (LL)(((long double) a * b + 0.5) / M) * M) % M + M) % M;
}
void exgcd(LL a, LL b, LL& x, LL& y, LL& g)
{
if (!b) {x = 1, y = 0, g = a; return;}
exgcd(b, a % b, y, x, g);
y -= a / b * x;
}
int main(int argc, char const* argv[])
{
IN(n);
for (int i = 1; i <= n; i += 1) IN(A[i]), IN(B[i]);
M = A[1], X = B[1]; LL g, x, y, t;
for (int i = 2; i <= n; i += 1)
{
exgcd(M, A[i], x, y, g);
if ((B[i] - X) % g) puts("No Solution"), exit(0);
t = A[i] / g, x = (mul(x, (B[i] - X) / g % t + t, t) + t) % t;
t = M, M = M / g * A[i], X += mul(t, x, M), X %= M;
}
printf("%lld\n", X);
return 0;
}
4 条评论
mqh · 2020年1月3日 9:40 下午
为什么 k 恒等 k1mod(p2/d) 呢?
不是 k1 恒等 kmod(p2/d) 吗? 上面有 k1=k+(p2/d)*i 呀
可以解释一下吗,谢谢啦
mqh · 2020年1月3日 8:51 下午
p1k1+a1 恒等 a2modp2
移项后不是
p1k1 恒等 -a1+a2modp2 吗?
mqh · 2020年1月3日 9:33 下午
我懂了 ^_^
Remmina · 2020年1月4日 12:31 下午
(* ̄︶ ̄) 那就好~