中国剩余定理解法

问题表述

 

《孙子算经》卷下第二十六问：“今有物，不知其数。三、三数之，剩二；五、五数之，剩三；七、七数之，剩二。问物几何？答曰：二十三。”

即 x%3=2；x%5=3；x%7=2；  》》  x=23+3*5*7*k

表述成定理为：

给m个两两互质的数，其乘积为P，若已知x<P，且知x与这m个数的模，则x唯一确定。

以m=2为例，设为a，b，令x=1~a*b；

x mod a 的结果为 (0,1,...,a-1) 的循环；

x mod b 的结果为 (0,1,...,b-1) 的循环；

则(x%a, x%b)的循环周期为P=a*b， 在P范围内，是唯一确定的。

 

 

问题求解

m=3，分别为a，b，c，如何求x？

brute-force：

     求出P内满足条件 i (x%ai=ci) 的x的集合，求交。

clever：

     构造出一个满足条件的x，不必小于P，x%P即为所求。具体的，构造t1为ab的倍数，且模c等于x模c，同理构造出t2，t3，则 x=(t1+t2+t3)%P。

这引出了下一个问题：a*x % b = c ，求x？

     事实上，只有当a与b互质，方程有解；因为其等价于 X=a*x，X%a=0，X%b=c，由中国剩余定理可知其有解。if (a*x)% b= 1, (a*c*x)%b = c，因此也等价于求 (a*x)% b= 1 。

     即求解 a的多少倍 - b的多少倍 = 1，用扩展的辗转相除法。

     如a=5，b=7：

          7/5 = 1 ...2，(5,7) 的系数为 (1,1)，剩下 5,2

          5/2 = 2 ...1，(5,2) 的系数为 (1,2)，于是(5,7)的系数为(3,2)

          5*3 - 7*2 = 1，求得 x = 3。

     该扩展辗转相除法还是相当精妙的！

     若 a和b的最大公约数d不为1，则只能求解 a*x % b = d 。

 

参考：Matrix67大牛的博文 ---“跨越千年的RSA算法”

 

 

代码实现

 

python写的代码，从上述思路到实现还是有一定距离的，^ ^

# solve a*x % b = c ;# a and b are co-prime ;# solve it equals to a*x % b = 1 ,# cause if (a*x)% b= 1, (a*c*x)%b = cdef example():    a = 115    b = 367    c = 2    for i in range(1,20):        print (a*i%b)# GCD(a,b) should be 1def extra_divide(a,b,c):    if not c<b:        print 'b should be greater than c'        return -1    # assume a > b    aa = a    bb = b    exchanged = False    if a<b:        exchanged = True        t = a        a = b        b = t        rec_array = []        r = 0    while r!=1:        r = a%b        d = a/b        rec_array.append( 1 )        rec_array.append( d )        a = b        b = r    #print rec_array    # this is the most important part, take 25,7 for example    # 25 % 7, 7 % 4, 4 % 3    #  1   3, 1   1, 1   1;  one 3 equals one 7 and one 4, so    #  1   3, 1   2;         two 4 equals two 25 and six 7, so    #  2   7;            for i in range(0,len(rec_array)-2,2):        j = len(rec_array) - i - 1        rec_array[j-2] *= rec_array[j]        rec_array[j-3] *= rec_array[j]        rec_array[j-2] += rec_array[j-1]    k1 = rec_array[0]    k2 = rec_array[1]    if exchanged:        t = k1        k1 = k2        k2 = t        # we want ax%b=1, now b%(ax)=1, so play a magic here    # k1*a + 1 = k2*b, also b*a = a*b, take a difference    # we get (b-k1)*a -1 = (a-k2)*b, bingo!    if aa*k1 < bb*k2:        k1 = bb-k1        k2 = aa-k2    print 'x='+str(k1), 'a*x='+str(k1*aa), 'b*y='+str(k2*bb)## testextra_divide(7,5,1)extra_divide(25,7,1)extra_divide(115,367,1)extra_divide(367,115,1)extra_divide(128,367,1)