中国古代数学算法的python实现(jupyter版本)

源文件

看了这篇《没有定理的中国古代数学，如何站在世界之巅？》。虽然我觉得题目很标题党，不过里面的内容很有趣啊，讲解了中国古代数学里的几个算法。由于我正在学Python，所以自然就拿来练手了。

重新用jupyter写一下 更相减损术 和 大衍求一术

更相减损术

更相减损术是用来求两个数的最大公约数.

“术曰：以少减多，更相减损，求其等也。”

这个很好写啦，读入两个数a和b，求其等也，就是一直要求到两个数相等为止。 所以用条件循环while，当a不等于b时就一直重复算。算什么呢，以少减多，就是判断一下谁大. 如果a>b, 就用a-b替换a，否则就用b-a替换b。 一直到两个数减到相等为止，就可以随便返回其中一个数作为最大公约数了。

貌似也可以扩展到有理数, 但为了防止死循环, 在循环中加入了下限的设定, 本来想是>0就可以, 但还是需要有一定的大小.

In [6]:
import numpy as np
def 更相减损术(a,b):
    while ((a != b) & (a>1e-10) & (b>1e-10)):
        if a > b:
            a=a-b
        else :
            b=b-a
    return a
In [14]:
if __name__=="__main__":
    print(更相减损术(2.5,np.pi))
    print(更相减损术(26,39))
      

3.9968028886505635e-11
13

大衍求一术

这个好帅，是求解ax≡1(mod b)。就是说a乘以x除以b所得的余数等于1。详细的解说还是看那篇公众号的文章吧。

因为打算用矩阵，所以要首先导入Numpy包，python很强大的一点就是有各种包，只要import一下就好像有了超能力。 大衍求一术要求先生成一个2*2的矩阵

$$ \begin{pmatrix} a & 1 \\ b & 0 \end{pmatrix} $$

这样的，所以先用np.matrix产生一个矩阵m，注意python的序数是从0开始的，所以m[0,0]=a, m[0,1]=1, m[1,0]=b, m[1,1]=0。

然后跟前面的更相减损术差不多，也是减来减去，区别是以行为单位来减，终点是把a的位置变成1，比大小的时候是用左边那列的元素比大小。

所以如果m[0,0]>m[1,0]，那就把上一行减去下一行m[0,:]-m[1,:]，再替换掉上面一行m[0,:]=m[0,:]-m[1,:]。反之亦然。一直重复，直到m[0,0] == 1。

通常是返回m[0,1]的数值就可以了。但有可能a,b互质，所以需要分情况讨论一下。

In [23]:
import numpy as np
def 大衍求一术(a,b):
    '''
    Solve a*x mod b ==1 with 大衍求一术
    置诸元数，两两连环求等，约奇弗约偶，遍约毕，乃变元数，皆曰定母，列右行，
    各立天元一为子，列左行，以定诸母，互乘左行之子，各得名曰衍数，
    次以各定母满去衍数，各余名曰奇数，以奇数与定母，用大衍术求一。
    '''
    m=np.matrix( [ [a,1], [b,0] ] )
    while (m[0,0] != 1 ):
        if m[0,0]>m[1,0]:
            m[0,:]=m[0,:]-m[1,:]
        else:
            m[1,:]=m[1,:]-m[0,:]
            
    if m[1,0]==1:
        return 1
    else:
        return m[0,1]
In [24]:
if __name__=="__main__":
    print(大衍求一术(65,83))

23

我的古文水平和编程水平都还不够高，不然把中国古代数学中的种种算法都写出程序也是一件很风雅的事情。

In [ ]: