月度归档:2013年03月

你们的都是我们的

想必看本博客的人,一定都经常混迹于网络,都听说过韩国高丽棒子的各种“无耻言论”。比如“韩国为汉字申遗”,“韩国为风水申遗”,“韩国教授考证孙中山是韩国人”,诸如此类。好笑的同时也不禁让我们产生疑问,韩国人真的这么无聊可耻么?

x稍加搜索其实就能发现,这些新闻都是中国的一家小报,两三个人编出来的。这家报纸叫《新快报》。这家报纸原本没什么人买,报纸主编日日为销量发愁,某天一编辑献计,编造新闻抹黑韩国,主编大悦。随后一则又一则非常雷同的新闻就这么炮制出来了,内容大概就是“中国人的都是韩国人的”,比如孙中山、孔子、风水、汉字。这些新闻的影响力大家都见识过了,各大门户媒体纷纷转载,网民们津津乐道。影响力大到连韩国报纸、韩国大学、韩国驻华使馆都出来纷纷辟谣,可悲的是,辟谣消息却没有媒体报道,也没什么人看到,反倒是谣言知道的人越来越多。原本以为只有愤青才会为这个新闻买单,结果连当初的小报主编都没想到,新闻变成网民心目中的“即成事实”了。平日茶余饭后必谈韩国人的笑话,大家欢欣鼓舞喜笑颜开。同样高兴的还有《新快报》报纸的主编,销量上去了,报纸影响力也提升了。如此被冤枉的韩国人的笑料变成了中国人的笑料,上演了一出让人啼笑皆非的人间闹剧。

殊不知,韩国人的“你们的都是我们的”都是中国人自己编造的。而某个历史更悠久的国家在几百年前就开始了“你们的都是我们的”的考证。

故事要从一个原本普普通通的意大利人说起。

这个人是意大利省长的长子,可以算是比官二代更厉害的太子党吧。从小就特别聪明,不管什么书看一遍就可以过目不忘。原本省长父亲是想让儿子好好念大学以后子承父业做官,但是这个人在朋友的引导下皈依了天主教,辍学去念了神学,投靠了伟大的主的怀抱,把父亲气的够呛。

某一个阳关灿烂的下午,教皇召见了他。然后大概发生了下面对话:

“你可愿为将吾主之光辉洒向世界的每一个阴暗角落而奉献一生?孩子”

“这是我毕生之心愿,殿下!”,他虔诚的匍匐在地上。

“地球的另一边,在遥远的东方,有一个庞大的帝国,无信而又愚昧,你可愿去帮助吾主驱除黑暗?”

“如您所愿!在吾主之光辉沐浴下,黑暗必将驱散!”

于是,这个年仅25岁的意大利小伙儿,为了毕生的信仰,义无反顾的离开了家乡,去了一个极其遥远又极其陌生的地方,进行极其艰难的传教事业。遥远到远在万里之外地球的另一端;陌生到需要面对一个语言不通没法交流的封闭群体;艰难到一旦说起传教就会变成过街老鼠人人喊打。

他去了遥远的东方一直到死都没有回过家乡意大利,在那里度过了后半辈子。

li这个人就是意大利传教士马蒂奥·利奇。他还有一个中国人极其熟悉的中文名字:利玛窦。

十六世纪的旅行,可不像今天那样,可以随便在某个周末早上跑到世界某个角落去喂鸽子。利玛窦和他的老师罗明坚以及另外十四位神父,先去了葡萄牙的里斯本,等了几个月等到前往东方的船只,然后在又小又闷的船舱里住了半年,好在没有发生传染病和瘟疫,但是等他们到了印度果阿,大家疲惫不堪,不得不又养病休息了几个月。

过了4年,利玛窦终于接到了命令,和老师罗明坚一起前往澳门。当时同行的神父都觉得在中国传教希望十分渺茫,因为明朝甚至规定所有西夷(洋鬼子)没有”路照”禁止进入广州。但是他和他的老师依然在澳门艰苦的学习着中文。他的老师当时已经36岁了,学习中文十分的吃力。当时可没有字典或者翻译!极其封闭的中国甚至找不到一个会外语的人,学习中文的唯一办法就是把东西画下来,找当地人问发音,然后一个词一个词的模仿。

又过了一年,师徒俩总算找到机会,联系上了肇庆知府,谎称自己是来自天竺的僧人(当时中国人觉得天竺就是了不得远的地方了),于是混进了广州。在广州肇庆,他们一边筹钱,一边兴建在中国的第一座天主教堂:仙花寺。当然,这是披着佛教外衣的天主教堂。然而好景不长,6年之后,广州总督经过人举报发现了他们传的不是佛教而是基督教,于是被驱逐出了肇庆,连好不容易造好的西式教堂也被广州总督霸占。利玛窦开始意识到,如果没有中国皇帝的允许,在中国传教是一件多么多么艰难的事情。

这时候他的老师罗明坚也回了意大利,只有他留下来继续传教。之后的12年,利玛窦认准了一个死理,就是必须去北京定居传教,从中枢影响这个国家。利玛窦找了种种机会去北京,期间遭遇了普通人难以想象的各种艰难险阻。第一次去,最远只到了南京,回来大病一场。第二次去,因为遇到了日本和朝鲜打仗(不可抗拒之力)又被迫回来了。第三次去,又被一个贪婪的太监抓了起来,打了一顿,然后关押了1年。后来利玛窦总算学会了中国官员的那套,给银子上下打点各种达官贵人,带着大量的贡品:自鸣钟、圣经、万国图志、大西洋琴,面见了皇帝。皇帝看到这么多新鲜玩意儿龙颜大悦,准许了利玛窦永久居住北京。

可见几百年前一个北京户口对外国友人的吸引力之大!12年光阴,无数金银财宝、稀奇玩意儿,各级官员上下打点,经历了强盗劫匪,疾病困苦,朝鲜战争,宦官发难,牢狱之灾,最后户口总算办下来了!

利玛窦期间为了能融入中国的仕族圈,不停的宣扬西方的科学知识,西方的文化。外加利玛窦的记忆力超群,他期间写下了很多书籍,同时还把中国的四书五经翻译成了拉丁文。

map利玛窦此举其实是无心,只是一个为了融入这个社会的手段。他带到东方的东西也不是什么稀罕玩意儿——世界地图:西方满大街都是了;天文学:哥白尼都死了几百年了;几何原本:西方的初级课本而已,是1800多年前的欧几里得写的。但是这些书籍和知识对当时的中国知识分子带来的冲击是非常巨大的,因为这些是第一次传进古老的中国。中国人第一次知道了原来地球是圆的而不是方的!地球是绕着太阳转的!原来阿尔朱巴尔(algebra 代数学的另一种翻译)这么高深复杂!原来几何是这么的深奥!

你一直在用诺基亚黑白屏手机,忽然有一天你发现邻居在用先进的iPhone5,比你的领先几个时代,你会作何反应?

  1. 努力奋斗好好赚钱卖肾买一个iPhone5去。
  2. 诺基亚就是好用,就是耐摔,我就喜欢用。你们爱用什么用什么,我不稀罕。
  3. 好贼子!这不就是我去年丢的那个破西门子么?你别以为你给染个颜色加个边框我就不认识了?

当时的中国知识分子从利玛窦那里看见了自鸣钟、世界地图、几何原本之后,大体也就是这三种反应。不过可悲的是,自尊心极强的古代知识分子们,大多数都选择了第三种办法来应对。

其中最正常而且明智的人莫过于徐光启了,他选择的是第一种办法。徐光启在看到这么多先进东西之后,喊出了震耳欲聋的口号:“如果想超越,一定得先学懂。如果要学懂,一定要先翻译!”。徐光启带领其他一些知识分子开始第一次主动的学起了外语,翻译起了外文书。徐光启带头翻译了利玛窦带来的18卷几何原本里的前6卷,还给了评价说:“所有的中国人都应该来学几何原本,如果这本书能学精通,那什么学问都可以精通,如果这本书可以学下来,那什么学问都可以学下来”。然后呢?没有然后了。徐光启的话被绝大多数中国穷酸儒们置之高阁,中国知识分子随后开始了轰轰烈烈的鸵鸟学潮。

0912291105_5让寒窗苦读的酸儒们抛弃四书五经,来看几何原本、代数、天文、农业、政治法律学,必定是极其困难的。但是他们最擅长的是研究古书,从古书海里找出一丝半句的牵强附会的依据来证明说明什么。其中黄宗羲、王锡阐无疑是个中高手。据他们的考证,这些个代数几何,都是从中国传过去的!黄宗羲说:“中国的古籍有《周髀算经》,里面有勾股之术,这不就是几何么,我们现在失传了,但是被西人窃取了而已”。王锡阐更有意思,说:“屈原的诗《天问》里就说过圜则九重,你们看了这个才开始研究天文学的”。还有更可怕的,比如“《诗经》里就提过民主的概念,你们西人不过是看了有启发而已”,“自行船千里镜不过是看了我们古书,仿照诸葛丞相的木牛流马而已”。各种奇怪的考证不一一列举,总之华夏文明依然是天下第一!你们现在虽然先进,不过是偷了我们老祖宗的!

这些知识分子不去学习一句半句外语,不去继续翻译外国人的书,而是一头钻进了中国浩如烟海的古书里,从里面找西人窃取古人的新证据!

这些人还包括了赫赫有名的方以智、顾炎武、王夫之、梅文鼎、颜元、毛奇龄、胡渭、阎若璩,以及后来的戴震。不过这些人加在一起也没有一个人的力量大,这个人就是康熙帝。

康熙很喜欢西学,他接触过很多传教士,甚至自己对代数有一点研究。某一次康熙下江南,接见了当时的学派领袖梅文鼎,两人在运河上泛舟三日,相谈甚欢。康熙给西学下了定论:“西洋的算法确实是厉害,不过以前也是中国的算法,后来传到西洋才叫阿尔朱巴尔(代数),其实是源自中国。”康熙还自己写了一本《三角形论》,大谈西学东源之说。

有了皇帝带头,手下当然更积极的找证据了。清朝期间找到的各种文献古籍无数,举不胜举。到了乾隆编著四库全书的时候,一个叫戴震的更加荒唐,写了一本《勾股割圆术》,这其实是一本西学的书,但是他把西学的词全部用古语替换,号称三角学都是从中国古书中推导出来的。

这个普通的虔诚传教士利玛窦,估计做梦都没想到,他最初为了传教带去中国的几样普通玩意儿,给东方这个腐朽的帝国带去如此大的震动。然而这潭水太深了,一些现代科学的火种很快就被淹没。世界地图偶尔看看图个乐子感慨一下原来世界这么大,自鸣钟被放进皇宫里没有维修很快坏了,那些科学著作都被考证原来是抄袭咱们老祖宗,没啥大不了的了。

是不是只有愚昧的清朝学者如此,到了近代这个“你们的都是我们的”理论就没法落脚了呢?其实不然。

中华文明自以为天下第一的不求上进,使得这个文明在之后的几百年一直在被欺负挨打,惨痛的教训让之后的中国人舔了几代人的伤口。在真的实现现代化接触了世界文明之后,才发现其实很多东西我们依然被欺骗了很久。

其中不得不谈到一个英国汉学家李约瑟。之前那么多古人编造的西学东源故事,恐怕还没这个人一个人加起来的多。1942年,中国对日战败,绝大多数中国领土沦陷,蒋介石逃到了重庆,全国一片哀号。这时候汉学家李约瑟出现了。在蒋介石的资助下,李约瑟在重庆1943年搞出了“四大发明”的成果,在抗战时期大大提升了广大军民的民族自豪感和斗志。张召忠将军评价李约瑟的这个成果”抵得上共和国卫队10个师”!

cover_715099_556_big1947年李约瑟回英国开始专心写书,1954年发布了一本奇书《中国科学技术史》,帮助中国人民总结了过去2000年中国所有的科学技术,让中国人大吃一惊,原来我们的老祖宗这么厉害过!

从此四大发明成为一个可笑的玩意儿,中国人关起国门来孤芳自赏引以为豪的东西。外国人从来也没有承认或者认可过,只有在中国官方肯定并且玩命儿似的宣传。

李约瑟的书以及四大发明的问题,请大家搜索一下必然能找到很多答案,这里略举几个例子,比如沈括《梦溪笔谈》只是一本奇闻录,活字印刷在中国从来只是传说故事。《中国科学技术史》说,中国人最先发明了机器人,因为汉朝就有关于陈平造木头人会跳舞的传说。还有,中国人最先发明了直升飞机,因为有古文记载某种竹蜻蜓能载人的故事。再比如二进制是易经里早有叙述(八卦,64卦能叫二进制么?),莱布尼茨可能是看了《易经》才得到启发。其中反复引用培根的话证明四大发明,但是培根和马克思也从来没说过四大发明是中国的,相反马克思说到中国是一个科学技术的贫瘠之地。

后来又过了很多年,另一个作者罗伯特·坦普尔在1986年,完全遵从李约瑟的大作,搞出来一本《中国的一百个世界第一》,成为中国儿童的科学启蒙教育书,这真是非常可悲的事情。

中华民族有着璀璨的历史,不用妄自菲薄,也不能随便拔高古人。现在被写进教科书里被吹嘘过头的各种古人的科学发明成果实在是数不胜数,比如指南车、地震仪、司南,等等(这些东西的作假证据就不一一给链接了)。中国古人确实是科技贫瘠,但是古人的手工业水平和制造业水平还是相当的高的。如果只是为了证明古代科技也曾经辉煌过,而胡编乱造,找一些野史传说来瞎凑自欺欺人,然后培养国人所谓自豪感,无疑是非常非常可悲的。

中华文明的未来需要我们去努力创造,不迷信古人,不吹嘘历史。

不要学被我们冤枉的韩国人——“你们的都是我们的”。

 

圆周率之随笔(3)

今天是圆周率节,赶工怒发博文!

在古代,为了计算高精度的圆周率,数学家们往往要用纸和笔辛苦计算几个月甚至几年。而今天我们点点鼠标就可以悠闲的分分秒把圆周率计算到几百万甚至上亿位。这一章我们来看看,在现代计算机上,如何在最短时间里计算出最长的圆周率。

不少人想到计算圆周率可能会联想到Superpi这个邪恶小程序。这个程序是日本人金田康正在1995年写的,可以在家用电脑Windows系统上把圆周率计算到3000万位。国内很多的装机商攒机商,都用这个程序来测试新装机的超频性能。这其实是很不科学的。首先这个程序是1995年开发的,距离今天已经有快20年,不光算法过时,编译器也已经过时,很多新CPU的特性都没有利用上。第二,这个程序是用单线程开发的,现代CPU动不动就是4核8核,单线程的计算程序早已经无法反应CPU真实能力。第三,单一的用Superpi做超频测试,也容易让硬件厂商进行有针对性的作弊。现在其实有比较新的 Hyper PI 可以做为Superpi的替代品。

Superpi 和 Hyper PI 远远不是现在在家用电脑上计算圆周率的主流和先锋。他们采用的高斯-勒让德算法也早已过时,有更好的可以替代。

由于圆周率的计算算法实在是五花八门,如果都总结下来的话,作者可能可以骗几十章博客更新,读者可能也懒得再看直接关窗口。那么这里只说一个目前为止最先进和最好的算法。这就是楚德诺夫斯基(Chudnovsky)算法。

楚德诺夫斯基是生于乌克兰的美国数学家,1992年纽约人杂志把楚德诺夫斯基评为世界上最厉害的数学家之一。1987年他发明的计算圆周率的楚德诺夫斯基算法,一直到今天仍然保持着圆周率的最佳计算记录。

该算法是如何来的,请看专业论文,这里直接给出公式:

826dc7788dba249ee86fc0135e06b035

看上去超级复杂是不是?其实不然,通过简化,代码只有短短的几行。该算法的最大特色就是,整个计算过程只用最后做一次除法和一次平方根运算。其余全都是加法和乘法!

a = lambda i: (1,1,13591409) if i==0 else ((6*i-5)*(2*i-1)*(6*i-1),i**3*10939058860032000,
        (6*i-5)*(2*i-1)*(6*i-1)*(13591409+545140134*i)*(-1 if i%2 else 1))
b = lambda x,y: (x[0]*y[0], x[1]*y[1], y[1]*x[2]+x[0]*y[2])
bs = lambda i,j: a(i) if j-i==1 else b(bs(i,(i+j)/2),bs((i+j)/2,j))
sqrtC = lambda digits: isqrt(1823176476672000*10**(2*digits))
def isqrt(n):
    x,d = 2**(n.bit_length()-1),1
    while d:
        d = (n/x-x)/2
        x += d
    return x

def pi(digits):
    N = int(digits/14.1816474627254776555+1)
    P, Q, T = bs(0, N)
    return (Q*sqrtC(digits)) / T

这段代码在Python 2.X上可执行,运行print pi(100)看看?圆周率的前100位就出来咯。

对于这段代码稍作一些解释:上文看似复杂的公式,其实被拆成了两个函数:a 和 b。整个圆周率计算被简化成了一个mapreduce操作:reduce(b, map(a, range(N)),也就是 N 次 a 运算和 N-1 次 b 运算。为了继续优化,我们设计了 bs 函数,做的是二分切割,将每次计算区域切割两份,然后分头计算,最后做归并。这个计算量其实是一样的,还是 N 次 a 运算和 N-1 次的 b 运算。那么二分切割为什么会运行的更快呢?其主要原因是进行二分切割可以将高精度的大数乘法尽量拖后,执行的层数越低,数字的长度才成指数倍提高。

楚德诺夫斯基算法号称是目前计算圆周率最快的算法,然而这段代码在笔者的电脑上运行计算 pi(10000) 居然要 60 秒之久。问题出在哪呢?经过简单的 benchmark 我们发现,其主要时间消耗在了 isqrt 函数里。由于 python 默认没有高精度整数的平方根操作,所以这个用牛顿拉普森法模拟的平方根操作很慢。

为了快速求平方根,下面我们祭出高精度运算的大杀器 gmp。同时采用python可以用的 gmpy 来加速平方根操作。经过修改的代码如下:

from gmpy2 import *

a = lambda i: (1,1,13591409) if i==0 else (mpz((6*i-5)*(2*i-1)*(6*i-1)),mpz(i)**3*10939058860032000,
        mpz((6*i-5)*(2*i-1)*(6*i-1))*(13591409+545140134*i)*(-1 if i&1 else 1))
b = lambda x,y: (x[0]*y[0], x[1]*y[1], y[1]*x[2]+x[0]*y[2])
bs = lambda i,j: a(i) if j-i==1 else b(bs(i,(i+j)/2),bs((i+j)/2,j))
sqrtC = lambda digits: isqrt(1823176476672000*mpz(10)**(2*digits))

def pi(digits):
    N = int(digits/14.1816474627254776555+1)
    P, Q, T = bs(0, N)
    return (Q*sqrtC(digits)) / T

经过修改的代码也更短更简洁了,只有短短的10行。让我们测试一下速度如何?(需要先安装 gmpy2 的 python 包才可执行)

位数 时间
10,000 0.04s
100,000 0.11s
1,000,000 1.07s
10,000,000 17.76s

20秒内,在笔记本上轻松将圆周率计算到了 1000 万位!1987年在巨型超级计算机上才创造的圆周率世界纪录,在今天的任何一台个人电脑上,都可以轻松碾压了。

然而我们还不满足。这个代码依然没有把楚德诺夫斯基算法的优势发挥出来。楚德诺夫斯基算法最大的优势是可以很好的并行化。在多核时代,单线程运行的代码实在没法体现出CPU的优势出来。

然而由于python的GIL的设计问题,多线程的代码无法很好的利用多核。因此我们用多进程的方法实现了楚德诺夫斯基算法的并行化。下面是在4核CPU的笔记本的运行环境下比较不同进程数量的程序性能。计算位数依然是1000万位:

计算核心(进程数) 时间
单核版本 17.76s
1 18.58s
2 12.66s
3 10.39s
4 9.88s
5 9.79s
6 10.17s

可以看到,通过并行化代码,利用CPU的多核,我们将性能提升了1倍!在计算核心继续增加超过CPU本身的4核之后,性能不升反降,这也是预料之中的事。

有人会问,为什么用4个计算核心提升的性能不是4倍?这里其实牵扯到并行化算法中比较复杂的问题。4核心最理想的性能当然是4倍的提升,但是由于算法本身限制,以及任务细分之后之间的依赖关系,计算核心之间的数据交换,导致最理想情况很难做到。算法中的任务划分也非常有学问,分的太粗会对多内核利用率不高,分的太细又会大大增加进程/线程间的通讯和同步的消耗。这里就不做更深入探讨了。

在这里我们不得不提到 gmp 官方给的一个 C 语言实现的楚德诺夫斯基算法的圆周率计算程序。这个程序对算法细节进行了大量的优化和微调,代码长度接近800行。其中 pgmp 是利用多线程做的并行楚德诺夫斯基算法。有兴趣的同学可以阅读一下代码了解更多优化细节。

单核C: gmp-chudnovsky
并行C: pgmp-chudnovsky
作者的山寨并行Python版: pi-chudnovsky

下面我们用作者自己的无优化山寨 Python 版本和 gmp 官方给的 C 语言版本的圆周率计算程序进行一下性能测试比较。测试环境是一台8核的Amazon c1.xlarge虚拟机。测试位数是圆周率的前1亿位(100,000,000)。

计算核心(进程数) 时间(C版本) 时间(Python版本)
单核版本 242.3s 397.0s
1 309.4s 419.4s
2 192.3s 226.9s
3 159.9s 168.1s
4 137.2s 135.1s
5 137.5s 116.7s
6 130.8s 105.3s
7 129.7s 102.4s
8 126.4s 99.8s

可见尽管作者的 python 实现的山寨程序尽管算法所用的数据结构优化不够,因此在计算核心较少情况下表现不佳,不过在并行化方面还是做的很好,充分利用到了多核CPU的每一个计算核心。100多行的代码用 python 操作符重载的特性来简化了原本应该很复杂的并行算法的编写。

当然,这些还远远没有到极限。关于圆周率计算的竞赛自古到今一直在进行,在互联网时代越发的激烈。计算竞争的平台也从之前的大型机小型机转到了个人计算机。这里介绍一个著名的圆周率项目叫 y-cruncher

y-cruncher是目前的圆周率计算世界纪录保持者,经作者测试,在刚才的同样测试环境下,计算圆周率1亿位只需要 36.9 秒。y-cruncher 的基本圆周率算法依然还是楚德诺夫斯基算法,但是其代码长度已经到达了17万行!由于优化的足够好,其主要的计算上限还是在内存的大小,y-cruncher的实现主要精力也是在底层内存交换、CPU底层特性的利用以及基础高精度算法的优化上。具体细节大家还是去看他的主页吧。

最后,这篇断断续续写了很久的博客总算完成了。由于篇幅关系,关于博大精深的圆周率以及圆周率文化只是浅尝辄止,希望大家看完以后也能和我一样从此喜欢上深邃的圆周率。

 

 

圆周率之随笔(2)

让我们继续上一章的话题。圆周率是不是祖冲之发明的呢?显然不是。圆周率是一个数学常量,是世界存在以后就亘古不变的一个数字。在人类诞生以前,它就存在了,就一直静静的在那里等着聪明人去发现。所以,很多电视问答题的题目都是错的。那么祖冲之做了什么呢?他是世界上第一位把圆周率算到小数点后第七位的人,同时领先了西方1200多年!

01300000214331124022076846267祖冲之是一个小官二代。祖父是修土木的小官儿,自己一辈子最大也就混了个县令。然而却喜爱数学和工程,在和自己儿子的共同研究之下,研究了一辈子出了一本数学著作,名叫《缀术》。这本书写的是啥,也没人知道了,因为后世就失传了。我们看看唐朝人是怎么评价这本书的。当时唐朝将这本书列为数学教科书,给出的参考学习年数是4年。就是说一般学生要学4年才能把这本书看明白。还有更夸张的,唐朝的算学第一高人李淳风说,“这本书太难了,我们都看不懂!所以大家慢慢也就不去看了”。唐朝的另一个算历博士王孝通说,“有人说这本书很精妙,但是其实里面全都是错的!”。据说宋朝还有人拿着此书练书法。传说是在南宋年间失传的。真相到底如何呢?书已经失传,所以就不得而知了。

《缀术》虽然已经失传了,那我们从何而知祖冲之的圆周率的呢?这个描述来自于《隋书》的《律历志》。作者就是唐朝的算学大师李淳风。原文如下:

古之九数,圆周率三,圆径率一,其术疏舛。自刘歆、张衡、刘徽、王蕃、皮延宗之徒,各设新率,未臻折衷。宋末,南徐州从事史祖冲之,更开密法,以圆径一亿为一丈,圆周盈数三丈一尺四寸一分五厘九毫二秒七忽,朒数三丈一尺四寸一分五厘九毫二秒六忽,正数在盈朒二限之间。密率,圆径一百一十三,圆周三百五十五。约率,圆径七,周二十二。又设开差冪,开差立,兼以正圆参之。指要精密,算氏之最者也。所著之书,名为《缀术》,学官莫能究其深奥,是故废而不理。

翻译成白话文就是说,南朝的徐州从事史,祖冲之,非常厉害,把圆直径一丈分为一亿份,周长在三丈一尺四寸一分五厘九毫二秒六忽到七忽之间 (3.1415926 < π < 3.1415927),等等。

这里略让人产生疑惑的一点是,除了唐朝李淳风在《隋书》里提到这么一句之外,翻遍整个中国古籍,都找不到其它任何关于祖冲之圆周率的描述。甚至在祖冲之之后1000多年的所有古代数学家们,也都全部不知道祖冲之能有如此之高的圆周率精度。看下表:

年代  朝代 数学家 圆周率 精度
公元前100 战国 《周髀算经》 3 0
公元前50 刘歆 3.1547 1
公元78 张衡 730/232=3.1467
92/29=3.1724
sqrt(10)=3.162
1-2
公元133 蔡邕 3.125 1
公元228 三国 王蕃 142/45=3.155 1
公元300 魏晋 刘徽 157/50=3.14
3927/1250=3.1416
2-3
公元370 何承天 3.1429 3.1432 2
公元429 南朝 祖冲之 3.1415926~3.1415927
22/7 355/113
7
公元602 李淳风 22/7=3.1429 2
公元1208 秦九韶 sqrt(10)=3.162 1
陈荩谟 3.1525 1
公元1536 朱载堉 sqrt(2)/0.45=3.1427 2
邢云路 3.1213 3.126 1
公元1533 程大位 25/8=3.125 1
公元1611 方以智 52/17=3.0588 0
公元1714 王元启 sqrt(10)=3.162 1
顾长发 3.125 1
公元1728 钱大昕 3.16 1

祖冲之在中国古代2000多年的圆周率史上,如皓日当空那般耀眼,不仅仅是前无古人,而且当真做到了后无来者(西方数学家算出来的不算)。

为何唐朝人写的《隋书》中,明明写清楚了祖冲之的高精度的圆周率,但是后世1000多年,这些个算学大师们一个一个都闭口不提,而且反而用精度更差的圆周率呢?《梦溪笔谈》的作者沈括也因为使用了很差的圆周率,导致做出来的工具很快就坏了。可见后世人完全不知道祖冲之有了如此了不起的研究成果。

本人有下面几个猜测:ri08_01735_0164_p0030

  1. 可能古代算学家根本不重视圆周率。算个大概就不管了。古人的成果也压根没有人去看过。古代的科学基本没有传承。只有现代中国人才从古老典籍里挖掘出一点点古人的辉煌聊以慰藉。如果是这样,真的是很可悲的事。
  2. 《隋书》可能是造假的。唐朝人写的隋书原本早就在历史的尘嚣中烟飞云散了。宋朝刻印的《隋书》留下的又不全,只有残破的几章。目前能找到的最早版本有记载祖冲之圆周率的《隋书》,是明末的毛氏汲古阁本。是不是后人改的?如果是的话,恐怕更加的可怕。

让我们继续推测,祖冲之是怎么计算圆周率的?这已经成为了一个千古之迷。很多人猜测也是和阿基米德一样的割圆法。是的,这是最合理的猜测了,因为在公元400年的古代,割圆法只需要开方根和勾股定理的知识就可以计算。

借助现代的计算机,我们很容易可以模拟当年的先辈是怎么用外接割圆和内接割圆计算圆周率的:

s = 0.5
for i in xrange(n):
    s = (s*s+(1-(1-s*s)**0.5)**2)**0.5/2
t = (s*2/3)**0.5
for i in xrange(n):
    t = (t**2-((t**2+1)**0.5-1)**2)/t/2
r = 6*2**n
print r, s*r, t*r

这段简单的小程序不另作说明了,通过计算内接和外接割圆边长,我们分分秒就能计算出每一次切割能达到的圆周率精度,和计算过程中所需要保证的小数点后的精度个数:

n边形 到达小数点后精度 过程需要保持精度 结果
6 0 2 3.0<pi<3.4
12 0 3 3.11<pi<3.22
24 1 4 3.132<pi<3.158
48 1 5 3.1392<pi<3.1462
96 2 6 3.14102<pi<3.14268
192 3 7 3.141454<pi<3.141873
384 3 8 3.1415576<pi<3.1416626
768 3 9 3.14158389<pi<3.14161021
1536 5 10 3.141590466<pi<3.141597034
3072 6 12 3.1415921056<pi<3.1415937485
6144 6 12 3.14159251670<pi<3.14159292743
12288 6 13 3.141592619364<pi<3.141592722051
24576 7 14 3.1415926450333<pi<3.1415926707020

可以看到,阿基米德在公元前250年,通过切割96边形,那么他需要在计算过程中保持6位的小数点后精度,才可以计算出 3.14102<pi<3.14268 的结果。祖冲之在公元400多年,如果需要保留小数点后7位精度的圆周率精度,需要切割到24576边形,同时计算过程中要保持14位的小数点后精度!

下面再具体列举一下我们的革命好同志祖冲之计算圆周率之艰辛:

  1. 进行12次迭代(也有猜测是从4边形开始隔,这样只要11次迭代),每次迭代都需要计算3次平方根,若干次乘法除法,另外计算过程要保证14位的小数点后精度。
  2. 没有阿拉伯数字。无论计算过程用的是一二三还是壹贰叁,都麻烦的不得了。
  3. 没有算盘可以打。用的是算筹。把很多小树枝放在地上,拨过来拨过去的,拨错一个就完蛋。
  4. 没法合作完成。算法不能并行操作,没法团队完成,必须一个人自己默默的算。
  5. 至少得验算2-3次吧?历史上算错圆周率贻笑大方的人很多。能一次性算对的很难很难。
  6. 用的是毛笔。没有铅笔钢笔圆珠笔。不管是记录中间计算过程还是结果,都得用掉成百上千纸张。

好同志祖冲之如果真的算到第7位,目测按照当时的条件,怎么也得算个三年五载的。为了一个在当时没什么太大意义的高精度圆周率,耗费如此多的光阴,让人觉得有点匪夷所思。

如果是后人伪造《隋书》编造的 3.1415926 ,那么为何要伪造呢?这里也有一些背景知识需要阐述。

c2cec3fdfc039245010022028794a4c27d1e2502明末清初的年代,一个很著名的学派,是西学东源学派。随着越来越多的西方传教士不远万里来到中国,为了传教学会了中文,然后把一些西方科学著作翻译成了中文,其中最著名的是西儒利玛窦,让当时的学者们认识到了中西的差距已经大到可怕了。西学把日历能推算到秒,可以算出地球是个球体,知道地球的大概半径,能够制造精确的钟表。工艺技术的差别,天文地理的差别,基础数理的差别,让当时的学者们无法面对。

但是总是需要去面对的。面对的方式不是说努力学习西文典籍迎头赶上,而是很猥琐的鸵鸟策略。西学东源的主要思想是,你们西学研究的那些东西,我们老祖宗典籍里全都有了,不足为奇!这些学者依靠各种牵强附会的考证,想要证明西算(数学)和西历(天文学)都是剽窃的我们的老祖宗!这些可耻和可悲的考证成果就不一一说了。这个学派一直到清朝继续发扬光大,清朝皇帝康熙和乾隆一边焚书坑儒,编纂删改古籍,最后出了删改后的《四库全书》,一边组织文人搞西学东源运动,接见学派首脑,让这个学派继续发扬光大。恐怕直到北洋舰队被日本打的个落花流水,以及八国联军打进了北京城,才让当时的中国人真的觉醒,才让他们能够发现,中西方的差距已经实在是太大了,赶紧老老实实的“师夷长技以制夷”吧!(虽然最后也没能制得了夷)

在明末刻印的《隋书》版本,是不是有经过算家的修改,也不得而知了。真相早已掩埋在历史长河中。本文的目的也不是要考证祖冲之是不是算了圆周率,只是把一些证据陈列出来,信者自信疑者自疑吧。即使祖冲之确实算到了第7位,但是后世数学家们的表现也让人不禁感慨,古代对算学这些奇技淫巧的忽视。伟大的人物昙花一现,只有后世西学传入之后,才让人真正了解他们的伟大的地方。

随着电子计算机的诞生这几十年,圆周率的计算长度已经开始几百倍几千倍的增长了,每天都会有新纪录诞生。下一章我们回到正题,谈谈圆周率的常用算法。

 

 

圆周率之随笔(1)

有人说,爱像圆周率,无限不循。那么,读懂了圆周率,你就读懂了爱:简单而又深邃。(这其实是标题党,第一句总要吸引人吧)

圆周率的定义很简单:一个圆的周长除以直径,就是圆周率。然而,圆周率已在人类世界存在了千年,至今还有很多关于它的谜团存在。

Pi-unrolled-720从人类文明有记载开始,远古文明就开始用 π = 3来开始大致计算圆的周长了。如果你不幸被穿越到异世界,怎么确定这个世界的第一常量,圆周率呢?很简单,拿一根藤条,先在地上靠一个支点用固定半径画出一个圆,然后用同一根藤条对着周长比划一下,看看是不是 3 根藤条的长度然后多一点点。如果是,那么恭喜你,你至少在一个常量一致的异世界里!

622762d0f703918f1e9863ba513d269758ee3d6d55fbe449圆周率比3多一点点。多这么的一点点到底是多少?远古文明从人类诞生就开始了探索。公元前1900-前1600的古巴比伦人,在一块后来被发现的破泥板上刻下了圆周率大概是 25/8 = 3.125。这大概是有文字记载的最早的圆周率。与此年代悠久不相上下的是古埃及人,曾经在一张后世被拍卖成天价的破草纸上写下的公式:(16/9)= 3.1605。这个时代我们华夏的古老祖先,还在黄河流域治水灾,没有留下任何文字的烙印,只有一些口口相传的神话故事。中华文明直到大概是战国时期(年代不明,可能是公元前300左右吧)的《周髀算经》里,第一次提到过圆周率:“径一周三”。

Domenico-Fetti_Archimedes_1620

之后的某一天,地球上一个在树下被苹果或者梨砸中的天才忽然想,为什么不能用圆内接的多边形来计算圆的周长呢?最开始这位天才挑中了六边形,因为六边形是边长和半径相等的特殊多边形。然后通过分割,把六边形切成十二边形,然后是24、48、96边形。随着边越多,多边形的周长也就和圆的周长越来越接近了。这位天才就是阿基米德。这个计算方法就是割圆法。

通过计算了96边形的周长,阿基米德给出了 π 的范围:223/71 < π < 22/7 (3.1408 < π < 3.1429)。在公元前250年时代,这是一个极其惊人的成就!可以说,亚里士多德、阿基米德等先贤辈出的古希腊科学家,在距离现今的2000多年前(比我们想象的早很多的时候)就奠定了西方的科学基础,中西差距就此拉开,然后越拉越远。从那个远古时代,西方科学家的地位就变得无比崇高,科学的火种就此被点燃了,然后越燃越旺,最后迸发出难以想象的活力。

而相比之下,中国同时期战国时代的先贤的追求就很不一样。诸子百家都是周游列国为君王出谋划策,以实现他们的政治理想,不管他们是不是什么思想家、教育家,他们首先都是政治家。了却君王天下事,赢得生前身后名,是他们最大的心愿。孔子对国王说的是:“君君,臣臣,父父,子子”;相比之下,同一时期,地球另一边的欧几里得在给国王讲几何学,同时嘲讽国王太笨学的太慢:“在几何学里,没有专为国王铺设的大道”。其追求不同如云壤之别。

除了计算圆周率,牛人阿基米德还在各个领域都有跨时代的突破,以及留给后世诗人们传唱的数不尽的传说故事。其中最著名的一个,就是他在浴缸里洗澡,忽然想出了浮力定律,激动的冲出浴池去裸奔的故事。

300px-Archimedes_Heat_Ray_conceptual_diagram.svg做为皇亲国戚阿基米德,据说为人疯疯癫癫,可以从裸奔一事里窥见一斑。阿基米德的另一个著名传说故事是,为了防守祖国叙拉古城邦,让士兵用了很多很多面镜子,将罗马的海上无敌舰队全部点燃,导致罗马落荒而逃。这个故事是在400多年后(公元2世纪)出现在罗马诗人卢坎的诗集里。然而充满质疑精神的西方,后世有很多人做过各种实验怀疑这个故事的真实性。在1973年,希腊科学家组织过一次实验想要重复阿基米德的故事,之后在2005年麻省理工学院的学生,甚至在2010年12月,美国总统奥巴马都赞助组织过几次类似实验。所有实验结论证明,这个传说故事很可能是假的。古代的铜镜由于反射度不好,外加隔着这么远的大海,除非船只都在那里一动不动,否则很难被镜子点燃。更可能的是阿基米德用他发明的投掷车,远程投掷燃烧弹的办法打败了罗马舰队。

这里不禁感慨一下中西文化的差异之大。西方文化里,质疑先贤是时代进步的象征,阿基米德的故事都可以拿来推翻;而东方文化里,在孔子儒家的法先祖的教育下,崇拜先贤才是思想上正确的,质疑古人被认为是不敬不尊重。

可悲可叹这位伟大的天才人物,最后死在了一个率先攻入叙拉古的罗马士兵手中。阿基米德在看到很多罗马士兵冲进房间之后,依然疯疯癫癫的说:别动我的圆!随后,阿基米德死在了听不懂希腊语而恼羞成怒的罗马士兵手里(或者他听懂了希腊语,但觉得阿基米德是在表现科学家的高傲)。可见阿基米德临死前可能还在研究圆的问题!

一直到公元1400年以后,西方进入了文艺复兴时期,科学技术忽然开始了蓬勃发展,随之而来的是 π 精度日新月异的提高。奥地利人Christoph Grienberger在1630年用多边形法将 π 计算到了38位,这也是阿基米德发明的多边形法的终极记录。随后便进入了级数展开的圆周率新时代。

在中国,街头采访问祖冲之是做什么的?大多数人都会说,祖冲之发明了圆周率。还有某一次电视知识竞赛里,电视节目也给了类似问题:谁发明了圆周率?答案是祖冲之。这个问题应该怎么答?祖冲之和圆周率有什么关系?咱们下一篇再见!

 

谈谈Zopfli

image01

最近 Google 出了一个新的开源项目 Zopfli。Zopfli是什么呢?简单说是一个 Deflate 压缩算法的另一种实现。推出之后国内国外媒体纷纷报道转载。昨天看到国内媒体的报道 (搜狐IT)中说道:“据悉,Zopfli的压缩率比现有的Zlib高3-8倍。”当时看到了就吓了一大跳,3-8倍这是要逆天啊!赶紧去Zopfli主页看一眼,原来只是3%-8%的提升。国内IT编辑真是吓死人不偿命。

Zopfli到底表现怎么样呢?先看看 Zopfli 自己提供的数据:

测试集 样本大小 gzip -9 7-zip kzip Zopfli
Alexa-top-10k 693,108,837 128,498,665 125,599,259 125,163,521 123,755,118
Calgary 3,141,622 1,017,624 980,674 978,993 974,579
Canterbury 2,818,976 730,732 675,163 674,321 669,933
enwik8 100,000,000 36,445,248 35,102,976 35,025,767 34,955,756

大概是比标准的gzip -9要小 3.7%-8.3% 的样子。但是压缩所需要的时间非常的惊人:

压缩算法 压缩时间
gzip -9 5.60 s
7-zip -mm=Defalte -mx=9 128 s
kzip 336 s
Zopfli 454 s

可以说如果用 gzip 压缩要喝口水的时间的话,用 Zopfli 就可以出去吃个便饭了。

我们再看看这些算法的出生年份:Deflate最早诞生于1994年。7-zip最早诞生于1999年,也是一个开源项目,虽然7-zip自己默认用的格式是LZMA和LZMA2算法,但是也支持对Deflate算法的更好的压缩。KZip 是 Pngout 作者写的一个小工具,诞生于2006年,是一个把现有Zip再压缩的工具。

这里不得不提几个Zopfli文档里没有提到的其它 Deflate 兼容的压缩工具:DeflOpt,是2007年诞生的一个Zip再压缩工具。AdvZip,是一个利用7-zip Deflate的一个Zip再压缩工具。

可以从表中看出,Zopfli比1994年原版Deflate确实有3%-8%的提升,但是对比比较新的Deflate实现,提升实在是很有限,比如比KZip只提升了大概1%不到。压缩时间对比KZip也增加了相当多。

发布没多久,就有人发现,利用 KZip+DeflOpt,压缩的结果可以比Zopfli压缩的更小,甚至Zopfli的几个样本的压缩结果,依然可以用 DeflOpt 再压缩一点(1%左右)。所以说Zopfli其实并不是他在文档里所说的,“所有已知Deflate算法实现里压缩比最高的”。

当然,Zopfli的意义在于它是开源算法,而KZip+DeflOpt这俩都不开源。关于KZip,DeflOpt是如何实现的,之前还有很多人在猜测。Zopfli的出现给大家提供了很好的解答。

另一点值得探讨的是,在今天研究Deflate算法更好实现是否还有价值。Zlib的作者本人Mark Adler在听说Zopfli之后说:“这很酷,不过看上去是一个付出了很多努力,但只取得了很小提升的一个糟糕结果。也许到了给HTTP的accept-encoding加上更好算法的时候了”。是的,诞生于1994年的Deflate确实太老了。现有的解压更快、压缩比更高的开源算法有很多很多,比如 bzip2, LZMA/XZ等等。LZMA/XZ在解压速度上,以及压缩比上,都完胜Deflate。Deflate对于很多 UTF-8 3 bytes的网页压缩效果也很不理想。直到今天,支持bzip2, LZMA的浏览器还寥寥无几。也许更有前途值得关注的是HTTP/2.0协议(一部分基于Google的SPDY协议)。

重新回到Deflate算法。为何Deflate有如此多的压缩实现呢?我们得详细的看看Deflate算法的具体内容。Deflate算法其实就是LZ77算法加Huffman算法。先经过LZ77的字典找重,然后用Huffman树进行降比特。不同的Deflate压缩的实现其关键在于LZ77的搜索重复单词,以及选择分块来进行Huffman。早期算法例如7-zip对分块都不是很重视,更多的考虑的是LZ77算法的优化。而KZip另辟蹊径对分块也进行了优化,使得最终比特流长度变得更短。

LZ77算法优化是一个有向图上的最短路径的搜索问题。对于字节流每个字节建立单边节点,重复单词序列建立更短的边,形成一个有向图。LZ77算法的目的就是如何找到有向图从起点到终点的最短路径。对于边的长度确定的情况下,用动态规划找最优解是很简单的。然而LZ77在搜索时,如果要考虑下一步Huffman过滤之后的长度,则边长度就是不固定的。Zopfli采用迭代的办法,先用一次贪心法拿到第一个次好的结果,然后通过使用结果字节的熵值(就是出现频率的N-log2n)来给出下一次迭代的每条边的边长,也就是每个字母的比特长度。通过反复迭代,来逐步逼近最好结果。当然,理论上也可能跌入一个次好结果的低谷。

Huffman_tree_2.svg然后是最关键的分块问题。分块为何会影响Huffman的压缩结果?其主要问题还是因为动态Huffman算法的随机性。如果都采用对字节使用频率统计完毕之后的静态Huffman,那么不管如何分块都不会对结果有影响,反而因为分块产生额外的比特,使得结果变大。动态Huffman的树是在字节流的处理过程中动态创建的,因此其字节流的开始片段不规律性往往使得结果不优。如何分块才能更优?因为随机性太大,也很难进行判断。所以KZip和Zopfli采用的都是尽可能的穷举。KZip就号称用了“重复单词的穷举(LZ77)”和“更高效率的分块(Huffman)”来实现,可以增加分块的个数进行进一步优化。而Zopfli,是不停的在一个最大块内找9个点,穷举判断哪一个最优,然后进行反复切割的办法。

Deflate压缩算法有没有最优解?这其实是个NP问题。对于短短32字节大小的数据,都有2**31 = 2147483648 种不同的分块方案。当然绝大多数分块都没有意义而产生更差结果。分块之后依然还有需要进行迭代的边长会变化的最短路径问题需要解决。

最后,尽管 Zopfli 的结果不是很令人满意,不过确实给众多不开源的 Deflate 压缩工具树立了标杆。那些想靠着 Deflate 算法做收费 PNG 压缩的软件可以洗洗睡了。