中国剩余定理(孙子定理)

“中国剩余定理”是公元5-6世纪、我国南北朝时期的一部著名算术著作《孙子算经》中的一个“物不知数”的解法问题:今有物不知其数,三三数之剩二,

五五数之剩三,七七数之剩二。问物几何?答曰:二十三。

 

 1.除数两两互质的情况

根据上面我们可以得到一组式子:

 

X ≡ 2  (mod 3)

X ≡ 3  (mod 5) 

X ≡ 2  (mod 7)

即:

X % 3 = 2

X % 5 = 3

X % 7 = 2

那么我们设   X = t1 * a + t2 * b + t3 * c

继续设

a ≡ 1  (mod 3)                          b ≡ 0  (mod 3)                         c ≡ 0  (mod 3) 

a ≡ 0  (mod 5)                          b ≡ 1  (mod 5)                         c ≡ 0  (mod 5) 

a ≡ 0  (mod 7)                          b ≡ 0  (mod 7)                         c ≡ 1  (mod 7)

a = LCM(5,7) * p(a%3=1)     b = LCM(3,7)*p(b%5=1)       c = LCM(3,5)*p(c%7=1)

   = 35 * 2 = 70                           = 21 * 1 = 21                           = 15 * 1 = 15

那么我们就可以根据相对应的颜色得到t1 = 2, t2 = 3, t3 = 2

所以 X = 2a + 3b + 2c = 2 * 70 + 3 * 21 + 2 * 15 = 233

233是X的一个解,通解是 X = 233 + gcd(3,5,7)*k = 233 + 105*k

这些a、b、c(下面同称为Ni)该如何求呢(除数3、5、7统称为mi)

X ≡ y1  (mod m1)

X ≡ y2  (mod m2) 

X ≡ y3  (mod m3)

...

X ≡ an  (mod mn)

设X = N1 * y1 + N2 *y2 + N3 * y3 + ... + Nn * yn;

设  Ni能够被m1、m2、m3、.... 、mi-1、 mi + 1、... 、mn整除,但就是不能被mi整除且Ni%mi = 1

设  Ni = LCM(m1, m2, m3, ..., mn) / mi * x  = mi * y + 1, 要求Ni就必须得求x

上面式子可以转化为  LCM(m1, m2, m3, ..., mn) / mi * x + (- mi) * y = 1

转化之后的式子就是扩展欧几里德了, 这样就可以求出x

 

模板:

int CRT(int a[], int b[], int n)
{//a[i]表示除数,b[i]表示余数
    int M = 1, N, ans = 0, x, y;
    for(int i = 0 ; i < n ; i++)
        M *= a[i];//因为除数两两互质,所以最小公倍数就是其乘积
    for(int i = 0 ; i < n ; i++)
    {
        N = M / a[i];
        gcd(N, a[i], x, y);//利用扩展欧几里德求x
        ans += N * x * b[i];
    }
    return ans % M;
}

 

2.除数不两两互质

 当除数不两两互质时上面的做法就不正确了,那么应该怎们做呢

我们可以两个的合并方程,如下面两个方程:

X ≡ b1  (mod n1)

X ≡ b2  (mod n2) 

 

X = b1 + n1 * k1;  (1)

X = b2 + n2 * k2;

可将上面两个方程画等: b1 + n1 * k1 = b2 + n2 * k2  ==> n1 * k1 = b2 - b1 + n2 * k2      (2)

设r = gcd(n1, n2),左右两边同时除以r得:

b1/r + n1*k1/r = b2/r + n2*k2/r

转化为 n1*k1/r = b2/r - b1/r + n2*k2/r    

所以 n1*k1/r = (b2 - b1)/r(mod n2/r)

==> n1*k1 = (b2 - b1)(mod n2/r)

k1 = (b2-b1)/n1(mod n2/r)

令k’ = (b2 - b1) / n1

所以 k1 = k’(mod n2/r)

k1 = k’ + n2/r * C

将k1代入(1)式中:

X = b1 + n1 * k1 = b1 + n1 * (k’  + n2/r * C ) = b1 + n1 * k’ + n1 * n2 / r * C

由(2)得:

       n1 * k1 = (b2 - b1)(mod n2)

所以 n1 * k1 + n2*k2 = b2 - b1

可用扩展欧几里德定理来求k1

t = n2 / r

k’ = (k1 % t + t) % t

其实k与k'是等价的

详解看下面链接:

http://972169909-qq-com.iteye.com/blog/1266328

 

模板:

 

ll CRT2(ll n[], ll b[], ll m)
{
    int f = 0;
    ll n1 = n[0], n2, b1 = b[0], b2, c, t, k, x, y;
    for(ll i = 1 ; i < m ; i++)
    {
        n2 = n[i];
        b2 = b[i];
        c = b2 - b1;
        gcd(n1, n2, x, y);//扩展欧几里德
        if(c % r != 0)//无解
        {
            f = 1;
            break;
        }
        k = c / r * x;//扩展欧几里德求得k
        t = n2 / r;
        k = (k % t + t) % t;//这个k就是上面的k’
        b1 = b1 + n1 * k;
        n1 = n1 * t;
    }
    if(f == 1)
        return -1;
    return b1;
}

 

posted @ 2015-11-12 16:25  午夜阳光~  阅读(713)  评论(0编辑  收藏  举报