中国剩余定理（CRT 孙子定理）——Biorhythms（POJ 1006）

定义：中国古代求解一次同余式组（见同余）的方法。

表述：

设正整数

两两互素，则同余方程组



有整数解。并且在模

下的解是唯一的，解为



其中

，而

为

模

的逆元。

代码：

int CRT(int a[],int m[],int n)
{
int M = 1;
int ans = 0;
for(int i=1; i<=n; i++)
M *= m[i];
for(int i=1; i<=n; i++)
{
int x, y;
int Mi = M / m[i];
extend_Euclid(Mi, m[i], x, y);
int x0 = x*(1/d)%m[i] + m[i];
ans = (ans + Mi * x0 * a[i]) % M;
}
if(ans < 0) ans += M;
return ans;
}

练习：http://poj.org/problem?id=1006

-AC代码：

#include <iostream>
#include <cstring>
#include <cstdio>
#include <cmath>

using namespace std;

int a[4];
int m[4];

int Exgcd(int a, int b, int &x ,int &y)
{
if(b == 0)
{
x = 1, y = 0;
return a;
}
int ret = Exgcd(b, a%b, x, y);
int tmp = x;
x = y;
y = tmp - a/b*y;
return ret;
}

int CRT (int a[], int m [], int n)
{
int M = 1;
int ans = 0;
for(int i=1; i<=n;i++)
M*=m[i];
for(int i=1;i<=n;i++)
{
int x,y;
int Mi = M / m[i];
int d = Exgcd(Mi, m[i], x, y);
int x0 = x*( 1/d)%m[i]+m[i];
ans = (ans + Mi * x0* a[i]) % M;
}
if(ans < 0) ans +=M;
return ans;
}

int main()
{
int p, e, i, d, t = 1;
while(cin>>p>>e>>i>>d)
{
if(p == -1 && e == -1 && i == -1 && d == -1)
break;
a[1] = p;
a[2] = e;
a[3] = i;
m[1] = 23;
m[2] = 28;
m[3] = 33;
int ans = CRT(a, m, 3);
if(ans <= d)
ans += 21252;
cout<<"Case "<<t++<<": the next triple peak occurs in "<<ans - d<<" days."<<endl;
}
return 0;
}

非互质情况下的中国剩余定理（合并求解多个模线性方程组）

（以下叙述及证明大部分转自：http://972169909-qq-com.iteye.com/blog/1266328 作者 KIDx）

给出同余方程组：



其中，bi, ni的值已知，且n1,n2.n3, … , ni两两之间不一定互质，求Res的值。

采用合并方程的做法：





合并方程求同余方程解的CRT模板：

LL CRT2 (LL b[], LL n[], int num) //b[]为余数数组， n[]为模数组
{
int i;
bool flag = false;
LL n1 = n[0], n2, b1 = b[0], b2, bb, d, t, k, x, y;
for (i = 1; i < num; i++)
{
n2 = n[i], b2 = b[i];
bb = b2 - b1;
d = Egcd (n1, n2, x, y);
if (bb % d)     //模线性解k1时发现无解
{
flag = true;
break;
}
k = bb / d * x;    //相当于求上面所说的k1【模线性方程】
t = n2 / d;
if (t < 0) t = -t;
k = (k % t + t) % t;    //相当于求上面的K`
b1 = b1 + n1*k;  //b1相当于最小正整数解
n1 = n1 / d * n2;
}
if (flag)
return 0;           //无解
/******************求正整数解******************/
if (b1 == 0)    //如果解为0，而题目要正整数解，显然不行
b1 = n1;    //n1刚好为所有ni的最小公倍数，就是解了
/******************求正整数解******************/
if (b1 > N)
return 0;
return (N-b1)/n1+1;    //形成的解：b1, b1+n1, b1+2n1,..., b1+xni...
}

练习：HDU 1573 http://acm.hdu.edu.cn/showproblem.php?pid=1573

AC 代码：

#include <iostream>
#include <cstring>
#include <cstdio>
#include <cmath>

using namespace std;
typedef long long ll;

int a[11];
int b[11];
int T[11];
ll N;

ll Exgcd( ll a, ll b, ll &x, ll &y)
{
if(b == 0)
{
x =1, y=0;
return a;
}
else
{
ll ret = Exgcd(b, a%b, x, y);
ll tmp = x;
x = y;
y = tmp -a/b*y;
return ret;
}
}

ll CRT(int b [], int a[], int n)
{
bool flag = false;
ll b1 = b[1], a1 = a[1];
ll x,y;
for(int i=2; i<=n;i++)
{
ll b2 = b[i];
ll a2 = a[i];
ll bb = b2 - b1;
ll d = Exgcd(a1, a2, x, y);
if(bb%d)//模线性解K1 = (b2-b1)/n1 (mode n2/d)时无解
{
flag = true;
break;
}
ll k = bb/d *x;//解上述模线性方程得到K的值
ll t = a2/d;
if(t<0) t = -t;
k = ( k % t + t ) % t; // 求解上述 k`
b1 = b1 + a1 *k;
a1 = a1 /d * a2;
}
if( flag )    //无解情况
return 0;
if( b1 == 0)  //如果解为0，不符合正整数解，换成n1刚好是所有ni的最小公倍数
b1 = a1;
if(b1>N)      //求出所有小于N的正整数解
return 0;
return (N-b1)/a1+1;
}

int main()
{
int t;
cin >> t;
while(t--)
{
int m;
scanf("%lld%d", &N, &m);

for(int i=1;i<=m;i++)
cin >> a[i];
for(int i=1;i<=m;i++)
cin >> b[i];

ll ans = CRT(b, a, m);

cout << ans << endl;
}

return 0;
}