从中国剩余定理谈POJ 1006

POJ 1006 – Biorhythms

试题链接：POJ 1006

一、费马小定理

若$p$为素数，且$a,p$互素，则
$$a^{p-1}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$$

二、欧拉定理

若$n,a$为整数，且$n,a$互素，则
$$a^{\varphi (n)}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}n)$$
其中，$\varphi (n)$为欧拉函数，表示小于$n$的正整数中与$n$互质的数的数目。欧拉函数的通式为
$$\varphi (n)=n\prod _{i=1}^k\left(1-\frac{1}{p_i}\right)$$
其中，$p_1,p_2,\cdots ,p_k$为$n$的所有素因数。

三、逆元

1、定义

若
$$ax\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}m)$$
则称$x=a^{-1}$为$a$在模$m$意义下的逆元（数论倒数）。

2、求法

求解数论逆元的方法有很多，下面先介绍根据费马小定理来推出计算公式。

（1）根据费马小定理（当模为素数）

根据费马小定理
$$a^{p-1}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$$
则有
$$a{a}^{p-2}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$$
所以，$a$在模$p$（$p$为素数）意义下的逆元
$$a^{-1}\equiv a^{p-2}(\mathrm{m}\mathrm{o}\mathrm{d}p)$$
因而，我们可以使用快速幂取模来计算逆元。

（2）根据欧拉定理（对于任意模）

根据欧拉定理，不难推出
$$a^{\varphi (n)-1}a\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}n)$$
所以，$a$在模$n$意义下的逆元
$$a^{-1}\equiv a^{\varphi (n)-1}(\mathrm{m}\mathrm{o}\mathrm{d}n)$$
同样可以使用快速幂取模求得。

四、快速幂取模

1、原理

对于一个十进制正整数$n=\overline{b_{k-1}{b}_{k-2}\cdots b_1{b}_0}$可以表示为$k$位二进制数，即
$$n=b_{k-1}{2}^{k-1}+b_{k-2}{2}^{k-2}+\cdots +b_1{2}^1+b_0{2}^0$$
其中， b i ∈ { 0 , 1 } ( i = 0 , 1 , ⋯   , k − 1 ) b_i \in \{0, 1\} (i = 0, 1, \cdots, k-1) bi​∈{0,1}(i=0,1,⋯,k−1)。

那么对于指数
$$\begin{array}{rl}{x}^n& =x^{b_{k-1}{2}^{k-1}+b_{k-2}{2}^{k-2}+\cdots +b_1{2}^1+b_0{2}^0}\\ & =x^{b_{k-1}{2}^{k-1}}\cdot x^{b_{k-2}{2}^{k-2}}\cdots x^{b_1{2}^1}\cdot x^{b_0{2}^0}\end{array}$$

其中，有一些$b_i$为0，即对应的指数项为1，不会影响最后的结果。如
$$x^{11}=x^{{\overline{1011}}_2}=x^8\cdot x^2\cdot x^1$$
$x^4$没有出现在最后的乘积中，因为$b_2=0$。

而
$$x^{2^i}=x^{2^{i-1}\cdot 2}=x^{2^{i-1}+2^{i-1}}=x^{2^{i-1}}\cdot x^{2^{i-1}}={\left(x^{2^{i-1}}\right)}^2$$

所以，就可以使用迭代来计算指数，时间复杂度为$O(n)$。

2、伪代码

// calculate (x^n) % M
int exp(int x, int n, int M)
{
    int res = 1;            // the answer
    
    while (n)               // for each bit of n, start from b_0
    {
        if (n & 1)          // if b_i = 1
            res = (res * x) % M;    // multiply into residue
        x = (x * x) % M;    // x^{2^i} = (x^{2^{i-1}})^2
        n >>= 1;            // move to next bit b_{i+1}
    }
    
    return res;
}

五、中国剩余定理（孙子定理）

假设整数$m_1,m_2,\cdots ,m_n$两两互素，则对于任意的整数$a_1,a_2,\cdots ,a_n$，同余方程组
$$\{\begin{array}{rl}x& \equiv a_1(\mathrm{m}\mathrm{o}\mathrm{d}{m}_1)\\ x& \equiv a_2(\mathrm{m}\mathrm{o}\mathrm{d}{m}_2)\\ & ⋮\\ x& \equiv a_n(\mathrm{m}\mathrm{o}\mathrm{d}{m}_n)\end{array}$$
有解。

设
$$M=m_1\times m_2\times \cdots \times m_n=\prod _{i=1}^n{m}_i$$
是整数$m_1,m_2,\cdots ,m_n$的乘积，并设
$$M_i=\frac{M}{m_i}$$
是除$m_i$以为其他$n-1$个模的乘积，其中$i=1,2,\cdots ,n$。且记$M_i^{-1}$为$M_i$模$m_i$的逆元，即
$$M_i{M}_i^{-1}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}{m}_i)$$
则解为
$$x\equiv a_1{M}_1^{-1}{M}_1+a_2{M}_2^{-1}{M}_2+\cdots +a_n{M}_n^{-1}{M}_n(\mathrm{m}\mathrm{o}\mathrm{d}M)$$

六、本题解答

1、推导

根据题意，设$x$为从当前年第一天开始到三个高峰同时到来之间的天数，我们可以得到同余方程组
$$\{\begin{array}{rl}x& \equiv p(\mathrm{m}\mathrm{o}\mathrm{d}23)\\ x& \equiv e(\mathrm{m}\mathrm{o}\mathrm{d}28)\\ x& \equiv i(\mathrm{m}\mathrm{o}\mathrm{d}33)\end{array}$$
于是题目要求计算的从当前年给定天$d$到三个高峰同时到来之间的天数$s$为
$$s=(x-d)\text{ mod }21252$$

由于上述方程组形式已定，所以可以通过中国剩余定理人工求出解的公式，然后将特定的$p,e,i$值代入即可。

本题中，$M=23\times 28\times 33=21252$， $M_1=28\times 33=924$， $M_2=23\times 33=759$， $M_3=23\times 28=644$。

由于$m_1=23$为素数，通过快速幂取模求出$M_1$的逆元，即
$$M_1^{-1}\equiv M_1^{m_1-2}(\mathrm{m}\mathrm{o}\mathrm{d}{m}_1)$$
得到$M_1^{-1}\equiv 6(\mathrm{m}\mathrm{o}\mathrm{d}23)$。

又$m_2=28=2^2\times 7$，其素因数为2和7，则$\varphi (m_2)=28\times \left(1-\frac{1}{2}\right)\times \left(1-\frac{1}{7}\right)=12$，于是$M_2^{-1}\equiv 19(\mathrm{m}\mathrm{o}\mathrm{d}28)$。

同理$m_3=33=3\times 11$，其素因数为3和11，则$\varphi (m_3)=33\times \left(1-\frac{1}{3}\right)\times \left(1-\frac{1}{11}\right)=20$，于是$M_3^{-1}\equiv 2(\mathrm{m}\mathrm{o}\mathrm{d}33)$。

所以，方程组的解为
$$\begin{array}{rl}x& \equiv p\times M_1\times M_1^{-1}+e\times M_2\times M_2^{-1}+m\times M_3\times M_3^{-1}(\mathrm{m}\mathrm{o}\mathrm{d}M)\\ & \equiv 5544p+14421e+1288i(\mathrm{m}\mathrm{o}\mathrm{d}21252)\end{array}$$

2、代码

#include <stdio.h>

int main()
{
    int p, e, i, d, x, s;
    int cnt = 0;

    while (scanf("%d %d %d %d", &p, &e, &i, &d) && p >= 0)
    {
        ++cnt;
        x = (5544 * p + 14421 * e + 1288 * i) % 21252;
        s = (21252 + x - d) % 21252;
        if (s == 0)
            s = 21252;
        printf("Case %d: the next triple peak occurs in %d days.\n", cnt, s);
    }

    return 0;
}