中国剩余定理及扩展中国剩余定理
中国剩余定理CRT
中国剩余定理是用来求线性同于方程组的。
\[\begin{aligned}
\left \{ \begin{matrix}
x \equiv c_1 (mod \,\,m_1 )\\
x \equiv c_2 (mod \,\,m_2 )\\
...\\
x \equiv c_n(mod \,\, m_n)
\end{matrix} \right.
\end{aligned}
\]
中国剩余定理是这样来的。
我们先考虑如下几个方程组:
\[\begin{aligned}\left \{ \begin{matrix} x_1 \equiv 1 (mod \,\,m_1 )\\ x_1 \equiv 0 (mod \,\,m_2 )\\...\\x_1 \equiv 0(mod \,\, m_n)\end{matrix} \right.,\left \{ \begin{matrix}x_2 \equiv 0 (mod \,\,m_1 )\\x_2 \equiv 1 (mod \,\,m_2 )\\...\\x_2 \equiv 0(mod \,\, m_n)\end{matrix} \right.,\left \{ \begin{matrix}x_n \equiv 0 (mod \,\,m_1 )\\x_n \equiv 0 (mod \,\,m_2 )\\...\\x_n \equiv 1(mod \,\, m_n)\end{matrix} \right.\end{aligned}
\]
那么,在\(gcd(m1,m2,...,m_n)\)时,就显然有如下结论:
若
\[\begin{aligned}
\prod_{i \neq 1} m_i x_1' & \equiv 1(mod \,\, m_1)\\
\prod_{i \neq 2} m_i x_2' & \equiv 1(mod \,\, m_2)\\
\prod_{i \neq n} m_i x_n' & \equiv 1(mod \,\, m_n)\\
\end{aligned}
\]
那么
\[\begin{aligned}
x_1 & =x_1'\prod_{i\neq 1} m_i \\
x_2 & =x_2'\prod_{i\neq 2} m_i \\
x_n & =x_n'\prod_{i\neq n} m_i \\
\end{aligned}
\]
然后最终的结果就是
\[x=\sum_{i=1}^nc_ix_i+k\prod_{i=1}^nm_i
\]
那么程序就比较好写啦。
#include <bits/stdc++.h>
#define LL long long
using namespace std;
LL n;
LL A[ 20 ], B[ 20 ];
LL M, Ans, T[ 20 ];
LL QM( LL x, LL y ) {
LL Ans = 0;
for( ; y; y >>= 1, x = x * 2 % M )
if( y & 1 ) Ans = ( Ans + x ) % M;
return Ans;
}
void Expower( LL a, LL b, LL &x, LL &y ) {
if( b == 0 ) {
x = 1; y = 0; return;
}
Expower( b, a % b, y, x );
y -= a / b * x;
return;
}
LL INV( LL a, LL b ) {
LL x, y;
Expower( a, b, x, y );
if( x < 0 ) x += b;
return x;
}
int main() {
scanf( "%lld", &n );
for( LL i = 1; i <= n; ++i ) scanf( "%lld", &A[ i ] );
for( LL i = 1; i <= n; ++i ) scanf( "%lld", &B[ i ] );
for( LL i = 1; i <= n; ++i ) A[ i ] %= B[ i ];
M = 1;
for( LL i = 1; i <= n; ++i ) M *= B[ i ];
for( LL i = 1; i <= n; ++i ) T[ i ] = QM( INV( M / B[ i ], B[ i ] ), ( M / B[ i ] ) );
Ans = 0;
for( LL i = 1; i <= n; ++i ) Ans = ( Ans + QM( A[ i ], T[ i ] ) ) % M;
printf( "%lld\n", Ans );
return 0;
}
扩展中国剩余定理ExCRT
刚才提到,中国剩余定理适用于模数互质的时候。要是模数不互质,那么就需要用到扩展中国剩余定理。
ExCRT是这样工作的:
我们先观察两个线性同余方程组:
\[\begin{aligned}
x & \equiv c_1 (mod\,\,m_1)\\
x & \equiv c_2 (mod\,\,m_2)
\end{aligned}
\]
我们将它写成这种形式:
\[\begin{aligned}
x & =c_1+k_1m_1 \\
x & =c_2+k_2m_2
\end{aligned}
\]
联立后得到:
\[\begin{aligned}
c_1+k_1m_1 &=c_2+k_2m_2\\
\Rightarrow k_1m_1-k_2m_2&= c_2-c_1
\end{aligned}
\]
由裴蜀定理得,方程有解的充要条件是\(gcd(m_1,m_2)|(c_2-c_1)\)。
这样的话,我们又可以得到:
\[\begin{aligned}
&k_1\frac{m_1}{gcd(m_1,m_2)}-k_2\frac{m_2}{gcd(m_1,m_2)}=\frac{c_2-c_1}{gcd(m_1,m_2)}\\
\Rightarrow &k_1\frac{m_1}{gcd(m_1,m_2)}\equiv \frac{c_2-c_1}{gcd(m_1,m_2)}\,\,(mod\,\, \frac{m_2}{gcd(m_1,m_2)})\\
\Rightarrow & k_1\equiv\frac{c_2-c_1}{gcd(m_1,m_2)}\times(\frac{m_1}{gcd(m_1,m_2)})^{-1} \,\,(mod\,\,\frac{m_2}{gcd(m_1,m_2)})
\end{aligned}
\]
然后将\(k_1\)代回\(x=c_1+k_1m_1\)中,得到
\[x\equiv\frac{c_2-c_1}{gcd(m_1,m_2)}\times (\frac{m_1}{gcd(m_1,m_2)})^{-1}\times m_1+c_2\,\,(mod \,\, \frac{m_2}{gcd(m_1,m_2)})
\]
我们又得到了一个形如\(x\equiv c \,\,(mod\,\,m)\)的线性同余方程。所以迭代求解即可。
#include <bits/stdc++.h>
#define LL long long
using namespace std;
const int Maxm = 20;
void Work();
int main() {
int TestCases;
scanf( "%d", &TestCases );
for( ; TestCases; --TestCases ) Work();
return 0;
}
int N, M, A[ Maxm ], B[ Maxm ];
struct equation {
LL A, B;
};
equation T1, T2;
LL GCD( LL x, LL y ) {
LL m = x % y;
while( m ) {
x = y; y = m; m = x % y;
}
return y;
}
void ExGCD( LL a, LL b, LL &x, LL &y ) {
if( b == 0 ) {
x = 1; y = 0; return;
}
ExGCD( b, a % b, y, x );
y -= a / b * x;
return;
}
LL Inv( LL a, LL b ) {
LL x, y;
ExGCD( a, b, x, y );
if( x < 0 ) x += b;
return x;
}
equation ExCRT( equation X, equation Y ) {
LL Gcd = GCD( X.A, Y.A );
if( ( Y.B - X.B ) % Gcd ) return ( equation ) { 0, 0 };
LL A = X.A * Y.A / Gcd;
LL B = Inv( X.A / Gcd, Y.A / Gcd ) * ( Y.B - X.B ) / Gcd % ( Y.A / Gcd ) * X.A + X.B;
return ( equation ) { A, B };
}
void Work() {
scanf( "%d%d", &N, &M );
for( int i = 1; i <= M; ++i ) scanf( "%d", &A[ i ] );
for( int i = 1; i <= M; ++i ) scanf( "%d", &B[ i ] );
T1 = ( equation ) { A[ 1 ], B[ 1 ] };
for( int i = 2; i <= M; ++i ) {
T2 = ( equation ) { A[ i ], B[ i ] };
T1 = ExCRT( T1, T2 );
if( !T1.A ) {
printf( "0\n" );
return;
}
}
if( T1.B < 0 ) T1.B += T1.A;
if( T1.B > N ) {
printf( "0\n" );
return;
}
if( T1.B ) printf( "%d\n", ( int ) ( ( N - T1.B ) / T1.A + 1 ) );
else printf( "%d\n", ( int ) ( N / T1.A ) );
return;
}