【数论】C++实现中国剩余定理求一次同余方程组

很久很久以前，《孙子算经》里，有一道“物不知数”题：

有物不知其数，三三数之剩二，五五数之剩三，七七数之剩二，问物几何？

就是说，有一堆东西，三个三个数剩余两个，五个五个数剩余三个，七个七个数剩余两个，请问有多少个？

粗涉数论的人应该知道“同余”和“同余式（组）”的概念。设这堆东西的数量为x，则原问题可表示为同余式组

　　　　x ≡ 2  (mod 3),

　　　　x ≡ 3  (mod 5),

　　　　x ≡ 2  (mod 7).

原问题在《孙子算经》里的解释为：“三三数之剩二，置一百四十；五五数之剩三，置六十三；七七数之剩二，置三十。并之，得二百三十三，以二百十减之，即得。答曰：二十三。”（引用自百度百科：物不知数）

即：先求被3除余2，并能同时被5、7整除的数，这样的数最小是35；

再求被5除余3，并能同时被3、7整除的数，这样的数最小是63；

然后求被7除余2，并能同时被3、5整除的数，这样的数最小是30。（引用自百度百科）

他们相加得128，但是它不是最小的。可以求它除以（3*5*7=105）的余数，也就是答案23.

同理，所有模105余23的数都是原问题的答案：

　　　　x ≡ 23  (mod 105).

这种问题有一般的推广形式，就是很重要的中国剩余定理：

设m1,m2,...,mk是k个两两互素的正整数，则对任意的整数b₁,b₂,...,b_k,同余式组

　　x ≡ b₁  (mod m₁),

　　x ≡ b₂  (mod m₂),

　　......

　　x ≡ b_k  (mod m_k)

一定有解，且解是唯一的。

若令

　　m = m₁m₂...m_k ,  m / m_i = M_i , M_i^'M_i ≡ 1  (mod m_i) 【即M_i^' 是M_i的模m_i逆元】

则

　　x ≡ b₁M_1^'M₁ + b₂M_2^'M₂ + ... + b_kM_k^'M_k  (mod m).

利用万能的C++，我们可以设计出运用中国剩余定理求x的程序。

首先是开头部分：

#include<iostream>
#include<cstdio>
#include<cmath>
#define MAXN 10000
using namespace std;
int TotalNums;//同余方程组中同余方程的数量
int b[MAXN],m[MAXN],M[MAXN],revM[MAXN],kM;//b[i]是第i个同余式右边，m[i]是第i个同余式的模，kM是所有模的乘积，M[i]是kM与m[i]的商，revM[i]是M[i]模m[i]的逆元。

对于程序，这几个函数是必要的：

int gcd(int a,int b){return (a%b==0)?b:gcd(b,a%b);}//欧几里得递归求最大公约数
int lcm(int a,int b){return a*b/gcd(a,b);}//根据两个数的最大公约数和最小公倍数的乘积等于这两个数的乘积，求出最小公倍数
void exgcd(int,int,int&,int&);//广义欧几里得算法
int getRev(int);//获得逆元，参数表示这是第几个同余式的情况。
int getProduct(int[]);//数组求积
bool isEachPrime();//判断数列中的数是否互素
int _x[MAXN],_y[MAXN];//临时变量，供exgcd函数后两个引用参数使用

void exgcd(int a,int b,int& x,int& y)//求满足xa+yb=(a,b)的整数x,y.由于不好返回两个值，因此存储在参数x和y中，因此要使用引用。
{
    if(b==0)
    {
        x=1;y=0;
        return;
    }
    exgcd(b,a%b,x,y);//具体请参阅“Bezout恒等式”
    int p=x;
    x=y;
    y=p-a/b*y;
}//其实这段代码我是背下来的，我也不太理解o(╥﹏╥)o
int getRev(int i)
{
    exgcd(M[i],m[i],_x[i],_y[i]);//对于满足xa+yb=(a,b)的同余式，如果b=1，则上述x是满足ax≡1(mod b)的解。可以证明当(a,b)=1时，这个同余式有唯一解。
    return (_x[i]+m[i])%m[i];//得到最小值
}
int getProduct(int a[])//获得数组a的元素乘积
{
    int p=1;
    for(int i=1;i<=TotalNums;i++)
    {
        p*=a[i];
    }
    return p;
}
bool isEachPrime()//判断数组m的元素是否两两互素。如果一个数组的元素都是两两互素，那么它们的最小公倍数是它们的乘积。求数组最小公倍数的方法：先求出前两个元素的最小公倍数，再求出这个最小公倍数和第三个元素的最小公倍数，以此类推。
{
    int p=m[1];
    for(int i=2;i<=TotalNums;i++)
    {
        p=lcm(m[i],p);
    }
    if(p!=getProduct(m)) return false;
    else return true;
}

最后是主程序：

int main()
{
    cout<<"请输入同余式数量：";
    cin>>TotalNums;
    cout<<"格式："<<endl;
    for(int i=1;i<=TotalNums;i++)
    {
        cout<<"x≡b_"<<i<<"(mod m_"<<i<<")"<<endl;
    }
    for(int i=1;i<=TotalNums;i++)
    {
        cout<<"请输入b_"<<i<<":";
        cin>>b[i];
        cout<<"请输入m_"<<i<<":";
        cin>>m[i];
　　　　 b[i]%=m[i];//防止出现同余式右边比模大的情况
    }
    if(!isEachPrime())//无解
    {
        cout<<"m没有两两互素，原方程组无解。"<<endl;
        return 0;
    }
    else//有解
    {
        kM=getProduct(m);//获得模的乘积
        for(int i=1;i<=TotalNums;i++)
        {
            M[i]=kM/m[i];//求得M[i]
            revM[i]=getRev(i);//获得逆元
        }
        int result=0;
        for(int i=1;i<=TotalNums;i++)
        {
            result+=(b[i]*revM[i]*M[i]);//根据中国剩余定理的定义
        }
        result%=kM;//获得最终结果
        cout<<"结果：
x≡"<<result<<"(mod "<<kM<<")";
        return 0;
    }
}