中國剩余定理

【定理概述】

　　中國剩余定理（孫子定理）是中國古代求解一次同余式組的方法。是數論中一個重要定理。一元線性同余方程組問題最早可見於中國南北朝時期（公元5世紀）的數學著作《孫子算經》卷下第二十六題，叫做“物不知數”問題，原文如下：

有物不知其數，三三數之剩二，五五數之剩三，七七數之剩二。問物幾何？即，一個整數除以三余二，除以五余三，除以七余二，求這個整數。《孫子算經》中首次提到了同余方程組問題，以及以上具體問題的解法，因此在中文數學文獻中也會將中國剩余定理稱為孫子定理。

 

 

【求逆元】

逆元的含義：模p意義下，1個數a如果有逆元x，那么除以a相當於乘以x。ax≡1(mod p)。

一個數有逆元的充分必要條件是gcd(a,p)=1，此時逆元唯一存在，注意這里的唯一是指在群中唯一 。

其實如果求出一個逆元x0，那么x0 + p*k都會滿足上面的等式，但是我只取p內的正整數x0.

 【證明】 由ax≡1(mod p)等價於這樣一個方程

a*x + p*y = 1 ，或者說這個方程x有解的話x必然滿足 ax≡1(mod p)

這個方程什么時候有解呢？很顯然，當gcd(a,p) | 1時有解，所以gcd(a,p)只能是1，即a,p互質，證明完畢。

由此還可以得到一個結論，如果要求逆元，可以用擴展歐幾里得求一組解（x,y），再求出x的最小正整數（x+p）%p,

x就是a的唯一逆元。

方法1：費馬小定理求逆元，p是質數，且gcd(a,p)=1

在模為素數p的情況下，有費馬小定理 
a^p-1 ≡ 1（mod p） 

則a * a^p-2 ≡ 1（mod p）

所以a的逆元就是a^p-2，用快速冪求即可。

#include<iostream>
using namespace std;

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

long long qPow(long long a ,long long n,long long mod){
    long long ans = 1;
    //如果n的二進制最后一位是1 結果參與運算 
    //因為如果是0，那么冪運算之后該項是1,1乘任何數還是那個數，對結果不影響
    while(n > 0){
        if(n & 1)  
            ans = (ans* a) % mod;
        a = (a*a) % mod;//底數加倍 
        n >>= 1;//移位 
    }
    return ans;    
}

//
long long invEle(long long a, long long mod){
     
    //如果a 和 模數不互質則必然不存在逆元 
    if(gcd(a,mod) != 1 || mod < 2)    return -1;
    return qPow(a,mod-2,mod);
}

int main(){
    long long a,b;
    int x,y;
    while(cin>>a>>b){
        cout<<invEle(a,b)<<endl;
    }
}

 

方法2：擴展歐幾里得求逆元（高效）

typedef  long long ll;
void extgcd(ll a,ll b,ll& d,ll& x,ll& y){
    if(!b){ d=a; x=1; y=0;}
    else{ extgcd(b,a%b,d,y,x); y-=x*(a/b); }
}
ll inverse(ll a,ll n){
    ll d,x,y;
    extgcd(a,n,d,x,y);
    return d==1?(x+n)%n:-1;
}

 

方法3：歐拉定理求逆元（很少用到）

模p不是素數的時候需要用到歐拉定理

逆元打表：

typedef  long long ll;
const int N = 1e5 + 5;
int inv[N];
 
void inverse(int n, int p) {
    inv[1] = 1;
    for (int i=2; i<=n; ++i) {
        inv[i] = (ll) (p - p / i) * inv[p%i] % p;
    }
}

 

【解方程組】

根據定理概述以及解法，得到以下方法

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);
        ans = (ans + Mi * x * a[i]) % M;
    }
    if(ans < 0) ans += M;
    return ans;
}

 

【擴展中國剩余定理】

當模數mi兩兩互質時有以上解法，當模數不確定是否兩兩互質呢？

 

摘自博客:https://blog.csdn.net/acdreamers/article/details/8050018

這種情況就采用兩兩合並的思想，假設要合並如下兩個方程

     

那么得到

      

在利用擴展歐幾里得算法解出的最小正整數解，再帶入

      

得到后合並為一個方程的結果為

       

這樣一直合並下去，最終可以求得同余方程組的解。

 

#include <iostream>
#include <string.h>
#include <stdio.h>
using namespace std;
typedef long long LL;
const int N = 1005;
LL a[N], m[N];
LL gcd(LL a,LL b)
{
    return b? gcd(b, a % b) : a;
}
void extend_Euclid(LL a, LL b, LL &x, LL &y)
{
    if(b == 0)
    {
        x = 1;
        y = 0;
        return;
    }
    extend_Euclid(b, a % b, x, y);
    LL tmp = x;
    x = y;
    y = tmp - (a / b) * y;
}
LL Inv(LL a, LL b)
{
    LL d = gcd(a, b);
    if(d != 1) return -1;
    LL x, y;
    extend_Euclid(a, b, x, y);
    return (x % b + b) % b;
}
bool merge(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;
    c = (c % m2 + m2) % m2;
    m1 /= d;
    m2 /= d;
    c /= d;
    c *= Inv(m1, m2);
    c %= m2;
    c *= m1 * d;
    c += a1;
    m3 = m1 * m2 * d;
    a3 = (c % m3 + m3) % m3;
    return true;
}
LL CRT(LL a[], LL m[], int n)
{
    LL a1 = a[1];
    LL m1 = m[1];
    for(int i=2; i<=n; i++)
    {
        LL a2 = a[i];
        LL m2 = m[i];
        LL m3, a3;
        if(!merge(a1, m1, a2, m2, a3, m3))
            return -1;
        a1 = a3;
        m1 = m3;
    }
    return (a1 % m1 + m1) % m1;
}
int main()
{
    int n;
    while(scanf("%d",&n)!=EOF)
    {
        for(int i=1; i<=n; i++)
            scanf("%I64d%I64d",&m[i], &a[i]);
        LL ans = CRT(a, m, n);
        printf("%I64d\n",ans);
    }
    return 0;
}