阿基米德计算π方法的改进

第一部份：我为什么要编算本文

 

2013年8月，我写过一篇题为《重走古希腊阿基米德、古中国刘徽的路，用割圆法计算圆周率π》的文章，以自己独立推导的割圆法公式，分步、循环递归，从正六边形出发，计算得正十二、二十四、四十八边形 … 的一系列π，可以说是重走了古希腊阿基米德、古中国刘徽的路，验证了π在223/71和22/7之间。三年来，此文在新浪博客上有4600多人阅读过，其中有一条评论还说：“这真是极好的”。

但阿基米德原著我没有见到过，所以他是怎样具体计算的，我还不得而知。最近，我化50元在卓越网上买到了一本《阿基米德全集》，得以见到《圆的度量》的原文。

研读之后，一，我懂得了他计算π的思路。二，钦佩他将起算数据1、√3、2转化为三个整数153、265、306，或780、1351、1561，这就使计算方便得多。三，明白了不相似直角三角形之间，边长的比例关系，要靠一个命题来支持，这就是《几何原本》第六卷中的一个命题6-3。如果没有这个命题，就无法推算有关边长的比例关系，也就无法计算π。因为它是一个理论基础。所以我有意补进了这一命题的内容，使全部计算有根有据。四，有了阿基米德的原文，就不必重复他的计算细节了。但阿基米德的算法，还可以改进一下，使它适合编程电算。这就是我编算本文的原因。

 

第二部份：阿基米德原文《圆的度量》的特点

 

一  计算π的思路。阿基米德用外切法与内接法求π。

外切法求π 见图二

1         以一个30⁰、60⁰、90⁰的外切直角三角形⊿AOC起算。

2         已知：直角长边OA=B=265，作为圆半径。外切的直角短边AC=S=153，对应着未知的AC弧。

圆心角∠AOC=30⁰ ，是圆周的1/12。则圆周率π= 圆周/直径＝AC弧*12/2B。但AC弧不知，便以AC代替，算得近似的π=AC*12/2B。代入S、B后，得=153*12/530=3.46415。

3         为了进一步获得较精确的π，将圆心角中分，得到外切直角三角形⊿AOD，∠AOD=15⁰，同

理，圆周率π=圆周/直径＝Ad弧*24/2B。但Ad弧不知，便以AD代替，如果知道了AD，则可算得近似的π=AD*24/2B。但AD也不知，所以关键是要在△AOC、△AOD中算出AD。

3  以后又中分出直角三角形⊿AOE，∠AOE=7.5⁰，而π= Ae弧*48/2B≒AE*48/2B。而要算π，关键是要在△AOD、△AOE中算出AE。

4  这样一次次中分下去，分到圆心角极小，则算出的π越精确。

5  《圆的度量》一文，就是在推导AO与AC、 AD、、AE、AF等等的比式。实际上是计算π的最大渐近值。阿基米德算到圆心角为3.75⁰时，即圆周的1/96时，得π=3+（667.5/4673.5）=3.1428266。大于π。

                        内接法求π，见图五、图六

1         阿基米德又用一个30⁰、60⁰、90⁰的内接直角三角形作起算。

2         已知：斜边AB为直径2R=B=1560，直角短边为弦BC=S=780，对应着未知的BC弧。圆周角

∠BAC=30⁰ ，圆心角∠BOC=60⁰ 为圆周的1/6，则圆周率π= 圆周/直径＝BC弧*6/B。但BC弧不知，便以BC弦代替，算得近似的π=BC*6/B。代入S、B后，得近似的π=780*6/1560=3

3  为了进一步获得较精确的π，将圆周角一次次中分，逐级算得各弦BD、BE、BF  …按上述公式，算得各次近似的π 。《圆的度量》一文，就是在推导AB与BD、BE、BF、BG 等等的比式。实际上是在计算π的最小渐近值。阿基米德算到圆周角为7⁰时，即圆周的1/96时，得π=6336/2017.25，大于π。

阿基米德由此作出结论：圆周与直径的比，小于22/7大于223/71。

细细体味，阿基米德求π，与多边形逼近圆的方法，原理上相同，但还是小有区别。他仅仅用一个三角形中分、中分、再中分，就达到了目的。并且关键仅是计算直角短边。所以推算公式既简单，又容易理解。

 

二 将起算数据由1、√3、2转化为三个整数的比，√3的误差仅为0.000025或0.00000047。

一个30⁰、60⁰、90⁰的直角三角形，作为推算π的起算，其对应的边长之比为1、√3、2 。这个道理在古希腊时代早已知道。但3的开方比较麻烦。阿基米德没有直接写出√3的具体数值，而是将√3表示为两个整数的比，即265:153 =√3。反算可得√3=1.73202614379，它与现代

√3=1.732050807568877 … 相比，仅差0.000025，真是了不起啊。这样，一个30⁰、60⁰、90⁰的直角三角形，对应的边长比例便为153、265、306，就作为推算π的实用起算数据。

更有高招。阿基米德在内接三角形求π时，却又将√3表示为1351:780即√3=1.73205128205，它与现代√3相比，误差更小，仅为0.00000047。

写到这里，我顺便要自我点赞一下。原来，还在我只读到265:153，而没有读到1351:780之前，我便停下来，思考了一下。想，有没有更精确的两个整数之比，更接近√3呢？就在电子表格上凑算，很快发现，还有6对整数之比，更接近√3，其中1351:780与2340:1351的误差最小。见下表：

 

S=直角短边	B=直角长边	√3 ≒ B / S	误差
阿采用  153	阿采用  265	1.73202614379	0.00002500
209	362	1.73205741626	0.00000670
306	530	1.73205741626	0.00000661
362	627	1.73204419889	-0.00000661
571	989	1.73204903677	-0.00000177
780	1351	1.73205128205	0.00000047
1351	2340	1.73205033308	-0.00000047
		精确√3=1.73205080756887

 

我当时就用780、1351、1560为起算数，又算了一遍π，结果当然小有变化，但不影响π应小于22/7的结论。后来我翻过几页，这才发现，阿基米德在内接法求π时，早已采用780、1351、1560为起算数了。我禁不住既高兴我的发现，又叹息我的迟到。因2000年的前人早已捷足登先了。但是我还是禁不住要点赞一下自已，毕竟79岁了，而思维、计算能力还在呢。

 

三  计算不相似直角三角形边与边比例关系的理论基础：《几何原本》的一个命题6-3。

在整个计算过程中，要知道三角形中边与边的比例关系。一般而言，只有在三角形相似的条件下，边与边才有比例关系。但是，当一个角度被中分后，分出来的两个三角形不相似，那么边与边有没有比例关系呢？如果没有比例关系，那将无法再计算下去。幸好，阿基米德的前辈欧几里德，在他的《几何原本》第六卷中，有一个关键命题6-3。说：当三角形的一个角度被中分后，存在一个比例关系X/Y=C/B。阿基米德一开始就引用了这个命题，因而推算出了三角形中边与边的各种比例关系，才得以算出π。查《几何原本》，得命题 6-3。见图一。

http://s2/bmiddle/0033AirQgy71zHcPpGVd1&690

命题 6-3  △FOA中，OD中分∠FOA，即∠1=∠2，则有X/Y=C/B

证明：

1  由A引平行于D O的直线，与F O的延长线交于K。

2  由于D O∕∕A K，所以∠3=∠2、∠1=∠4、则∠3=∠4、OA=OK=B

3  △FOA与△FKA相似，所以C/X=(C+B)/(X+Y)   → C(X+Y)=X(C+B)   → CX+CY=CX+BX   →

→ CY=BX   → X/Y=C/B   证毕。

4  由X/Y=C/B也得B/Y=C/X   注意，其中C/X=(C+B)/(X+Y)特别有用处。

 

四  阿基米德算法可以改进为循环计算

我在阅读过程中觉得，阿基米德在找边与边的比例关系时，虽然严密正确，但似乎也繁了一些，转弯多了一些。每中分一个三角形，都要重新去找这些比例关系。其实，在整个计算过程中，关键是怎样由上一个直角三角形的三条边S、B、C，计算出下一个直角三角形的三条边Y、B、E。并将Y、B、E作为又下一级三角形的S、B、C，这就构成了一个递归。如果整理出它们的转换递归公式，就不必每次去找边与边的比例关系，则整个计算就有规律可循，乃至可以编电算程序。于是我就整理出“由上一个直角三角形的三条边，计算出下一个直角三角形的三条边”的公式，并一一计算，得到正确验证。现在，我就按照阿基米德的思路、图形，但按照我的理解，将阿基米德的算法，改动一下，也可说改进一下，并把我的计算过程记录一下、发表一下，算是我学习阿基米德算法的一个心得吧。

 

第三部份：阿基米德计算π方法的改进

 

一 外切法求π。见图二、图三。

 

1  设OA为任意圆的半径。AC是过A的切线，⊿AOC中，∠AOC=30⁰

2  则 直角三角形短边AC=S=153  长边OA=B=265  斜边OC=C=306 。 即S：B：C=1：3√：2

3  二等分∠AOD，中分线交AC于D。∠AOD=15⁰ 。设 AD=Y  DC=X   AC =X+Y =S   OD=K

4  根据《几何原本》命题 6-3，则有X/Y=C/B 、B/Y=C/X 、C/X= B/Y =(C+B)/(X+Y) →

    B/Y =(C+B)/S  

5  于是Y = (S*B) / (C+B)=153*265/(306+265)=40545/571=71.007  (注意：Y就是下一个三角形⊿AOD的短边S )

6  K =√(Y * Y + B*B)=√ (71.007²+265² )=274.348  (注意：K就是下一个三角形⊿AOD的斜边C )

7  ∠AOD=15⁰ ，它所对的圆弧为圆周360⁰的1/24，设T=24。

8  现用AD=Y作为近似圆弧，则 π=Y*T/(2B)=71.007*24/530=3.215412

http://s7/bmiddle/0033AirQgy71zHsSrfE26&690    http://s9/bmiddle/0033AirQgy71zHDzaI008&690

第一个三角形算完。逼近一步，将⊿AOD再中分，见图四。

1         ⊿AOD中，∠AOD=15⁰

2          短边AD = S=71.007  (就是上一级的Y)

    长边OA=B=265 不动

    斜边OD =C=274.348  (就是上一级的K)

3  二等分∠AOD，中分线交AD于E ，∠AOE=7.5⁰  设AE=Y   ED=X  A D =X+Y =S   OE=K

4  根据《几何原本》命题 6-3，则有X/Y=C/B 、B/Y=C/X 、C/X= B/Y =(C+B)/(X+Y) →

  B/Y =(C+B)/S  

5  于是Y = (S*B) / (C+B)= 71.007*265 / (274.348+265)=18816.856/539.348=34.888  (注意：Y就是下一个三角形⊿AOE的短边S )

6  K =√(Y * Y + B*B)= √(34.888² + 265²)=267.287 (注意：K就是下一个三角形⊿AOE的斜边C )

7  ∠AOE=7.5⁰ ，它所对的圆弧为圆周360⁰的1/48，设T=48。

8  现用AE=Y作为近似圆弧，则 π=Y*T/(2B)= 34.888*48/530=3.159680

http://s12/bmiddle/0033AirQgy71zHTltDd9b&690

再步步逼近得：

次  圆心角  长边B  短边S   斜边C     求 Y=          求K=           求π=

I    ∠      B       S        C       (S*B) / (C+B)   √(Y * Y + B*B)    Y*T/(2B)  

1    15      265     153         306           71.007        274.348           3.215412

2   7.5     265     71.007     274.34      34.888        267.287           3.159680

3   3.75    265     34.888    267.287     17.369        265.569           3.146106

4   1.875   265     17.369    265.569     8.675         265.142           3.142734

…   …     …      …       …          …           …                …

8   0.1171875  265   1.084    265.002     0.542         265.0006          3.141617

16                                                                                   3.1416124796971

24                                                                                  3.14161247963026

36                                                                                   3.14161247963026

可见，到24 次后，π 就基本隐定在3.14161247963026，近似为3+1/7.06。

   我顺手编一个内接求π的Visual Basic程序，完成上述计算：

Private Sub Form_Click() ：Form1.Width = 11520：Form1.Height = 15360

 K = 306 ： B = 265： T = 24 ：G = 15 ： S = 153

  For I = 1 To 8 ：C = K： Y = S * B / (C + B)：K = Sqr(Y * Y + B*B) ：P = Y * T / 2/B

  Print "I="; I; "  G="; G; "  C="; C; "  S="; S; "  Y="; Y; "  K="; K; "  P="; P

   S = Y： T = T * 2: ：G= G / 2

 Next I  Print Spc(4); "成功"  

End Sub

二  内接法求π。见图五、图六。

 

1  设AB=2R是圆的直径。且设AC交圆于C，圆周角∠CAB=30⁰。连接BC。

2  则 直角三角形短边BC=S=780  长边AC=C=1351  斜边AB=B=1560 ，即S：C：B=1：3√：2

3  二等分∠CAB，中分线交BC于d ，交圆于D。连接BD。∠CAd=∠DAB=15⁰ 。即分出了两个相似三角形△CAd与△DAB。

 设Cd=X   dB=Y   C B=X+Y =S  DB=F   Ad=E

4  根据《几何原本》命题 6-3，则有X/Y=C/B 、B/Y=C/X 、C/X= B/Y =(C+B)/(X+Y) →

    C/X= B/Y =(C+B)/(S)   → C/X =(C+B)/S  

5  于是X = (S*C) / (C+B)=780*1351/(1351+1560)=105378/2911=362 

6  E =√(X * X + C*C)=√(131044²+1351²) =1398.658

7 △CAd与△DAB相似，所以X/F = E/B  有

F = X * B / E=362*1560/1398.658=403.758  (注意：F就是下一个三角形△DAB的短边S )

8  ∠DAB=15⁰ ，与它同弧的圆心角∠DOB=30⁰ ，它是圆周360⁰的1/12，设T=12。

现用BD=F作为近似圆弧，则 π=  F * T / B=403.758*12/1560=3.105828

9  为了后续的计算，预算AD=√(B*B－F*F)=√(1560²－403.758²)=1506.844

(注意：AD就是下一个三角形△DAB的长边C )

10  至此，由第一个内接三角形来计算π的任务已经完成，并为下一个三角形的计算安排好了起算数据，即△DAB的短边S=403.758  长边C=1506.844 ，斜边B=1560

 

http://s2/bmiddle/0033AirQgy71zI5UGWt21&690      http://s5/bmiddle/0033AirQgy71zIn1Yfac4&690

 

第一个三角形算完。逼近一步，将△DAB再中分，见图七。

1  己知△DAB的短边S=403.758  长边C=1506.844 ，斜边B=1560

3  二等分∠DAB后，中分线交BD于e ，交圆于E。连接BE。∠DAe=∠EAB=7.5⁰ 。又分出了两个相似三角形△DAe与△EAB。

设De=X   eB=Y  D B=X+Y =S  EB=F   Ae=E

4  根据《几何原本》命题 6-3，则有X/Y=C/B 、B/Y=C/X 、C/X= B/Y =(C+B)/(X+Y) →

    C/X= B/Y =(C+B)/(S)   → C/X =(C+B)/S  

5  于是X = (S*C) / (C+B)= 403.758*1506.844/(1506.844+1560)=198.380 

6  E =√(X * X + C*C)=√(198.380²+1506.844²) =1519.847

7 △CAd与△DAB相似，所以X/F = E/B  有

F = X * B / E=198.380*1560/1519.847=203.621  (注意：F就是下一个三角形的短边S )

8  ∠EAB=7.5⁰ ，与它同弧的圆心角∠EOB=15⁰ ，它是圆周360⁰的1/12，设T=24。

现用BE=F作为近似圆弧，则 π= F * T / B=203.621*24/1560=3.132628

9  为了后续的计算，预算AD=√(B*B－F*F)=√(1560*1560－203.6208*203.6208)=1546.65

(注意：AD就是下一个三角形的长边C )

10        至此，由第二个内接三角形来计算π的任务已经完成，并为下一个三角形的计算安排好了起算数据。就不再类似计算下去了。又编了一个内接求π的Visual Basic程序：

 

Private Sub Form_Click() ：Form1.Width = 11520 ：Form1.Height = 15360

  B = 1560 ： C = 1351 ：S = 780 ：T = 12 ：G = 30

  For I = 1 To 8 ：X = S * C / (C + B) ： E = Sqr(X * X + C * C) ： F = X * B / E ： P = F * T / B

    Print "I="; I; " G="; G; "  C="; C; "  X="; X; "  E="; E; "  F="; F; "  P="; P

    T = T * 2 ：G = G / 2 ：C = Sqr(B * B - F * F) ：S = F

  Next I ：Print Spc(4); "成功" ：Print

End Sub

http://s2/bmiddle/0033AirQgy71zIun9fPe1&690

运算结果如下

 

次 圆心角   斜边    短边         长边      求 X=              求E=          求F=        求π=

I     G       B          S              C     (S*C) / (C+B)   √(X² + C²)    X * B / E    F * T / B

1    30     1560   780           1351       362           1398.658       403.758    3.105828

2    15     1560   403.758   1506.84    198.380      1519.847      203.621    3.132628

3    7.5    1560   203.621    1546.65   101.373      1549.973      102.029    3.139350

4    3.75   1560   102.029    1556.65   50.960       1557.494      51.042     3.141032

…   …    …      …         …      …           …          …         …

8 0.1171875  1560   3.1907    1559.9869   3.191      1559.9902     3.191     3.14159008

…   …    …      …

16                                                                                                 3.14159227214488

24                                                                                                3.14159227217831

到24 次后，π就基本隐定。现取π=3.14159227217831，近似为222/70

阿基米德逼近4次，相当于96正多边形求π=6336/2017.25＞3+(10/71)。或＞223/71

 

三  结 论

 

经外切、内接逼近，证实π在3.14159 22721 7831与3.14161 24796 3026与之间。而历代有：

阿基米德推算的   223/71＜π＜22/7

祖冲之推算的   疏率(或作约率)为22/7，密率为355/113

刘徽推算的   157/50  ＜π＜  3927/1250

 

现在，圆周率有2061亿位精度。实用上取π=3.14159  26  53589 已足够了。有谐音帮你记忆：

山巅一寺一壶酒（3.14159），儿乐（26），我三壶不够（53589），… … …