线性同余方程和中国剩余定理学习笔记

线性同余方程介绍

形如 a x ≡ c ( m o d b ) ax \equiv c \pmod b ax≡c(modb) 的方程被称为 线性同余方程(Congruence Equation)。

求解方法

根据以下两个定理，我们可以求出同余方程 a x ≡ c ( m o d b ) ax \equiv c \pmod b ax≡c(modb) 的解。

定理 1：方程 a x + b y = c ax+by=c ax+by=c 与方程 a x ≡ c ( m o d b ) ax \equiv c \pmod b ax≡c(modb) 是等价的，有整数解的充要条件为 gcd ⁡ ( a , b ) ∣ c \gcd(a,b) \mid c gcd(a,b)∣c。

根据定理 1，方程 a x + b y = c ax+by=c ax+by=c，我们可以先用扩展欧几里得算法求出一组 x 0 , y 0 x_0,y_0 x0,y0，也就是 a x 0 + b y 0 = gcd ⁡ ( a , b ) ax_0+by_0=\gcd(a,b) ax0+by0=gcd(a,b)，然后两边同时除以 gcd ⁡ ( a , b ) \gcd(a,b) gcd(a,b)，再乘 c c c。然后就得到了方程 a c gcd ⁡ ( a , b ) x 0 + b c gcd ⁡ ( a , b ) y 0 = c a\dfrac{c}{\gcd(a,b)}x_0+b\dfrac{c}{\gcd(a,b)}y_0=c agcd(a,b)cx0+bgcd(a,b)cy0=c，然后我们就找到了方程的一个解。

定理 2：若 gcd ⁡ ( a , b ) = 1 \gcd(a,b)=1 gcd(a,b)=1，且 x 0 x_0 x0、 y 0 y_0 y0 为方程 a x + b y = c ax+by=c ax+by=c 的一组解，则该方程的任意解可表示为： x = x 0 + b t x=x_0+bt x=x0+bt， y = y 0 − a t y=y_0-at y=y0−at, 且对任意整数 t t t 都成立。

根据定理 2，可以求出方程的所有解。但在实际问题中，我们往往被要求求出一个最小整数解，也就是一个特解 x = ( x m o d t + t ) m o d t x=(x \bmod t+t) \bmod t x=(xmodt+t)modt，其中 t = b gcd ⁡ ( a , b ) t=\dfrac{b}{\gcd(a,b)} t=gcd(a,b)b。

中国剩余定理介绍

「物不知数」问题

有物不知其数，三三数之剩二，五五数之剩三，七七数之剩二。问物几何？

即求满足以下条件的整数：除以 3 3 3 余 2 2 2，除以 5 5 5 余 3 3 3，除以 7 7 7 余 2 2 2。

该问题最早见于《孙子算经》中，并有该问题的具体解法。宋朝数学家秦九韶于 1247 年《数书九章》卷一、二《大衍类》对「物不知数」问题做出了完整系统的解答。上面具体问题的解答口诀由明朝数学家程大位在《算法统宗》中给出：

三人同行七十希，五树梅花廿一支，七子团圆正半月，除百零五便得知。

2 × 70 + 3 × 21 + 2 × 15 = 233 = 2 × 105 + 23 2\times 70+3\times 21+2\times 15=233=2\times 105+23 2×70+3×21+2×15=233=2×105+23，故答案为 23 23 23。

算法简介及过程

中国剩余定理 (Chinese Remainder Theorem, CRT) 可求解如下形式的一元线性同余方程组（其中 p 1 , p 2 , ⋯ , p n p_1, p_2, \cdots, p_n p1,p2,⋯,pn 两两互质）：

{ x ≡ a 1 ( m o d p 1 ) x ≡ a 2 ( m o d p 2 ) ⋮ x ≡ a k ( m o d p k ) \begin{cases} x &\equiv a_1 \pmod {p_1} \\ x &\equiv a_2 \pmod {p_2} \\ &\vdots \\ x &\equiv a_k \pmod {p_k} \\ \end{cases} ⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧xxx≡a1(modp1)≡a2(modp2)⋮≡ak(modpk)

上面的「物不知数」问题就是一元线性同余方程组的一个实例。

算法流程

计算所有模数的积 c n t cnt cnt；
对于第 i i i 个方程：
1. 计算 m i = c n t p i m_i=\frac{cnt}{p_i} mi=picnt；
2. 计算 m i m_i mi 在模 n i n_i ni 意义下的逆元 m i − 1 m_i^{-1} mi−1；
3. 计算 c i = m i m i − 1 c_i=m_im_i^{-1} ci=mimi−1（不要对 a i a_i ai 取模）。
方程组的唯一解为： x = ∑ i = 1 k a i c i ( m o d c n t ) x=\sum_{i=1}^k a_ic_i \pmod {cnt} x=∑i=1kaici(modcnt)。

ll extend_gcd(ll a, ll b, ll& x, ll& y) {if (!b) {x = 1; y = 0;return a;}ll d = extend_gcd(b, a % b, y, x);y -= (a / b) * x;return d;
}
ll CRT(int n, vector<ll> p, vector<ll> a) {ll cnt = 1, ans = 0;for (int i = 1; i <= n; i++) cnt = cnt * p[i];for (int i = 1; i <= n; i++) {ll m = cnt / p[i], x, y;extend_gcd(m, p[i], x, y);ans = (ans + a[i] * m * x % cnt) % cnt;}return (ans % cnt + cnt) % cnt;
}

算法的证明

我们需要证明上面算法计算所得的 x x x 对于任意 i = 1 , 2 , ⋯ , n i=1,2,\cdots,n i=1,2,⋯,n 满足 x ≡ a i ( m o d p i ) x\equiv a_i \pmod {p_i} x≡ai(modpi)。

当 i ≠ j i\neq j i=j 时，有 m j ≡ 0 ( m o d p i ) m_j \equiv 0 \pmod {p_i} mj≡0(modpi)，故 c j ≡ m j ≡ 0 ( m o d p i ) c_j \equiv m_j \equiv 0 \pmod {p_i} cj≡mj≡0(modpi)。又有 c i ≡ m i ⋅ ( m i − 1 m o d p i ) ≡ 1 ( m o d p i ) c_i \equiv m_i \cdot (m_i^{-1} \bmod {p_i}) \equiv 1 \pmod {p_i} ci≡mi⋅(mi−1modpi)≡1(modpi)，所以我们有：

x ≡ ∑ j = 1 n a j c j ( m o d p i ) ≡ a i c i ( m o d p i ) ≡ a i ⋅ m i ⋅ ( m i − 1 m o d p i ) ( m o d p i ) ≡ a i ( m o d p i ) \begin{aligned} x&\equiv \sum_{j=1}^n a_jc_j &\pmod {p_i} \\ &\equiv a_ic_i &\pmod {p_i} \\ &\equiv a_i \cdot m_i \cdot (m^{-1}_i \bmod p_i) &\pmod {p_i} \\ &\equiv a_i &\pmod {p_i} \end{aligned} x≡j=1∑najcj≡aici≡ai⋅mi⋅(mi−1modpi)≡ai(modpi)(modpi)(modpi)(modpi)

即对于任意 i = 1 , 2 , ⋯ , n i=1,2,\cdots,n i=1,2,⋯,n，上面算法得到的 x x x 总是满足 x ≡ a i ( m o d p i ) x\equiv a_i \pmod{p_i} x≡ai(modpi)，即证明了解同余方程组的算法的正确性。

因为我们没有对输入的 p i p_i pi 作特殊限制，所以任何一组输入 { a i } \{a_i\} {ai} 都对应一个解 x x x。

另外，若 x ≠ y x\neq y x=y，则总存在 i i i 使得 x x x 和 y y y 在模 a i a_i ai 下不同余。

故系数列表 { a i } \{a_i\} {ai} 与解 x x x 之间是一一映射关系，方程组总是有唯一解。

Garner 算法

CRT 的另一个用途是用一组比较小的质数表示一个大的整数。

例如，若 a a a 满足如下线性方程组，且 a < ∏ i = 1 k p i a < \prod_{i=1}^k p_i a<∏i=1kpi（其中 p i p_i pi 为质数）：

{ a ≡ a 1 ( m o d p 1 ) a ≡ a 2 ( m o d p 2 ) ⋮ a ≡ a k ( m o d p k ) \begin{cases} a &\equiv a_1 \pmod {p_1} \\ a &\equiv a_2 \pmod {p_2} \\ &\vdots \\ a &\equiv a_k \pmod {p_k} \\ \end{cases} ⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧aaa≡a1(modp1)≡a2(modp2)⋮≡ak(modpk)

我们可以用以下形式的式子（称作 a a a 的混合基数表示）表示 a a a：

a = x 1 + x 2 p 1 + x 3 p 1 p 2 + … + x k p 1 … p k − 1 a = x_1 + x_2 p_1 + x_3 p_1 p_2 + \ldots + x_k p_1 \ldots p_{k-1} a=x1+x2p1+x3p1p2+…+xkp1…pk−1

Garner 算法 将用来计算系数 x 1 , … , x k x_1, \ldots, x_k x1,…,xk。

令 r i j r_{ij} rij 为 p i p_i pi 在模 p j p_j pj 意义下的逆：

p i ⋅ r i , j ≡ 1 ( m o d p j ) p_i \cdot r_{i,j} \equiv 1 \pmod{p_j} pi⋅ri,j≡1(modpj)

把 a a a 代入我们得到的第一个方程：

a 1 ≡ x 1 ( m o d p 1 ) a_1 \equiv x_1 \pmod{p_1} a1≡x1(modp1)

代入第二个方程得出：

a 2 ≡ x 1 + x 2 p 1 ( m o d p 2 ) a_2 \equiv x_1 + x_2 p_1 \pmod{p_2} a2≡x1+x2p1(modp2)

方程两边减 x 1 x_1 x1，除 p 1 p_1 p1 后得

a 2 − x 1 ≡ x 2 p 1 ( m o d p 2 ) ( a 2 − x 1 ) r 1 , 2 ≡ x 2 ( m o d p 2 ) x 2 ≡ ( a 2 − x 1 ) r 1 , 2 ( m o d p 2 ) \begin{array}{rclr} a_2 - x_1 &\equiv& x_2 p_1 &\pmod{p_2} \\ (a_2 - x_1) r_{1,2} &\equiv& x_2 &\pmod{p_2} \\ x_2 &\equiv& (a_2 - x_1) r_{1,2} &\pmod{p_2} \end{array} a2−x1(a2−x1)r1,2x2≡≡≡x2p1x2(a2−x1)r1,2(modp2)(modp2)(modp2)

类似地，我们可以得到：

x k = ( . . . ( ( a k − x 1 ) r 1 , k − x 2 ) r 2 , k ) − . . . ) r k − 1 , k m o d p k x_k=(...((a_k-x_1)r_{1,k}-x_2)r_{2,k})-...)r_{k-1,k} \bmod p_k xk=(...((ak−x1)r1,k−x2)r2,k)−...)rk−1,kmodpk

for (int i = 0; i < k; ++i) {x[i] = a[i];for (int j = 0; j < i; ++j) {x[i] = r[j][i] * (x[i] - x[j]);x[i] = x[i] % p[i];if (x[i] < 0) x[i] += p[i];}
}

该算法的时间复杂度为 O ( k 2 ) O(k^2) O(k2)。

应用

某些计数问题或数论问题出于加长代码、增加难度、或者是一些其他原因，给出的模数：不是质数！

但是对其质因数分解会发现它没有平方因子，也就是该模数是由一些不重复的质数相乘得到。

那么我们可以分别对这些模数进行计算，最后用 CRT 合并答案。

下面这道题就是一个不错的例子。

“洛谷 P2480 [SDOI2010]古代猪文”
给出 G , n G,n G,n（ 1 ≤ G , n ≤ 1 0 9 1 \leq G,n \leq 10^9 1≤G,n≤109），求：
G ∑ k ∣ n ( n k ) m o d 999 911 659 G^{\sum_{k\mid n}\binom{n}{k}} \bmod 999~911~659 G∑k∣n(kn)mod999 911 659
首先，当 G = 999 911 659 G=999~911~659 G=999 911 659 时，所求显然为 0 0 0。

否则，根据欧拉定理，可知所求为：

G ∑ k ∣ n ( n k ) m o d 999 911 658 m o d 999 911 659 G^{\sum_{k\mid n}\binom{n}{k} \bmod 999~911~658} \bmod 999~911~659 G∑k∣n(kn)mod999 911 658mod999 911 659

现在考虑如何计算：

∑ k ∣ n ( n k ) m o d 999 911 658 \sum_{k\mid n}\binom{n}{k} \bmod 999~911~658 k∣n∑(kn)mod999 911 658

因为 999 911 658 999~911~658 999 911 658 不是质数，无法保证 ∀ x ∈ [ 1 , 999 911 657 ] \forall x \in [1,999~911~657] ∀x∈[1,999 911 657]， x x x 都有逆元存在，上面这个式子我们无法直接计算。

注意到 999 911 658 = 2 × 3 × 4679 × 35617 999~911~658=2 \times 3 \times 4679 \times 35617 999 911 658=2×3×4679×35617，其中每个质因子的最高次数均为一，我们可以考虑分别求出 ∑ k ∣ n ( n k ) \sum_{k\mid n}\binom{n}{k} ∑k∣n(kn) 在模 2 2 2， 3 3 3， 4679 4679 4679， 35617 35617 35617 这几个质数下的结果，最后用中国剩余定理来合并答案。

也就是说，我们实际上要求下面一个线性方程组的解：

{ x ≡ a 1 ( m o d 2 ) x ≡ a 2 ( m o d 3 ) x ≡ a 3 ( m o d 4679 ) x ≡ a 4 ( m o d 35617 ) \begin{cases} x \equiv a_1 \pmod 2\\ x \equiv a_2 \pmod 3\\ x \equiv a_3 \pmod {4679}\\ x \equiv a_4 \pmod {35617} \end{cases} ⎩⎪⎪⎪⎨⎪⎪⎪⎧x≡a1(mod2)x≡a2(mod3)x≡a3(mod4679)x≡a4(mod35617)

而计算一个组合数对较小的质数取模后的结果，可以利用卢卡斯定理。

扩展：模数不互质的情况

两个方程

设两个方程分别是 x ≡ a 1 ( m o d m 1 ) x\equiv a_1 \pmod {m_1} x≡a1(modm1)、 x ≡ a 2 ( m o d m 2 ) x\equiv a_2 \pmod {m_2} x≡a2(modm2)；

将它们转化为不定方程： x = m 1 p + a 1 = m 2 q + a 2 x=m_1p+a_1=m_2q+a_2 x=m1p+a1=m2q+a2，其中 p , q p, q p,q 是整数，则有 m 1 p − m 2 q = a 2 − a 1 m_1p-m_2q=a_2-a_1 m1p−m2q=a2−a1。

由裴蜀定理，当 a 2 − a 1 a_2-a_1 a2−a1 不能被 gcd ⁡ ( m 1 , m 2 ) \gcd(m_1,m_2) gcd(m1,m2) 整除时，无解；

其他情况下，可以通过扩展欧几里得算法解出来一组可行解 ( p , q ) (p, q) (p,q)；

则原来的两方程组成的模方程组的解为 x ≡ b ( m o d M ) x\equiv b\pmod M x≡b(modM)，其中 b = m 1 p + a 1 b=m_1p+a_1 b=m1p+a1， M = lcm ( m 1 , m 2 ) M=\text{lcm}(m_1, m_2) M=lcm(m1,m2)。

多个方程

用上面的方法两两合并即可。