中国剩余定理解法

最新推荐文章于 2021-08-10 00:02:37 发布

mandycool

最新推荐文章于 2021-08-10 00:02:37 发布

阅读量1.2k

点赞数

CC 4.0 BY-SA版权

分类专栏：算法题

本文链接：https://blog.youkuaiyun.com/mandycool/article/details/9214527

算法题专栏收录该内容

4 篇文章

订阅专栏

问题表述

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

即 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 = c

def example():
    a = 115
    b = 367
    c = 2
    for i in range(1,20):
        print (a*i%b)

# GCD(a,b) should be 1
def 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)

## test
extra_divide(7,5,1)
extra_divide(25,7,1)
extra_divide(115,367,1)
extra_divide(367,115,1)
extra_divide(128,367,1)