扩展 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;
}
分类: 文章

XZYQvQ

炒鸡辣鸡的制杖蒟蒻一枚QvQ

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 下午

    我懂了 ^_^

发表回复

Avatar placeholder

您的电子邮箱地址不会被公开。 必填项已用 * 标注