中国剩余定理 与 扩展中国剩余定理

中国剩余定理 与 扩展中国剩余定理

一、问题引入

先看一道poj上的题目:POJ1006Biorhythms

题意
  人自出生起就有体力,情感和智力三个生理周期,分别为232833天。一个周期内有一天为峰值,在这一天,人在对应的方面(体力,情感或智力)表现最好。通常这三个周期的峰值不会是同一天。现在给出三个日期,分别对应于体力,情感,智力出现峰值的日期。然后再给出一个 起始日期d,要求从这一天开始,算出最少再过多少天后三个峰值同时出现。

分析
  首先我们知道,任意两个峰值之间一定相距整数倍的周期。假设一年的第N天达到峰值,则 下次 达到峰值的时间为N+kT(T是周期,k是任意正整数)。

所以,三个峰值同时出现的那一天(S)应满足

{S=N1+T1k1S=N2+T2k2S=N3+T3k3

N1,N2,N3分别为为体力,情感,智力出现峰值的日期, T1,T2,T3 分别为体力,情感,智力周期。 需要求出k1,k2,k3三个非负整数使上面的等式成立。

想直接求出k1,k2,k3很难,但是目的是求出S, 可以考虑 从结果逆推。根据上面的等式,S 满足三个要求:

  • 除以T1余数为N1
  • 除以T2余数为N2
  • 除以T3余数为N3

找到最小的一个这样的数就可以了。

这就是著名的 中国剩余定理,老祖宗在几千年前已经对这个问题想出了一个精妙的解法。依据此解法的算法,时间复杂度可达到O(1)

二、中国剩余定理

在《孙子算经》中有这样一个问题:“今有物不知其数,三三数之剩二(除以32),五五数之剩三(除以53),七七数之剩二(除以72),问物几何?”这个问题称为 孙子问题 ,该问题的一般解法国际上称为 中国剩余定理

解法共三步

  • ① 找三个数:

    • 35的公倍数中找出被7除余1的最小数15
    • 37的公倍数中找出被5除余1的最小数21
    • 57的公倍数中找出被3除余1的最小数70
  • ② 找符合条件的一个数:

    • 15乘以22为最终结果除以7的余数)
    • 21乘以33为最终结果除以5的余数)
    • 70乘以22为最终结果除以3的余数)
      然后把三个乘积相加152+213+702=233
  • ③ 用233除以357三个数的最小公倍数105,得到余数,即233%105=23

答:23就是符合条件的最小数

Q:这有数学依据吗?
现在将 孙子问题 拆分成几个简单的小问题,从零开始,揣测古人是如何推导出这个解法的:

首先,我们假设:

  • n1是满足除以32的一个数,比如258等等,也就是满足3k1+2k1>=0 的一个任意数。
  • n2是满足除以53的一个数,比如3813等等,也就是满足5k2+3k2>=0 的一个任意数。
  • n3是满足除以72的一个数,比如2916等等,也就是满足7k3+2k3>=0 的一个任意数。

有了前面的假设,我们先从n1这个角度出发,已知n1满足除以32,能不能使得n1+n2的和仍然满足除以32,进而使得n1+n2+n3的和仍然满足除以32

 这就牵涉到一个 最基本数学定理,如果有a%b=c,则有(a+kb)%b=c(k为非零整数),换句话说,如果一个除法运算的余数为c,那么被除数与k倍的除数相加(或相减)的和(差)再与除数相除,余数不变。

以此定理为依据,如果n23的倍数,n1+n2就依然满足除以32。同理,如果n3也是3的倍数,那么n1+n2+n3的和就满足除以32。这是从n1的角度考虑的,再从n2n3的角度出发,我们可推导出以下三点:

  • 为使n1+n2+n3的和满足除以32n2n3必须是3的倍数

  • 为使n1+n2+n3的和满足除以53n1n3必须是5的倍数

  • 为使n1+n2+n3的和满足除以72n1n2必须是7的倍数

因此,为使n1+n2+n3的和作为 孙子问题 的一个最终解,需满足:

  • n1除以32,且是57的公倍数
  • n2除以53,且是37的公倍数
  • n3除以72,且是35的公倍数

所以,孙子问题解法的 本质 是:

  • 57的公倍数中找一个除以32的数n1
  • 37的公倍数中找一个除以53的数n2
  • 35的公倍数中找一个除以72的数n3

再将三个数相加得到解。在求n1n2n3时又用了一个 小技巧,以n1为例,并非从57的公倍数中直接找一个除以32的数,而是

先找一个 除以31的数,再乘以2。也就是先求出57的公倍数模3下的逆元 ,再用逆元去乘余数!

Q:为啥这样干呢?
:不这样干不好整啊,57的公倍数可不少,一个个暴力去找太累了。不这样暴力还能咋办?
(a / b) % p = (a%p / b%p) %p (错)
因为除法没法进行取模运算,需要找到乘法逆元。
Q:哪里来的除法,哪里来的取模?
:Mi=M/mi就是除法,需要对mi取模。

这里又有一个数学公式,如果a%b=c,那么(ak)%b=a%b+a%b++a%b=c+c++c=kc(k>0),也就是说,如果一个除法的余数为c,那么被除数的k倍与除数相除的余数为kc,展开式中已证明。

最后,我们还要清楚一点,n1+n2+n3只是问题的一个解,并不是最小的解。如何得到最小解?我们只需要从中最大限度的减掉357的公倍数105即可。道理就是前面讲过的定理 如果a%b=c,则有(akb)%b=c ,所以(n1+n2+n3)%105就是最终的最小解。

【中国剩余定理】

xZ,有

{xm1(mod a1)xm2(mod a2)xm3(mod a3)...xmn(mod an)

令$ \large \left{\begin{matrix}
A=a_1a_2a_3...a_n & \
A_i=A/a_i& \
t_i 是A_i关于a_i 的逆元,即A_i
t_i \equiv 1 (mod \ a_i)
\end{matrix}\right.$

res=i=1nmiAiti

【扩展欧几里得求逆元】

扩展欧几里得算法,顾名思义,是欧几里得算法的扩展,这意味着该算法不仅保留了欧几里得算法的核心——求最大公约数,还具备欧几里得算法没有的功能——求解贝祖等式,即可以求出方程ax+by=c的一组解(x0,y0)

特别的,当 gcd(a,b)=1 ,即a,b互素时,ax+by=1中的解

  • xab的乘法逆元,即ax1(mod b)
  • yba的乘法逆元,即bx1(mod a)

解释:方程ax+by=1写成同余方程,就是:ax1(mod b),也就是逆元的定义式

r = (r % a + a) % a

作用:将余数加上若干个循环节,使得结果落在0a1之间。

举栗子
5个一组划分,最后缺少2

本质:就是多了3个,r其实是3

Q:那就直接r+a=3不就行了吗,(r%a+a)%a是什么意思?

A:(1)如果r>=0,r+a>=a了,这时应该再多分出去一组才对,不能直接r+a ,但 (r%a+a)%a就没有问题了
$ (2)r<0,(r%a+a)%a$也一样正确

结论:公式 (r%a+a)%a可以把负的余数改为正的余数

#include <bits/stdc++.h>

using namespace std;
int main() {
 
    int r = -2, a = 5;
    r = (r % a + a) % a;
    cout << r << endl;
    // 输出3
    return 0;
}

三、扩展中国剩余定理

Q:中国剩余定理要求所有的mi互质,那么如果不互质,怎么求解同余方程组?
答:对于模数不互质的情况,我们可以通过合并相邻两组方程的方式来逐步得到解。

{xm1 (mod a1)xm2 (mod a2)...xmn (mod an)

先计算两个线性同余方程组的解,之后依此类推
① 可以写成x=k1a1+m1
② 可以写成x=k2a2+m2

③=④ k1a1+m1=k2a2+m2 k1a1k2a2=m2m1

根据 裴蜀定理,当
m2m1gcd(a1,a2)的倍数时,方程 ⑤ 有无穷多组整数解

k1a1k2a2=gcd(a1,a2)=d,可以用 拓展欧几里得算法 来求解出一组特解(k1,k2),即exgcd(a1,a2,k1,k2)

如何求解k1a1k2a2=gcd(a1,a2)=d中的通解k1,k2呢?

结论:k1,k2的通解

{k1=k1+ka2dk2=k2+ka1d

证明:
k1的通解为s1,k2的通解为s2,k1的特解为k1,k2的特解为k2

{a1s1a2s2=da1k1a2k2=d

d

{a1ds1a2ds2=1a1dk1a2dk2=1

a1d(s1k1)a2d(s2k2)=0

gcd(a1,a2)=da1da2d互质。

解释:最大公约数嘛,如果除完了还有约数,那还是什么最大公约数~

a2ds1k1,a1ds2k2

s1k1=ka2dk1=k1+ka2d整除式变形

s2k2=ka1dk2=k2+ka1d整除式变形

k1通解代入 ③ 得

x=k1a1+m1=(k1+ka2d)a1+m1=k1a1+m1m1+a2da1a1k

对照原方程x=a1k1+m1,x=a2k2+m2,我们发现,两个方程合并后得到的新方程,其实和原来的两个方程都是差不多的样子,只不过系数发生了变化,而合并后的方程肯定是即符合 ① 双符合 ②, 如果继续迭代下去,迭代n1次后,可以将n个线性同余方程合并为一个方程,最终方程就能满足所有的约束要求。为了能够实现代码循环迭代,那么就可以

{m1=k1a1+m1a1=a2da1k

这样反复循环迭代就可以找出最后一全兼容方程了!

Q:此处为什么要取绝对值呢
A:因为不知道a2d的正负性,我们在原基础上要尽量减多个abs(a2d),使其为正整数且最小。

四、练习题

模数互质

P1495 【模板】中国剩余定理(CRT)/曹冲养猪

Code 中国剩余定理解法

#include <bits/stdc++.h>

using namespace std;
#define int long long
#define endl "\n"
const int N = 20;
typedef __int128 INT128;

int n;     // 同余方程个数
int a[N];  // 除数数组
int m[N];  // 余数数组
int A = 1; // 除数累乘积
int x, y;  // (x,y)是方程 Mi * x + m[i] * y = 1的一个解,x是Mi关于m[i]的逆元
int res;   // 最小整数解
// 扩展欧几里得模板
int exgcd(int a, int b, int &x, int &y) {
    if (!b) {
        x = 1, y = 0;
        return a;
    }
    int d = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return d;
}

// 中国剩余定理模板
void CRT() {
    //  一、预处理
    for (int i = 1; i <= n; i++) {
        cin >> a[i] >> m[i];                // 读入除数数组和余数数组
        m[i] = (m[i] % a[i] + a[i]) % a[i]; // ① 预算余数为正数: r[i]可能为负数,预处理成正数,本题没有这个要求,但考虑到通用模板,加上了这块代码
        A *= a[i];                          // ② 预处理除数连乘积
    }
    // 二、计算
    // 功能: 求 x ≡ r_i mod m_i 方程组的解
    for (int i = 1; i <= n; i++) {
        int Ai = A / a[i];                                 // ①除数连乘积除去当前的除数,得到Mi
        exgcd(Ai, a[i], x, y);                             // ②扩欧求逆元
        x = (x % a[i] + a[i]) % a[i];                      // ③逆元防负数
        res = (res + INT128(m[i] * A / a[i] * x) % A) % A; // ④累加res(看公式) 快速乘,防止LL也在乘法过程中爆炸
    }
}

signed main() {
    cin >> n;
    CRT();
    cout << res << endl;
}

Code 扩展中国剩余定理解法

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define endl "\n"

const int N = 100010;
int a[N], m[N];
int n;
int exgcd(int a, int b, int &x, int &y) {
    if (!b) {
        x = 1, y = 0;
        return a;
    }
    int d = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return d;
}

int exCRT() { // m代表合并多少次,有n个方程,其实m=n-1
    int a1 = a[1], m1 = m[1], a2, m2, k1, k2;

    for (int i = 2; i <= n; i++) { // a1,m1:做为先遣军,准备与后续的方程合并
        a2 = a[i], m2 = m[i];
        int d = exgcd(a1, -a2, k1, k2); // 扩展欧几里得求方程特解,x有用,y无用
        if ((m2 - m1) % d) return -1;   // 拼出来的同余方程无解,继续不下去
        k1 *= (m2 - m1) / d;            // 同比扩大指定倍数
        int t = abs(a2 / d);            // 最小累计单元
        k1 = (k1 % t + t) % t;          // 取得最小正整数解,否则可以会在计算过程中爆long long
        m1 = k1 * a1 + m1;              // 准备下一轮的m[1],a[1]
        a1 = abs(a1 / d * a2);          // 两个顺序不能反
    }
    return (m1 % a1 + a1) % a1; // 取得最小正整数解
}

signed main() {
    cin >> n;
    for (int i = 1; i <= n; i++) cin >> a[i] >> m[i];
    cout << exCRT() << endl;
}

POJ1006Biorhythms

Code 中国剩余定理解法

#include <cstdio>
#include <iostream>
using namespace std;
#define endl "\r"
#define int long long

const int N = 10;

int n = 3;
int a[N] = {0, 23, 28, 33};
int m[N];
int A = 1;
int x, y;
int res;

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

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

void CRT() {
    for (int i = 1; i <= n; i++) {
        m[i] = (m[i] % a[i] + a[i]) % a[i];
        A *= a[i];
    }

    for (int i = 1; i <= n; i++) {
        int Ai = A / a[i];
        exgcd(Ai, a[i], x, y);
        x = (x % a[i] + a[i]) % a[i];
        res = (res + qmul(m[i], A / a[i] * x, A)) % A;
    }
}

int T;
int p, e, i, d;
signed main() {
    while (cin >> p >> e >> i >> d && ~p) {
        T++;
        A = 1;
        res = 0;
        m[1] = p, m[2] = e, m[3] = i;
        CRT();
        if (res <= d) res += A;
        printf("Case %lld: the next triple peak occurs in %lld days.\n", T, res - d);
    }
}

Code 扩展中国剩余定理

#include <cstdio>
#include <iostream>
#include <stdlib.h>
using namespace std;
#define int long long
#define endl "\r"
const int N = 10;

int n = 3;
int a[N] = {0, 23, 28, 33};
int m[N];
int exgcd(int a, int b, int &x, int &y) {
    if (!b) {
        x = 1, y = 0;
        return a;
    }
    int d = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return d;
}

int exCRT() { // m代表合并多少次,有n个方程,其实m=n-1
    int a1 = a[1], m1 = m[1], a2, m2, k1, k2;

    for (int i = 2; i <= n; i++) { // a1,m1:做为先遣军,准备与后续的方程合并
        a2 = a[i], m2 = m[i];
        int d = exgcd(a1, -a2, k1, k2); // 扩展欧几里得求方程特解,x有用,y无用
        if ((m2 - m1) % d) return -1;   // 拼出来的同余方程无解,继续不下去
        k1 *= (m2 - m1) / d;            // 同比扩大指定倍数
        int t = abs(a2 / d);            // 最小累计单元
        k1 = (k1 % t + t) % t;          // 取得最小正整数解,否则可以会在计算过程中爆long long
        m1 = k1 * a1 + m1;              // 准备下一轮的m[1],a[1]
        a1 = abs(a1 / d * a2);          // 两个顺序不能反
    }
    return (m1 % a1 + a1) % a1; // 取得最小正整数解
}

int T;
int m1, m2, m3, d;
signed main() {
    while (cin >> m[1] >> m[2] >> m[3] >> d && ~m[1]) {
        T++;
        int res = exCRT();
        int A = a[1] * a[2] * a[3];
        if (res <= d) res += A;
        printf("Case %lld: the next triple peak occurs in %lld days.\n", T, res - d);
    }
}

P3868 [TJOI2009] 猜数字

Code 中国剩余定理实现代码

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define endl "\n"
const int N = 20;

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

int n;
int a[N];
int m[N];
int A = 1;
int x, y;
int res;

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

void CRT() {
    for (int i = 1; i <= n; i++) cin >> m[i];
    for (int i = 1; i <= n; i++) cin >> a[i];

    for (int i = 1; i <= n; i++) {
        m[i] = (m[i] % a[i] + a[i]) % a[i];
        A *= a[i];
    }

    for (int i = 1; i <= n; i++) {
        int Mi = A / a[i];
        exgcd(Mi, a[i], x, y);
        x = (x % a[i] + a[i]) % a[i];
        res = (res + qmul(m[i], A / a[i] * x, A)) % A;
    }
}

signed main() {
    cin >> n;
    CRT();
    cout << res << endl;
}

Code 扩展中国剩余定理实现代码

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define endl "\n"
const int N = 100010;
int a[N], m[N];
int n;
int exgcd(int a, int b, int &x, int &y) {
    if (!b) {
        x = 1, y = 0;
        return a;
    }
    int d = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return d;
}

int exCRT() { // m代表合并多少次,有n个方程,其实m=n-1
    int a1 = a[1], m1 = m[1], a2, m2, k1, k2;

    for (int i = 2; i <= n; i++) { // a1,m1:做为先遣军,准备与后续的方程合并
        a2 = a[i], m2 = m[i];
        int d = exgcd(a1, -a2, k1, k2); // 扩展欧几里得求方程特解,x有用,y无用
        if ((m2 - m1) % d) return -1;   // 拼出来的同余方程无解,继续不下去
        k1 *= (m2 - m1) / d;            // 同比扩大指定倍数
        int t = abs(a2 / d);            // 最小累计单元
        k1 = (k1 % t + t) % t;          // 取得最小正整数解,否则可以会在计算过程中爆long long
        m1 = k1 * a1 + m1;              // 准备下一轮的m[1],a[1]
        a1 = abs(a1 / d * a2);          // 两个顺序不能反
    }
    return (m1 % a1 + a1) % a1; // 取得最小正整数解
}

signed main() {
    cin >> n;
    for (int i = 1; i <= n; i++) cin >> m[i];
    for (int i = 1; i <= n; i++) cin >> a[i];
    cout << exCRT() << endl;
}

模数不互质

洛谷P4777 【模板】扩展中国剩余定理(EXCRT.

AcWing 204. 表达整数的奇怪方式

上面两道题是同一道题,代码是一样的。

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define endl "\n"

const int N = 100010;
int a[N], m[N];
int n;
int exgcd(int a, int b, int &x, int &y) {
    if (!b) {
        x = 1, y = 0;
        return a;
    }
    int d = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return d;
}

int exCRT() { // m代表合并多少次,有n个方程,其实m=n-1
    int a1 = a[1], m1 = m[1], a2, m2, k1, k2;

    for (int i = 2; i <= n; i++) { // a1,m1:做为先遣军,准备与后续的方程合并
        a2 = a[i], m2 = m[i];
        int d = exgcd(a1, -a2, k1, k2); // 扩展欧几里得求方程特解,x有用,y无用
        if ((m2 - m1) % d) return -1;   // 拼出来的同余方程无解,继续不下去
        k1 *= (m2 - m1) / d;            // 同比扩大指定倍数
        int t = abs(a2 / d);            // 最小累计单元
        k1 = (k1 % t + t) % t;          // 取得最小正整数解,否则可以会在计算过程中爆long long
        m1 = k1 * a1 + m1;              // 准备下一轮的m[1],a[1]
        a1 = abs(a1 / d * a2);          // 两个顺序不能反
    }
    return (m1 % a1 + a1) % a1; // 取得最小正整数解
}
signed main() {
    cin >> n;
    for (int i = 1; i <= n; i++) cin >> a[i] >> m[i];
    cout << exCRT() << endl;
}
posted @   糖豆爸爸  阅读(165)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· 被坑几百块钱后,我竟然真的恢复了删除的微信聊天记录!
· 没有Manus邀请码?试试免邀请码的MGX或者开源的OpenManus吧
· 【自荐】一款简洁、开源的在线白板工具 Drawnix
· 园子的第一款AI主题卫衣上架——"HELLO! HOW CAN I ASSIST YOU TODAY
· Docker 太简单,K8s 太复杂?w7panel 让容器管理更轻松!
历史上的今天:
2022-01-29 AcWing 1055. 股票买卖 II
Live2D
点击右上角即可分享
微信分享提示