祖冲之是怎么计算圆周率的？

简介：如果我们按《九章算术》刘徽注文中的一个思路动手去尝试一下，就会惊奇地发现：筹算出8位圆周率并不神奇。如果祖冲之使用这个方法，那么他2天就能可靠地将圆周伟哥有效果吗率重算一遍，计算量及复杂程度都远低于大家的预期。

1500多年前祖冲之精确计算出圆周率在3.1415926～3.1415927之间，但由于祖冲之所著《缀术》失传，后人再也不知道他是怎么计算的，因此，既有人神化祖冲之，也有西方学者怀疑祖冲之圆周率是中国后人为美化自己的历史而造假。我在带小孩认识圆周率的过程中，偶然发现用算筹可以快速并可靠地进行高精度开方运算，翻阅《九章算术》后发现刘徽写的注文中也描述了计算开方的近似值的相似方法，可能祖冲之计算圆周率就这么简单，神化和造假都是多余。

用几何法计算圆周率的先驱是古希腊阿基米德，他用圆内接正多边形的面积及圆外切正多边形的面积双向逼近圆的面积。中国古代数学家刘徽用圆内接正多边形的面积及其外加多个矩形的面积双向逼近圆的面积。阿基米德和刘徽的思路都直观易懂，关键的运算也一样，就是如何快速并可靠地计算出精确的开方值。

祖冲之年代尚未出现算盘，运算仅靠小木棍（算筹），靠它算出8位圆周率，还真有点天方夜谭的感觉。虽然祖冲之年代很久之前就已用十进制计数，但小数还没出现，那时如果不能整除，是用分数运算的。祖冲之圆周率记为三丈一尺四寸一分五厘九毫二秒六忽，相当于今天的3.1415926。

《九章算术》注文中刘徽对割圆术及其几何计算进行了详细的描述，方田章节有约分运算的详细描述，少广章节有开方运算的详细描述，方程章节有负数及其运算的详细描述，等等，这些都可以是祖冲之计算圆周率的数学基础。

刘徽从圆的内接正6边形开始割圆。如上图，设圆的半径为1，圆内接正n边形边长 X_{n} = AB ，圆内接正2n边形边长 X_{2n} = AC 。对三角形BDO使用勾股定理，求得股长 OD = \sqrt{1 - ~X_{n}^{2}/4} ，小勾 CD=1-OD ，再对三角形CDA使用勾股定理，求得 X_{2n}^{2} =曹操盗墓 {AC}^{2} = {CD}^{2} + X_{n}^{2}/4 。因为 X_{6}=1 ，所以递推可计算出 X_{12}^{2}、X_{24}^{2}、X_{48}^{2}\ldots\ldots 每次割圆计算出 X_{2n}^{2} ，除了简单的加减运算之外，复杂的运算包括一次开方运算和一次平方运算（小勾CD的平方）。

刘徽割圆通过两次勾股定理计算出 X_{2n}^{2} 的几何算法，可以被简化如下： {AC}^{2} = {CD}^{2} + {AD}^{2} = \left( {~1 - OD~} \right)^{2} + {AD}^{2} = 2 - 2OD = 2 - 2\sqrt{1 - ~{AB}_{n}^{2海的儿子}/4} ，即 X_{2n}^{2} = 2 - \sqrt{4 - ~X_{n}^{2}} ，利用这个递推公式，每次割圆计算与刘徽的计算相比，可以减少一次平方运算。

刘徽用上图双向逼近圆的面积，扇形AOB的面积大于四边形AOBC的面积、小于五边形AOBEF的面积，而且矩形ABEF的面积为三角形ABC面积的两倍。假设 S_{n} 表示圆内接正n边形的面积， S 表示圆的面积，则有不等式 S_{2n} < S < S_{n} + 2*\left( {S_{2n} - S_{n}} \right) 。《九章算术》中说明了圆周率乘以半径为半周长、圆面积等于半周长乘以半径，因此上述不等式等价于 X_{n}*\frac{n}{2} < \pi < X_{n/2}*\frac{n}{4} + 2*\left( {X_{n}*\frac{n}{2} - X_{n/2}*\frac{n}{4}} \right) ，仅需计算出圆内接正n/2边形的边长 X_{n/2} 及圆内接正n边形的边长 X_{n} ，就可以计算出圆周率π的范围。刘徽在《九章算术》注文中使用了上述方法，但从上述计算过程可以看到，这些计算都绕不过令人头疼的开方运算。

开方对古代的人们而言是比较复杂的，否则阿基米德在公元前200多年就可将圆周率精确计算到小数点后10位或更高位了。古代中国和古代欧洲，使用的数都是整数和分数，用几何法计算圆周率无法绕过开方运算。这些数开方之后都是无理数，用整数和分数的加乘运算来表达无理数只能是近似值，表现出来的问题就是找不到很好的方法来精确计算开方值。精确计算圆周率需要高精度开方，这对古人是巨大的挑战，哪怕采用现代计算工具，实现一个高精度开方（例如精度为100位）也不是件易事。因为我家12位老外看点的计算器开方算出的结果精度不够，导致计算出的圆周率不在3.1415926～3.1415927之间而引来小孩怀疑的目光，我们才偶然遇到祖冲之计算圆周率这个问题并去研究它的。

《九章算术》对开方有详细的描述，今天部分教科书上手工开方的方法江水生就是来源于《九章算术》，然而这种开方方法，为了规避解二次方程需要采用试乘，多次试乘运算不但增加了运算量而且还增加了出错的可能。刘徽在《九章算术》注文中，仅将圆周率可靠地计算到3.14（史称徽率，另一有争议之说刘徽计算出圆周率为3.1416，但这个数值的计算方法不是严谨的），应该就是饱受筹算多位数开方之苦。

《九章算术》开方术中，刘徽写了一段注文“术或有以借算加定法而命分者，虽粗相近，不可用也。凡开积为方，方之自乘当还复有积分。令不加借算而命分，则常微少；其加借算而命分，则又微多。其数不可得而定。故惟以面命之，为不失耳。譬犹以三除十，以其余为三分之一，而复其数可以举。不以面命之，加定法如前，求其微数。微数无名者以为分子，其一退以十为母，其再退以百为母。退之弥下，其分弥细，则朱幂虽有所弃之数，不足言之也”。这段注文的大概意思是：对A开方开不尽，如果 A = B^{2} + ~\Delta ，那么 B + \frac{\Delta}{2B + 1}~ < 英语口头语~\sqrt{A}~ < ~B + \frac{\Delta}{2B} ，并且 \sqrt{A} \approx B + \frac{\lambda_{1}}{10} + \frac{\lambda_{2}}{100} + \cdot \cdot \cdot \cdot \cdot \cdot

刘徽计算开方近似值的不等式可以比较容易地证明（网上搜索说阿基米德当年也用了这个不等式），刘徽在开方近似值的表达式中引入了微数这一概念。《九章算术》常见的单位为丈尺寸分厘毫秒忽，这些单位用完之后，就会出现令现代人抓狂的表述，例如刘徽割调制圆计算圆内理性情绪疗法接正12边形的股长为“九寸六分五厘九毫二秒五忽五分忽之四”，粗看我们容易理解为股长0.96592554，正确的股长为0.9659258，单位“忽”用完之后出现了“五分忽之四”，即0.8忽。估计是因为计量单位不够用很不方便，所以刘徽扩展出了微数这一概念。因为每个微数没有具体的名称，所以不同位置上的微数的含义是不同的，但它的运算规则与丈尺寸分厘毫秒忽类同。刘徽的微数虽然采用十进制分数表示，但其实质完全等同于现在的十进制小数。刘徽这段注文的目的是对开方术做补充，而编辑器漏洞开方本来就是件令人头疼的事情，更何况开不尽的情况，所以有兴趣去思考它的人应该是少之又少。刘徽本人应该也没有认真仔细地研究过微数表达式与开方近似值范围相结合的妙用，否则刘徽当年就可以严格地将圆周率计算到小数点后四位或七位或更高了。

祖冲之熟读过《九章算术》，能准确理解刘徽注文的含义，采用刘徽计算圆周率的方法进一步提升圆周率的精度，唯一的障碍就是高精度开方。我推测祖冲之尝试使用了一下刘徽开方的近似计算思路，利用刘徽开方近似值的上下限范围计算出第一个微数 \frac{\lambda_{1}}{10} \approx \frac{A - B^{2}}{2B} ，然后验黄金棋牌算开方近似值的平方是否与A相等或与A还有多少差距， B_{1}^{2} = \left( {B + \fr网上买家具ac{\lambda_{1}}{10}} \right)^{2} = B^{2} + 2*\frac{\lambda_{1}}{10}*B + \left( \frac{\lambda_{1}}{10} \right)^{2} 。如果加了一个微数之后开方近似值的平方与A还有较大差距，很容易想到用同样的方法计算出第2个微数 \frac{\lambda_{2}}{100} \approx \frac{A - B_{1}^{2}}{2B_{1}} ，再次验算开方近似值的平方是否与A相等或者与A还有多少差距， B_{2}^{2} = \left( {B_{1} + \frac{\lambda_{2}}{100}} \right)^{2} = B_{1}^{2} + 2*\frac{\lambda_{2}}{100}*B_{1} + \left( \frac{\lambda_{2}}{100} \right)^{2} 。利用同样的方法，祖冲之想将开方的精度提升到多少个微数都可以，开方的精度不受限！我们可对上述开方每增加一个微数精度的运算步骤，总结如下：

第一步： 计算第k个微数 \frac{\lambda_{k}}{10^{k}} \approx \frac{A - B_{k - 1}^{2}}{2B_{k - 1}} ，由于 \lambda_{k}~ 仅是1位数，因此分母可四舍五入取整，分子四舍五入保留k个微数，计算结果四舍五入取1位，易于计算。

第二步：计算B_{k~}^{2} = \left( {B_{k - 1} + \frac{\lambda_{k}}{10^{k}}} \right)^{2} = B_{k - 1}^{2} + 2*\frac{\lambda_{k}}{10^{k}}*B_{k - 1} + \left( \fra红磷和白磷c{\lambda_{k}}{10^{k}} \right)^{2}，仅需将演艺厅已经算出的 B_{k - 1~}^{2}加上 2*\frac{\lambda_{k}}{10^{k}}*B_{k - 1} 及 \left( \frac{\lambda_{k}}{10^{k}} \right)^{2} 两项，后两项都易于计算。

第三步：多位数的加减运算。

上面三个步骤，使用算筹运算都很方便。因此，利用上述方法，可以将复杂的高精度开方运算转化为简单的加乘运算。

当被开方数A在区间 [ 3，25 ）的时候，上述开方仅会涉及到1～2位数除以1位数的除法（四舍五入取1位结果，容易心算出来）、多位数乘以1位数的乘法（乘除10因为很简单不计入2位数乘除法）及多位数加减法。实际运算过程中由于使用四舍五入，如果发现 λ_{k} 大于或等于6，则可在已算出的第（k－1）个微数的基础上，在不增加微数个数的前提下，按计算微数的方法将第（k－1）个微数增加或减少1，以便减小 \left( {A - B_{k - 1}^{2}} \right) 差值战胜心魔，避免计算 B_{k}^{2} 时需要用2位数乘以多位数的乘法。

为了简化计算，当被开方数A不在区间 [ 3，25 ）的时候，可先对A做预处理如下：区间（1，3）内的数值可以先乘以4之后再开方；区间（25，100）内的数值可以先除以4之后再开方；区间（0，1）内的数值可以先乘以10的偶次方之后再开方；100以上的数值可以先除以10的偶次方之后再开方。

下面举个例子说明这种开方的过程，例如对 A = 4 + \frac{6}{10} + \frac{7}{100} + \frac{8}{1000} + \frac{9}{10000} 进行开方运算，开方值B的微数逐个增加如下：

2的平方最接近A，所以取 B_{0} = 2

上述计算表面看起来有点吓人，其实用算筹计算时却很简单，分数中零的个数主要用来对准算筹的位置，例如8/1000表示8在整数位之后的第3个算筹位上。刘徽所称的微数，实质就是十进小军制小数，为了便于描述，本文后面的段落不法医学研究生再使用微数而直接更换成小数。将上述“计算第3个微数”的步骤改成“计算第3个小数位”如下，这样看就不会觉得它复杂了。

\lambda_{3} \approx \left( {A - B}_{2}^{2} \right)/2B_{2} \approx \left( {4.679 - 4.666} \right)/4 = 0.013/4 \approx 0.003

B_{3} = B_{2} + \lambda_{3} = 2.16 + 0.003 = 2.163

B_{3}^{2} = \left( {B_{2} + \lambda_{3}} \right)^{2} = B_{2}^{2} + 2\lambda_{3}B_{2} + \lambda_{3}^{2} = 4.6656 + 2*0.003*2.16 + 0.003^{2} = 4.678569

筹算加减法无需像我们常用阿拉伯数字竖式运算那样还要重复抄写一些数字、对齐小数位也很方便，因此熟练筹算之后，筹算速度有可能达到甚至超过我们用阿拉伯数字竖式运算。筹算多位数乘以一位数的乘法也比较容易。我试过筹算一次16位有效精度的开方运算，大约需要1小时。由于这种筹算开方的方法仅包含几种简单的运算，因此比较快、也不易出错。2016年7月我拍了个简单的视频，在爱奇艺网站 www.iqiyi/ 搜索“儿子和我一起用小棍子精确计算圆周率”，便能看到7岁的孩子就可以胜任16位高精度筹算开方。

如果祖冲之对刘徽割圆术的几何计算方法进行过简化，得到前述圆内接正多边形边长的递推公式 X_{2n}^{2} = 2 - \sqrt{4 - ~X_{n}^{2}} ，那么经过12次16位精度的开方可得到下表：

利用上表最后两行的结果，再次开方可计算出 X_{12288} = \sqrt{2 - ~1.999999738544777} = 0.000511326923797 和 X_{24576} = \sqrt{2 - ~1.999999934636193} = 0.000255贪婪算法663464343 ，总计14次开方。将 n=24576及 X_{n} 、 X_{n/2} 代入前述刘徽计算圆周率的公式 X_{n}*\frac{n}{2} < \pi < X_{n/2}*\frac{n}{4} + 2*\left( {X_{n}*\frac{n}{2} - X_{n/2}*\frac{n}{4}} \right) ，得到 3.141592649846784 < \pi < 3.141592679884800 ，因而可确定圆周率介于3什么是皮肤癌.1415926～3.1415927之间。其中开方要精确到小数点后15位或更高（读者可以试一下，将上述表格最后一行四舍五入减少1位，那么π的下限就超过3.1415927）。

将圆周率界定在3.1415926～3.1415927之间总计需要14次16位或以上的高精度开方运算。祖冲之筹算的速度一定比我快，因此他重算一次圆周率世界尽头与冷酷仙境仅需要2天时间，根本无需常年累月地计算。如果祖冲之未对刘徽割圆的几何计算进行简化，完全按照刘徽的割圆计算，祖冲之需要做的开方运算也是14次，但要多做12次平方运算，因此他重算一次圆周率的时间也不会超过3天时间。

汉字“缀”有“补上、修饰”的意思，前述逐位增加开方值的小数位数，在增加小数位数时正负双向逼近最优开方值，因此，这个开方的方法很符合“缀”字“补上、修饰”的意思，因此我简称前述开方的方法为“缀开方”。

祖冲之生活的年代，阿拉伯数字还没被传入中国、算盘也没被发明出来，仅依靠《九章算术》的开方术将圆周率界定在3.1415926～3.1415927之间是无法想象的，也正因为此，国人神化祖冲之、西方怀疑祖冲之的成就。《缀术》失传之后“祖冲之如何计算圆周率”成了历史悬案，上至数学家华罗庚下至广大普通的数学爱好者都思考过祖冲之计算圆周率的问题，《科技导报》2008年第5期还将“祖冲之究竟是怎样计算出圆周率π值的？”列为中国公众关注的18个未解科学难题之一。通过上文，读者看到筹算圆周率可以就这么简单。刘徽已经在《九章算术》的注文中示范了如何通过迭代的方法计算圆周率、无意中还指了一下高精度开方的方向，因此，在颂扬祖冲之圆周率的时候，不能忽略刘徽的巨大贡献。

上文假设祖冲之受到了刘徽开方近似计算思路的启发，但也不能排除祖冲之自己独创了上述的缀开方。我推测祖冲之在《缀术》中，可能详细地记载了缀开方。缀开方完全可以替换《九章算术》中的开方术，因为缀开方比《九章算术》的开方术更科学、更简洁、更容易计算、更不易出错。详细的筹算开方过程，可参看《数学通报》2018年第6期的“小孩和我偶遇组冲之计算圆周率之谜”。

上述内容是我们在发现高精度手工开方的方法之后，反推祖冲之是否可以利用《九章算术》已有的知识计算出他的圆周率，所以可能会有其他未考虑到的问题。

Zou@南京