22 KiB

Raw Permalink Blame History Unescape Escape

This file contains ambiguous Unicode characters that may be confused with others in your current locale. If your use case is intentional and legitimate, you can safely ignore this warning. Use the Escape button to highlight these characters.

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

一、问题引入

先看一道poj上的题目：【POJ1006】 Biorhythms

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

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

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


\large \left\{\begin{matrix}
S=N_1+T_1∗k_1 & ① \\ 
S=N_2+T_2∗k_2 & ② \\ 
S=N_3+T_3∗k_3 & ③ \\ 
\end{matrix}\right.

N_1,N_2,N_3分别为为体力，情感，智力出现峰值的日期， T_1,T_2,T_3 分别为体力，情感，智力周期。需要求出k_1,k_2,k_3三个非负整数使上面的等式成立。

想直接求出k_1,k_2,k_3很难，但是目的是求出S，可以考虑 从结果逆推。根据上面的等式，S 满足三个要求：

除以T_1余数为N_1
除以T_2余数为N_2
除以T_3余数为N_3

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

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

二、中国剩余定理

在《孙子算经》中有这样一个问题：“今有物不知其数，三三数之剩二（除以3余2），五五数之剩三（除以5余3），七七数之剩二（除以7余2），问物几何？”这个问题称为 孙子问题 ，该问题的一般解法国际上称为 中国剩余定理。

解法共三步

① 找三个数：
- 从3和5的公倍数中找出被7除余1的最小数15
- 从3和7的公倍数中找出被5除余1的最小数21
- 从5和7的公倍数中找出被3除余1的最小数70
② 找符合条件的一个数：
- 用15乘以2（2为最终结果除以7的余数）
- 用21乘以3（3为最终结果除以5的余数）
- 用70乘以2（2为最终结果除以3的余数）然后把三个乘积相加15∗2+21∗3+70∗2=233
③ 用233除以3，5，7三个数的最小公倍数105，得到余数，即233\%105=23

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

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

首先，我们假设:

n_1是满足除以3余2的一个数，比如2，5，8等等，也就是满足3∗k_1+2（k_1>=0） 的一个任意数。
n_2是满足除以5余3的一个数，比如3，8，13等等，也就是满足5∗k_2+3（k_2>=0） 的一个任意数。
n_3是满足除以7余2的一个数，比如2，9，16等等，也就是满足7∗k_3+2（k_3>=0） 的一个任意数。

有了前面的假设，我们先从n_1这个角度出发，已知n_1满足除以3余2，能不能使得n_1+n_2的和仍然满足除以3余2,进而使得n_1+n_2+n_3的和仍然满足除以3余2？

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

以此定理为依据，如果n_2是3的倍数，n_1+n_2就依然满足除以3余2。同理，如果n_3也是3的倍数，那么n_1+n_2+n_3的和就满足除以3余2。这是从n_1的角度考虑的，再从n_2，n_3的角度出发，我们可推导出以下三点：

为使n_1+n_2+n_3的和满足除以3余2，n_2和n_3必须是3的倍数
为使n_1+n_2+n_3的和满足除以5余3，n_1和n_3必须是5的倍数
为使n_1+n_2+n_3的和满足除以7余2，n_1和n_2必须是7的倍数

因此，为使n_1+n_2+n_3的和作为 孙子问题 的一个最终解，需满足：

n_1除以3余2，且是5和7的公倍数
n_2除以5余3，且是3和7的公倍数
n_3除以7余2，且是3和5的公倍数

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

从5和7的公倍数中找一个除以3余2的数n_1
从3和7的公倍数中找一个除以5余3的数n_2
从3和5的公倍数中找一个除以7余2的数n_3

再将三个数相加得到解。在求n_1，n_2，n_3时又用了一个 小技巧，以n_1为例，并非从5和7的公倍数中直接找一个除以3余2的数，而是 先找一个除以3余1的数，再乘以2。也就是先求出5和7的公倍数模3下的逆元，再用逆元去乘余数!

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

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

最后，我们还要清楚一点，n_1+n_2+n_3只是问题的一个解，并不是最小的解。如何得到最小解？我们只需要从中最大限度的减掉3，5，7的公倍数105即可。道理就是前面讲过的定理 如果a\%b=c，则有(a−k∗b)\%b=c ,所以(n_1+n_2+n_3)\%105就是最终的最小解。

【中国剩余定理】

设x \in Z,有


\large \left\{\begin{matrix}
 x \equiv m_1(mod \ a_1) & \\ 
 x \equiv m_2(mod \ a_2) & \\
 x \equiv m_3(mod \ a_3) & \\
 ... &  \\
 x \equiv m_n(mod \ a_n) &
\end{matrix}\right.

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

则

\large res=\sum_{i=1}^{n}m_iA_it_i

【扩展欧几里得求逆元】

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

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

x 是a模b的乘法逆元，即a * x \equiv 1(mod \ b)
y 是b模a的乘法逆元，即b * x \equiv 1(mod \ a)

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

【`r = (r % a + a) % a`】

作用：将余数加上若干个循环节，使得结果落在0 \sim a-1之间。

举栗子：按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:中国剩余定理要求所有的m_i互质,那么如果不互质，怎么求解同余方程组? 答：对于模数不互质的情况，我们可以通过合并相邻两组方程的方式来逐步得到解。


\large \left\{\begin{matrix}
 x \equiv m_1 \ (mod \ a_1) & ① \\ 
 x \equiv m_2 \ (mod \ a_2) & ② \\ 
 ... \\
 x \equiv m_n \ (mod \ a_n) &  \\ 
\end{matrix}\right.

先计算两个线性同余方程组的解，之后依此类推 ① 可以写成x=k_1*a_1+m_1 ③ ② 可以写成x=k_2*a_2+m_2 ④

③=④ \rightarrow k_1*a_1+m_1=k_2*a_2+m_2 \rightarrow k_1*a_1-k_2*a_2=m_2-m_1 ⑤

根据 裴蜀定理,当 m_2-m_1为 gcd(a_1,-a_2)的倍数时，方程 ⑤ 有无穷多组整数解

而k_1*a_1-k_2*a_2=gcd(a_1,-a_2)=d，可以用 拓展欧几里得算法 来求解出一组特解(k_1',k_2')，即exgcd(a_1,-a_2,k_1',k_2')

如何求解k_1*a_1-k_2*a_2=gcd(a_1,-a_2)=d中的通解k_1,k_2呢？

结论:k_1,k_2的通解


\large \left\{\begin{matrix}
 k_1=k_1'+k*\frac{a_2}{d} \\ 
 k_2=k_2'+k*\frac{a_1}{d} 
\end{matrix}\right.

证明: k_1的通解为s_1,k_2的通解为s_2,k_1的特解为k_1',k_2的特解为k_2'

则


\large \left\{\begin{matrix}
a_1*s_1-a_2*s_2=d & \\ 
a_1*k_1'-a_2*k_2'=d & 
\end{matrix}\right.

↓方程组同时除d


\large \left\{\begin{matrix}
\frac{a_1}{d}*s_1-\frac{a_2}{d}*s_2=1 &  \\ 
\frac{a_1}{d}*k_1'-\frac{a_2}{d}*k_2'=1 & 
\end{matrix}\right.

↓相减下

\frac{a_1}{d}*(s_1-k_1')-\frac{a_2}{d}*(s_2-k_2')=0

\because gcd(a_1,-a_2)=d，\therefore \frac{a_1}{d}与\frac{a_2}{d}互质。

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

\therefore \frac{a_2}{d} \mid s_1-k_1', \frac{a_1}{d} \mid s_2-k_2'

\therefore s_1-k_1'=k*\frac{a_2}{d} \rightarrow k_1=k_1'+k*\frac{a_2}{d}整除式变形

\therefore s_2-k_2'=k*\frac{a_1}{d} \rightarrow k_2=k_2'+k*\frac{a_1}{d}整除式变形

将k_1通解代入 ③ 得

 x=k_1*a_1+m_1=(k_1'+k*\frac{a_2}{d})*a_1+m_1= \underbrace{k_1'*a_1+m_1}_{更新后的m_1}+ \overbrace{\frac{a_2}{d}*a_1}^{更新后的a_1}*k

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


\large \left\{\begin{matrix}
m_1=k_1'*a_1+m_1 & \\ 
a_1=\frac{a_2}{d}*a_1*k & 
\end{matrix}\right.

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

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

四、练习题

模数互质

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;
}

【POJ1006】 Biorhythms

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;
}

22 KiB Raw Permalink Blame History Unescape Escape

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