中国剩余定理

一、中国剩余定理（互质版）

核心代码：

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

【例题】Biorhythms POJ - 1006

Some people believe that there are three cycles in a person’s life that start the day he or she is born. These three cycles are the physical, emotional, and intellectual cycles, and they have periods of lengths 23, 28, and 33 days, respectively. There is one peak in each period of a cycle. At the peak of a cycle, a person performs at his or her best in the corresponding field (physical, emotional or mental). For example, if it is the mental curve, thought processes will be sharper and concentration will be easier.
Since the three cycles have different periods, the peaks of the three cycles generally occur at different times. We would like to determine when a triple peak occurs (the peaks of all three cycles occur in the same day) for any person. For each cycle, you will be given the number of days from the beginning of the current year at which one of its peaks (not necessarily the first) occurs. You will also be given a date expressed as the number of days from the beginning of the current year. You task is to determine the number of days from the given date to the next triple peak. The given date is not counted. For example, if the given date is 10 and the next triple peak occurs on day 12, the answer is 2, not 3. If a triple peak occurs on the given date, you should give the number of days to the next occurrence of a triple peak.
Input
You will be given a number of cases. The input for each case consists of one line of four integers p, e, i, and d. The values p, e, and i are the number of days from the beginning of the current year at which the physical, emotional, and intellectual cycles peak, respectively. The value d is the given date and may be smaller than any of p, e, or i. All values are non-negative and at most 365, and you may assume that a triple peak will occur within 21252 days of the given date. The end of input is indicated by a line in which p = e = i = d = -1.
Output
For each test case, print the case number followed by a message indicating the number of days to the next triple peak, in the form:

Case 1: the next triple peak occurs in 1234 days.

Use the plural form ``days’’ even if the answer is 1.
Sample Input
0 0 0 0
0 0 0 100
5 20 34 325
4 5 6 7
283 102 23 320
203 301 203 40
-1 -1 -1 -1
Sample Output
Case 1: the next triple peak occurs in 21252 days.
Case 2: the next triple peak occurs in 21152 days.
Case 3: the next triple peak occurs in 19575 days.
Case 4: the next triple peak occurs in 16994 days.
Case 5: the next triple peak occurs in 8910 days.
Case 6: the next triple peak occurs in 10789 days.

解法一（暴力枚举）：
AC代码：

#include<iostream>

using namespace std;

int main()
{
    int p,e,i,d,kase=0;
    while(cin>>p>>e>>i>>d)
    {
        if(p==-1)break;
        while(p>=d){p-=23;}
        while(e>=d){e-=28;}
        while(i>=d){i-=33;}
        for(int j=0;;j++)
        {
            if(j>d && (j-p)%23==0 && (j-e)%28==0 && (j-i)%33==0)
            {
                cout<<"Case "<<++kase<<": the next triple peak occurs in "<<j-d<<" days."<<endl;
                break;
            }
        }
    }
    return 0;
}

解法二（中国剩余定理）：
AC代码：

#include<iostream>

using namespace std;

int a[5],m[5];

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

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

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

二、中国剩余定理（非互质版）

合并方程的具体推导过程如下。这里要特别注意，合并方程后要求的最小的x就是a%n，所以再继续合并方程时，用x作为新的a。

• 【例题1 Strange Way to Express Integers POJ - 2891】

Elina is reading a book written by Rujia Liu, which introduces a strange way to express non-negative integers. The way is described as following:
Choose k different positive integers a1, a2, …, ak. For some non-negative m, divide it by every ai (1 ≤ i ≤ k) to find the remainder ri. If a1, a2, …, ak are properly chosen, m can be determined, then the pairs (ai, ri) can be used to express m.

“It is easy to calculate the pairs from m, ” said Elina. “But how can I find m from the pairs?”

Since Elina is new to programming, this problem is too difficult for her. Can you help her?

Input
The input contains multiple test cases. Each test cases consists of some lines.

Line 1: Contains the integer k.
Lines 2 ~ k + 1: Each contains a pair of integers ai, ri (1 ≤ i ≤ k).
Output
Output the non-negative integer m on a separate line for each test case. If there are multiple possible values, output the smallest one. If there are no possible values, output -1.

Sample Input
2
8 7
11 9
Sample Output
31
Hint
All integers in the input and the output are non-negative and can be represented by 64-bit integral types.

AC代码：

#include<cstdio>
#define ll long long

ll gcd(ll a,ll b)
{
    if(b==0)return a;
    return gcd(b,a%b);
}

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

ll inv(ll a, ll n)
{
    ll x,y;
    ll t = exgcd(a,n,x,y);
    if(t != 1)
        return -1;
    return (x%n+n)%n;
}


bool united(ll a1,ll m1,ll a2,ll m2,ll &a3,ll &m3)//合并方程
{
    ll d=gcd(m1,m2);
    ll c=a2-a1;
    if(c%d)return false;
    m3=m1*m2/d;
    ll INV=inv(m1/d,m2/d);
    ll k1=(c/d*INV)%(m2/d);
    a3=k1*m1+a1;
    return true;
}

ll crt2(ll len,ll a[],ll n[])//中国剩余定理非互质版本
{
    ll a1=a[0],n1=n[0];
    ll a2,n2;
    for(int i = 1; i < len; i++)
    {
        ll aa,nn;
        a2 = a[i],n2=n[i];
        if(!united(a1,n1,a2,n2,aa,nn))
            return -1;
        a1 = aa;
        n1 = nn;
    }
    return (a1%n1+n1)%n1;
}

ll a[1000],m[1000];


int main()
{
    int i;
    int k;
    while(scanf("%d",&k)!=EOF)
    {
        for(i = 0; i < k; i++)
            scanf("%I64d %I64d",&a[i],&m[i]);
        printf("%I64d\n",crt2(k,m,a));
    }
    return 0;
}

• 【例题2 Two Arithmetic Progressions CodeForces - 710D】
  You are given two arithmetic progressions: a1k + b1 and a2l + b2. Find the number of integers x such that L ≤ x ≤ R and x = a1k’ + b1 = a2l’ + b2, for some integers k’, l’ ≥ 0.

Input
The only line contains six integers a1, b1, a2, b2, L, R (0 < a1, a2 ≤ 2·109,  - 2·109 ≤ b1, b2, L, R ≤ 2·109, L ≤ R).

Output
Print the desired number of integers x.

Example
Input
2 0 3 3 5 21
Output
3
Input
2 4 3 0 6 17
Output
2

本题思路：
利用中国剩余定理（非互质版）合并方程组，求出解x，加减m1,m2的最小公倍数的整数倍，并且在[L,R]这个区间内的所有x的个数即为所求。本题中，由于a1,a2,k,l均是非负数，所以x必然大于等于max(b1,b2)，并且有L=max(L,max(b1,b2))。上面的结论通过数形结合也很容易得到，因为x=a1k1+b1=a2l+b2，相当于两条直线，斜率都是非负数，且横坐标取值都是正数。本题的关键是确定在[L,R]区间内的x的个数，因为所有的数据都是整数，下面的代码中用加1来调整变换误差。对于下面的代码中x-=((x-L)/T+1)*T;数形结合比较容易理解，而对于ans=(R-x)/T-(L-1-x)/T;数形结合可以理解(R-x)/T，而(L-1-x)/T的理解我加以补充，数形结合分析知，L=x,T=1的时候该式的值有可能等于-1，这一点利用数形结合很容易理解-1的含义和减去-1的用意。

AC代码如下：

#include<cstdio>
#include<algorithm>
#define ll long long

using namespace std;

ll T;

ll gcd(ll a,ll b)
{
    if(b==0)return a;
    return gcd(b,a%b);
}

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

ll inv(ll a,ll mod)
{
    ll x,y;
    ll d=exgcd(a,mod,x,y);
    if(d!=1)return -1;
    return (x%mod+mod)%mod;
}

bool united(ll a1,ll m1,ll a2,ll m2,ll &a3,ll &m3)
{
    ll d=gcd(m1,m2);
    ll c=a2-a1;
    if(c%d)return false;
    m3=m1*m2/d;
    ll Inv=inv(m1/d,m2/d);
    ll k1=(c/d*Inv)%(m2/d);
    a3=k1*m1+a1;
    return true;
}

ll crt2(ll a[],ll m[],ll n)
{
    ll a1=a[0],m1=m[0];
    ll a2,m2;
    for(ll i=1;i<n;i++)
    {
        ll aa,mm;
        a2=a[i],m2=m[i];
        if(!united(a1,m1,a2,m2,aa,mm))
            return -1;
        a1=aa,m1=mm;
    }
    T=m1;
    return (a1%m1+m1)%m1;
}

int main()
{
    ll m1,a1,m2,a2,L,R;
    ll a[5],m[5];
    while(scanf("%lld%lld%lld%lld%lld%lld",&m1,&a1,&m2,&a2,&L,&R)!=EOF)
    {
        a[0]=a1,a[1]=a2;
        m[0]=m1,m[1]=m2;
        ll x=crt2(a,m,2);
        if(x==-1)printf("0\n");
        else
        {
            ll ans=0;
            L=max(L,max(a1,a2));
            if(x>=L)x-=((x-L)/T+1)*T;
            if(L<=R)
                ans=(R-x)/T-(L-1-x)/T;
            printf("%lld\n",ans);
        }
    }
    return 0;
}