中国剩余定理
问题引入
在《孙子算经》中有这样一个问题:“今有物不知其数,三三数之剩二(除以3余2),五五数之剩三(除以5余3),七七数之剩二(除以7余2),问物几何?”这个问题称为“孙子问题”,该问题的一般解法国际上称为“中国剩余定理”。
我们把上述问题抽象出来,就是求解同余方程组:
\begin{cases} x\equiv b_1(\operatorname{mod}a_1) \\ x\equiv b_2(\operatorname{mod}a_2) \\ x\equiv b_3(\operatorname{mod}a_3) \\ \cdots \\ x\equiv b_n(\operatorname{mod}a_n) \end{cases}
其中a_1,a_2,\cdots,a_n互质。
解法
设M=\prod_{i=1}^na_i,\frac{M}{a_i}|n_i且n_i\equiv b_i(\operatorname{mod}a_i).
那么方程的一个解为:
\sum_{i=1}^nn_i
虽然解法难以想到,但原理却容易验证;我们设x=\sum_{i=1}^nn_i,则:
x\operatorname{mod}a_i=\sum_{j=1}^n(n_j\operatorname{mod}a_i)\operatorname{mod}a_i
由于对k\not=i都有n_k=\frac{M}{a_k}=a_1a_2\cdots a_{k-1}a_{k+1}\cdots a_n,其包含因子a_i,则n_j\operatorname{mod}a_i=0.
故:
x\operatorname{mod}a_i=\sum_{j=1}^n(n_j\operatorname{mod}a_i)\operatorname{mod}a_i=n_i\operatorname{mod}a_i=b_i
那么我们如何寻找n_i呢?我们设n_i=t\frac{M}{a_i},则:
t\frac{M}{a_i}\equiv b_i(\operatorname{mod}a_i)
由欧拉定理:
\left(\frac{M}{a_i}\right)^{\phi(a_i)}\equiv1(\operatorname{mod}a_i)
故:
t=\left(\frac{M}{a_i}\right)^{\phi(a_i)-1}
即:
n_i=\left(\frac{M}{a_i}\right)^{\phi(a_i)}
而x的一个解为:
\sum_{i=1}^nn_i
其最小正整数解为:
\sum_{i=1}^nn_i\operatorname{mod} M
这是因为x-a_it\equiv b_i(\operatorname{mod}a_i).
例题
https://www.luogu.com.cn/problem/P1495
#include<cstdio>
using namespace std;
int n,cnt;
int prime[500001];
long long phi[2000001];
bool vis[2000001];
long long m;
long long a[11];
long long b[11];
long long multi(long long a,long long b,long long p)
{
return ((a * b - (long long)(((long double)a * b + 0.5) / p) * p) % p + p) % p;
}
long long quickpower(long long a,long long b,long long mod)
{
long long t = 1;
while (b)
{
if (b & 1)
t = multi(t,a,mod);
a = multi(a,a,mod);
b >>= 1;
}
return t;
}
int main()
{
phi[1] = 1;
for (int i = 2;i <= 2000000;i ++)
{
if (!vis[i])
{
prime[++cnt] = i;
phi[i] = i - 1;
}
for (int j = 1;j <= cnt && prime[j] * i <= 2000000;j ++)
{
vis[i * prime[j]] = true;
if (i % prime[j] == 0)
{
phi[i * prime[j]] = phi[i] * prime[j];
break;
}
phi[i * prime[j]] = phi[i] * phi[prime[j]];
}
}
scanf("%d",&n);
for (int i = 1;i <= n;i ++)
scanf("%lld%lld",&a[i],&b[i]);
m = 1;
for (int i = 1;i <= n;i ++)
m *= a[i];
long long ans = 0;
for (int i = 1;i <= n;i ++)
{
long long now = m / a[i];
ans += multi(multi(b[i],now,m),quickpower(now % a[i],phi[a[i]] - 1,a[i]),m);
ans %= m;
}
printf("%lld",ans);
return 0;
}
