祖冲之(429-500)确定圆周率至3.1415926~3.1415927，且得到密率355/113，乃是我国古代数学之骄傲，但具体的计算方法已不得而知；不过，其沿用和发展了刘徽（225 - 295）割圆术则无疑问，基于“以圆径一亿为一丈”或可复原其过程。 

1   刘徽的圆周率

《九章算术注·卷一》有刘徽割圆法论证“圆面积为半周与半径之积”，且记录了计算圆周率的具体过程。利用勾股定理逐次求内接正多变形边长的平方，再求边长计算多边形面积，并估算误差。抄录部分内容，并略作说明。

割六觚以为十二觚术曰：置圆径二尺，半之为一尺，即圆里觚之面也。令半径一尺为弦，半面五寸为句，为之求股。以句幂二十五寸减弦幂，余七十五寸， 开方除之，下至秒、忽。又一退法，求其微数。微数无名，知以为分子，以十为分母，约作五分忽之二。故得股八寸六分六厘二秒五忽五分忽之二①。以减半径，余一寸三分三厘九毫七秒四忽五分忽之三，谓之小句。觚之半面又谓之小股。为之求弦。其幂二千六百七十九亿四千九百一十九万三千四百四十五忽②，余分弃之。开方除之，即十二觚之一面也。

注：圆半径为1 尺，六边形边长1 尺。数字皆有长度单位，忽为最小单位，其后可用分数估计。① 0.8660254尺就是二分之根号三，毫没有写出，数字补零，精确值0.866025404，即计算值皆为有效数字。②正12边形的边长之平方“觚幂 2679亿4919万3945(平方)忽”，其精确值为二减根号三即0.2679 4919 2431平方尺。计算值后四位数字对开方求边长“觚面0.5176381尺”即精确到十分之一忽没有影响。因不需要求12边形、24边形的面积，其边长在《九章算术注》没有给出

（中略）

割四十八觚以为九十六觚术曰：亦令半径为弦，半面为句，为之求股。置次上弦幂，四而一，得四十二亿七千七百五十六万九千七百三忽，余分弃之，即句 幂也。以减弦幂，其余，开方除之，得股九寸九分七厘八毫五秒八忽十分忽之九。以减半径，余二厘一毫四秒一忽十分忽之一，谓之小句。觚之半面又谓之小股。为之求小弦。其幂四十二亿八千二百一十五万四千一十二忽，余分弃之。开方除之，得小弦六分五厘四毫三秒八忽，余分弃之，即九十六觚之一面。

注：利用公式(1L) 从48边形边长平方求96边形边长平方，最后开方得到边长0.065438尺，精确值为0.065438 17，即计算值的五位数字皆为有效。

以半径一尺乘之，又以四十八乘之，得幂三万一千四百一十亿二千四百万忽，以百亿除之， 得幂三百一十四寸六百二十五分寸之六十四，即一百九十二觚之幂也。以九十六觚之幂减之，余六百二十五分寸之一百五，谓之差幂。倍之，为分寸之二百一十， 即九十六觚之外弧田九十六所，谓以弦乘矢之凡幂也。加此幂于九十六觚之幂， 得三百一十四寸六百二十五分寸之一百六十九，则出圆之表矣。故还就一百九十二觚之全幂三百一十四寸以为圆幂之定率而弃其余分。

注：基于正多边形边长求面积，并估算误差，

S_2n= nRb_n/2 = nR[sqrt(b_n²)]/2                                        (2L)

S_2n< πR² < S_n+2(S_2n – S_n) = S_2n+ (S_2n – S_n)               (3L)

半径一尺即十寸，割圆至96边形，确定其面积在314 又625分之64 至314 又625分之169 平方寸之间，即圆周率在3.141024~3.142064, 取3.14为偏小值。面积为平方寸，但仍写为寸。 

评论：刘徽的以上工作是严谨的。不过，没有将公式(1L) 化简而逐项计算不仅增加了工作量，还可能使舍入误差积累；因没有零和小数点的概念，几何长度的数据皆有长度单位，限制了计算精度的提高。  

2   祖冲之的圆周率

《南齐书·卷五十二》 和《南史·卷七十二》 皆有祖冲之的小传，未涉及圆周率；《隋书·志十一·律历上》有：

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

2.1   至关重要的是“以圆径一亿为一丈”，割圆想来会取半径R为一亿。这固然是将刘徽的半径一尺增加到一丈，且不必考虑各位数字的长度单位，而1为刘徽的十分之一忽。祖冲之想来利用式(1Z)计算，开方运算可以达到亿以下而不必如刘徽止于十分之一忽，即有效数字的位数可以增加。

上式与式(1L)形式等价而含义不同，且减少了加法、除法、平方的计算各一次。仅一次开方计算产生舍入误差。

2.2  若割圆至96边形以上，边长精确到1即十分之一忽，利用公式(2L)求面积就会将边长的舍入误差（最大0.05忽）放大48倍或以上，那么无论采用何种修正方法都不能得到精度一忽的祖率。请注意，祖冲之并不使用小数。容我猜测，祖冲之会设法提高精度：不求边长而直接开方得到面积

S_2n= nRb_n/2  = sqrt(n²b_n²/4)    亿亿                                                (2Z)

半径R为一亿，算筹的首位为亿亿。又，式(1Z)和(2Z)两次开方运算的结果单调增加且有极限，计算量因之而减小。

2.3   既然以多边形面积逼近圆面积，那么考察多边形面积的增加过程进而估算剩余弓形的面积，是自然而然的。增加的比率D 逐渐减小而趋于1/4，是一眼可见：割圆新增三角形的底和腰都略大于原三角形的一半，而底角为一半，故而面积稍大于八分之一。因而圆面积满足

S_2n+ (S_2n – S_n)/3< πR² < S_n+ (S_2n – S_n) * D_n /(1 – D_n)              (3Z) 

祖冲之利用冬至、夏至前后大致相同日子的日影准确判断交节时刻，并为元代郭守敬采用，见《元史·历一》，表明祖冲之知道插值和外推，提出公式(3Z)是完全可能的。又，下限就是以抛物线替换圆弧。

2.4   “以圆径一亿为一丈”以及上列三个公就是祖冲之对刘徽割圆术的发展。算筹首位为亿亿，其后有十一位数字，则割圆五次至192边形即可得到祖率。现在有计算器和计算机，演算是容易的。在以十二位数字计算时，增加割圆次数，边长平方值变小、有效位数减少，割圆六次不能提高圆周率精度，且下限已经超过真值。

如果注意 D 与0.25 的差异是以1/4 的比率趋于0，那么以十位数字割圆四次至96边形，以D₁₉₂ = 0.25005 估算192边形的S_2n，也可以得到祖率。

至于密率355/113可逐次逼近：1/8=0.125，1/6=0.166667，而1/7=0.142857…， 3又1/7 即22/7是圆周率的偏大估计。

将1/7 修改为 P = 1/(7+1/n) 即 n/(7n+1) 在1/8~1/7 之间，而

n=10，P =0.1408451;  n=20，P =0.1418440

n=15，P =0.1415094;  n= 17，P =0.1416667

最后可确定n=16，P =16/113=0.1415929 

3    三十年前写作拙稿[1]时只是参考文[2]，未能阅读《九章算术》，部分叙述稍欠准确，但结论似仍能成立。拙稿发表于二十年前，也曾在博客介绍，或许尚没有得到学界的注意，不揣浅陋，再作申述：祖冲之不会使用小数点和零，《隋书》所记“以圆径一亿为一丈”，意味着直接利用刘徽的割圆术不可能得到祖率——割圆11次至12288边形并非事实；祖冲之能够利用插值和外推确定冬至和夏至的交节时刻，因而能够完成前述计算。

[1]  尤明庆．割圆术确定圆周率方法的改进——祖冲之确定圆周率过程之猜测. 安阳师范学院学报，2003,(2)：11~12, 14

[2]  夏道行. π和e. 上海教育出版社，1964.

[3]  王振东，姜楠. 祖冲之与圆周率. 力学与实践，2015, 37(3):404-408

事非经过不知难——对祖冲之确定圆周率之揣测.2014.3.14

https://blog.sciencenet.cn/blog-275648-776036.html

放大和缩小的计算技巧：以确定圆周率为例. 2017.9.28

https://blog.sciencenet.cn/blog-275648-1078100.html

割圆术或为阿基美德（BC287-BC212）首创，刘徽和祖冲之计算圆周率固然是杰出的工作，但真不必强调祖冲之领先于西方千年以上——历史长河中不知堙灭凡几。

德国人鲁道夫·科伊伦(1540-1610)仍以割圆法确定圆周率到小数点后35位，或许真是穷毕生之力。其后牛顿(1643-1727)和莱布尼茨(1646-1716)发明微积分，引入级数求解圆周率的公式。如反正切函数可展开为

atan x = x– x³ /3+ x⁵/5– x⁷/7 + x⁹/9 – ……=A(x)

π= 4 A(1) 收敛较慢是众所周知；不过，tan α=1/2，tanβ=1/3，则 tan (α+β) =1，故而 π=4A(1/2) +4A(1/3);又因 A(1/2) = A(1/3)+A(1/7), 有 π = 8A(1/3) +4A(1/7) 可以提高计算效率。类似的公式很多。

Machin 在1706年以 公式 π = 16 A(1/5) – 4A(1/239) 将圆周率计算至百位，或许最为著名。笔者觉得利用 A(1/5) = 2A(1/10) –A(1/515) 将其转换为

π = 32A(1/10) – 4 A(1/239) – 16 A(1/515)

手工计算时具有更高的效率，或许优于上表所列公式。式中1/10 大于Gauss 公式的1/18，需要计算更多的位数，但实际工作量并不增加。这样的改进想来早有人写出，只是Machin 没有注意到有些奇怪。

现在计算机推求圆周率则会用更复杂的公式，如努马金(1887-1920)公式，其首项有七位准确数字，其后每计算一项增加八位准确数字，网上多有介绍。

4     四十四年前所购书籍和30年前所作计算。“夏道行.π和e” 1964年第1版，1978年重印，3万6千字，价0.15元，恰是二两猪肉钱——其时含骨统肉一斤0.73元；以十位数字的计算器Casio fx-4500p 计算。从最后一列的面积增大的比值d可以看出有限数字的丧失。