中国剩余定理CRT
中国剩余定理是用来求线性同于方程组的。
[egin{aligned}
left { egin{matrix}
x equiv c_1 (mod \,\,m_1 )\
x equiv c_2 (mod \,\,m_2 )\
...\
x equiv c_n(mod \,\, m_n)
end{matrix}
ight.
end{aligned}
]
中国剩余定理是这样来的。
我们先考虑如下几个方程组：
[egin{aligned}left { egin{matrix} x_1 equiv 1 (mod \,\,m_1 )\ x_1 equiv 0 (mod \,\,m_2 )\...\x_1 equiv 0(mod \,\, m_n)end{matrix}
ight.,left { egin{matrix}x_2 equiv 0 (mod \,\,m_1 )\x_2 equiv 1 (mod \,\,m_2 )\...\x_2 equiv 0(mod \,\, m_n)end{matrix}
ight.,left { egin{matrix}x_n equiv 0 (mod \,\,m_1 )\x_n equiv 0 (mod \,\,m_2 )\...\x_n equiv 1(mod \,\, m_n)end{matrix}
ight.end{aligned}
]
那么，在(gcd(m1,m2,...,m_n))时，就显然有如下结论：
若
[egin{aligned}
prod_{i
eq 1} m_i x_1' & equiv 1(mod \,\, m_1)\
prod_{i
eq 2} m_i x_2' & equiv 1(mod \,\, m_2)\
prod_{i
eq n} m_i x_n' & equiv 1(mod \,\, m_n)\
end{aligned}
]
那么
[egin{aligned}
x_1 & =x_1'prod_{i
eq 1} m_i \
x_2 & =x_2'prod_{i
eq 2} m_i \
x_n & =x_n'prod_{i
eq n} m_i \
end{aligned}
]
然后最终的结果就是
[x=sum_{i=1}^nc_ix_i+kprod_{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
", Ans );
return 0;
}
扩展中国剩余定理ExCRT
刚才提到，中国剩余定理适用于模数互质的时候。要是模数不互质，那么就需要用到扩展中国剩余定理。
ExCRT是这样工作的：
我们先观察两个线性同余方程组：
[egin{aligned}
x & equiv c_1 (mod\,\,m_1)\
x & equiv c_2 (mod\,\,m_2)
end{aligned}
]
我们将它写成这种形式：
[egin{aligned}
x & =c_1+k_1m_1 \
x & =c_2+k_2m_2
end{aligned}
]
联立后得到：
[egin{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))。
这样的话，我们又可以得到：
[egin{aligned}
&k_1frac{m_1}{gcd(m_1,m_2)}-k_2frac{m_2}{gcd(m_1,m_2)}=frac{c_2-c_1}{gcd(m_1,m_2)}\
Rightarrow &k_1frac{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_1equivfrac{c_2-c_1}{gcd(m_1,m_2)} imes(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)中，得到
[xequivfrac{c_2-c_1}{gcd(m_1,m_2)} imes (frac{m_1}{gcd(m_1,m_2)})^{-1} imes m_1+c_2\,\,(mod \,\, frac{m_2}{gcd(m_1,m_2)})
]
我们又得到了一个形如(xequiv 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
" );
return;
}
}
if( T1.B < 0 ) T1.B += T1.A;
if( T1.B > N ) {
printf( "0
" );
return;
}
if( T1.B ) printf( "%d
", ( int ) ( ( N - T1.B ) / T1.A + 1 ) );
else printf( "%d
", ( int ) ( N / T1.A ) );
return;
}