鄒暢根　鄒源清

(南京市北京東路小學　210008)

兒子幼兒園畢業(yè)后，帶他到南京東郊風景區(qū)去玩，在景區(qū)的歷史文化展覽館，我們看到了祖沖之的雕像及關(guān)于他的介紹.我開玩笑讓小孩拜祖沖之為師，記住圓周率π在3.1415926到3.1415927之間，并要求他在7歲之前搞清π的含義.

為了搞清楚π，我?guī)『⒁黄鹩眉毦€量了樹干的周長，并在陽臺上畫了個半徑為1米的半圓如上圖，用細線測量π的值(誤差巨大).

原本要求小孩能理解圓周率π為圓的周長除以直徑，但后來發(fā)現(xiàn)這個要求太低，我改變計劃，要求他能用計算器算出π.

在小孩上小學一年級的時間里，教了他簡單的四則運算，以及三角形特征、勾股定理等簡單的幾何學知識，重點是要小孩初步理解這些基礎(chǔ)知識，以便能理解祖沖之用割圓術(shù)計算圓周率的方法.

祖沖之采用更早的數(shù)學家劉徽的割圓術(shù)如右圖，以圓內(nèi)接正多邊形的面積及其外加一個矩形面積之和，雙向逼近圓的面積.

假設(shè)半徑為1的圓，它的內(nèi)接正n邊形的邊長AB為Xn、面積為Sn，內(nèi)接正2n邊形的邊長BC為X2n、面積為S2n，假設(shè)三角形ABC的面積為S△ABC，矩形ABEF的面積為S四邊形ABEF，圓的面積為S，那么

n·S△ABC=S2n-Sn，

n·S四邊形ABEF=2·(S2n-Sn)

由于圓的面積大于圓內(nèi)接正2n邊形面積、小于圓內(nèi)接正n邊形面積加n個矩形ABEF的面積，所以S2n

利用三角形面積公式，可以計算出四邊形AOBC的面積為AB·OC/2，所以

S2n=Xn·n/2 ，Sn=Xn/2·n/4 ，

因此圓周率π的范圍如下：

Xn·n/2<π

由于籌算比較費時費力，生于數(shù)學世家的祖沖之很清楚節(jié)省運算量的重要性.上述推理僅用到了勾股定理，祖沖之利用數(shù)學家劉徽的割圓術(shù)計算圓周率，應(yīng)該很容易推出上述運算量最少的計算方法，但越不過多位數(shù)的開方.

我們準備用計算器開方來算出π，但是失敗了，用計算器居然無法計算出祖沖之算出的π，經(jīng)仔細分析后發(fā)現(xiàn)，使用的計算器是12位的(例如3.12345678901，總計12位)、精度不夠，因此決定在網(wǎng)上買個16位的計算器.

網(wǎng)購的計算器需要3天后才能到貨，期間正好夾了個星期六和星期天，我決定帶兒子一起分析祖沖之是如何做運算的.

經(jīng)互聯(lián)網(wǎng)搜索發(fā)現(xiàn)：祖沖之生活的年代已經(jīng)用十進制計數(shù)，但阿拉伯數(shù)字還沒被傳入中國，算盤也沒有被發(fā)明出來，計算工具是小棍子，這些小棍子稱為算籌，用算籌做運算的過程稱為籌算.祖沖之寫了一本書叫《綴術(shù)》，其中記錄了對圓周率的研究及成果，但后來這本書失傳了，祖沖之如何精確計算出圓周率成為未解之謎.小孩和我用計算器都算不出祖沖之用小棍子算出的π，我們更加佩服祖沖之！中國科協(xié)主辦的《科技導(dǎo)報》2008年第5期上，將“祖沖之究竟是怎樣計算出圓周率π值的？”列為公眾關(guān)注的18個未解科學難題之一.互聯(lián)網(wǎng)上有不少分析祖沖之是如何計算π的內(nèi)容，絕大部分都認為計算很困難.祖沖之計算π真像有些文章描述的那樣“一雙手被磨出了厚厚的老繭”嗎？我決定帶兒子用小棍子來試一試.

我們在擺弄小棍子的過程中，16位計算器到貨了，借助于它幾分鐘就能將π算到祖沖之算出的3.1415926.經(jīng)過幾天的探索，我們發(fā)現(xiàn)用割圓術(shù)計算π的工作量遠遠沒有人們想象的那么大，如果不考慮探索研究的時間，在計算不出錯、動作快等前提下，籌算一兩天就能將π計算到祖沖之所達到的精度.

一個小于1的小數(shù)△的平方，例如0.052比0.05小很多，將它除以一個大于1的數(shù)之后則變得更小，假設(shè)對A開方，有個大于1的數(shù)B，B2與A接近、差值為△，則可以用如下逼近方法計算：

≈B+△/2B=B+λ.

利用上面的逼近方法，每次將B的精度多提升1個小數(shù)位、步長為λ，A的逼近值每次都可以精確計算出來，其值為B2+2·λ·B+λ2，每次逼近僅需對原值做2·λ·B+λ2修正，運算量很小、而且很適合籌算.古代算籌可以有些變異，例如有人將小棍子染成不同的顏色、或采用不同的放置方法避免出錯，但實質(zhì)都是一樣的.在本文中，算籌采用如下方法表示0 ～ 9這10個數(shù)字.

為了便于用算籌做開方運算，需要對開方過程中用的數(shù)值描述如下：

1.被開方值A(chǔ)；

2.A的逼近值，“A的逼近值_前”為前一次的逼近值，“A的逼近值_后”為本次逼近運算后的值；

3.開方值B，“開方值B_前”為前一次逼近的開方值，“開方值B_后”為本次逼近運算后的值，它的平方精確等于“A的逼近值_后”；

4.除數(shù)D，對2·B(用開方值B_前)做四舍五入取整(或取2·B的整數(shù)位加1)得到除數(shù)D，取值范圍為4、5、6、7、8、9、10這7個數(shù)值之一(請看下面第7小點說明)；

5.差值△，為被開方值A(chǔ)減去“A的逼近值_前”，可能為正或負，以“開方值B_后”小數(shù)點最后一位為參考位，將被開方值及其逼近值做四舍五入取整后，再計算它們的差值；

6.步長λ，為△除以D后再做四舍五入取整，一般情況下λ小于等于5.由于D、△、λ都做了取整，因此有可能出現(xiàn)少量的情況下λ大于5，此時需要在不增加開方值B小數(shù)位的基礎(chǔ)上再做一次逼近，降低被開方值與其逼近值的差值，以避免λ超過5而引入2位數(shù)的乘法運算(當然，λ大于5也沒有關(guān)系，只是籌算2·λ·B略有不便)；

7.上述開方運算，對于區(qū)間 [ 3，25 )內(nèi)的任何數(shù)，開方時僅會涉及到2位數(shù)除以1位數(shù)的除法以及多位數(shù)乘以1位數(shù)的乘法(乘除10因為很簡單不計入2位數(shù)乘除法).對于其他情況：區(qū)間(1，3)內(nèi)的數(shù)值可以先乘以4之后再開方；區(qū)間(25，100)內(nèi)的數(shù)值可以先除以4之后再開方；區(qū)間(0，1)內(nèi)的數(shù)值可以先乘以10的偶次方之后再開方；100以上的數(shù)值可以除以10的偶次方之后再開方.

只需要用2位數(shù)除以1位數(shù)的除法計算出步長λ及多位數(shù)乘以1位數(shù)的乘法、多位數(shù)的加減法，即可以按需要將開方的精度提升到任何所需的位數(shù)，這樣的運算對籌算而言比較簡單.下面舉例利用割圓術(shù)計算π過程中需要開方的一個數(shù)，例如3.732050807568877，利用籌算對它進行開方運算，這是一個16位精度的數(shù)，對它開方后也要精確到小數(shù)點后15位.

·開方值整數(shù)位計算

B取2最接近；A的逼近值為4.

·開方值小數(shù)點后第1位計算

如下圖，逼近值大于被開方值，步長λ為負，除數(shù)D=2·B=4，

差值△≈40-37=3，

λ=△/D≈3/4，取λ=1，

2·λ·B=2·1·2=4，

λ2=1，

A的逼近值變?yōu)?4-4·0.1+1·0.01=3.61，

B變?yōu)?2-1·0.1 = 1.9.

·開方值小數(shù)點后第2位計算

如下圖，逼近值小于被開方值，步長λ為正，除數(shù)D=2·B≈4，

差值△≈73-61=12，

λ=△/D≈12/4， 取λ=3，

2·λ·B=2·3·1.9=11.4，

λ2=9，

A的逼近值變?yōu)?3.61+11.4×0.01+9×0.0001=3.7249，

B變?yōu)?.9+3×0.01=1.93.

·開方值小數(shù)點后第3位計算

如下圖，逼近值小于被開方值，步長λ為正，除數(shù)D=2·B≈4，

差值△≈32-25=7，

λ= △/D≈7/4,取λ= 2,

2·λ·B= 2×2×1.93=7.72，

λ2=4，

A的逼近值變?yōu)?.7249+7.72×0.001+4×0.000001 = 3.732624，

B變?yōu)?.93+2×0.001 = 1.932.

·開方值小數(shù)點后第4位計算

如下圖，逼近值大于被開方值，步長λ為負，除數(shù)D=2×B≈ 4，

差值△≈26-21=5，

λ= △/D≈5/4，取λ=1，

2·λ·B=2×1×1.932=3.864，

λ2=1，

A的逼近值變?yōu)?.732624-3.864×0.0001+1×0.00000001=3.73223761,

B變?yōu)?.932-1×0.0001=1.9319.

從上述籌算過程可以看到，開方過程僅會涉及到如下計算：2位數(shù)除以1位數(shù)的除法及多位數(shù)乘以1位數(shù)的乘法(乘除10因為很簡單不計入2位數(shù)乘除法)，多位數(shù)的加減法.重復(fù)上面的運算很容易開方到第15位小數(shù)的精度，特別是籌算加減法無需要像我們常用的豎式運算那樣還要重復(fù)抄寫一些數(shù)字、對齊小數(shù)位也比較方便，因此熟練籌算之后，籌算速度有可能達到甚至超過我們用阿拉伯數(shù)字豎式運算.我試過籌算一次16位精度的開方大約需要1小時，祖沖之籌算的速度一定比我快很多，因此如果不考慮任何探索性時間開支，僅考慮最小運算量的情況下，祖沖之最快一天就可以將圓周率重新算一次，確定它在3.1415926與3.1415927之間.

下圖用阿拉伯數(shù)字更符合現(xiàn)代人們的理解習慣，從圖中可以很明顯地看到開方的正負雙向逼近過程.

3.000000000000000 (計算

3.732050807568877

3.931851652578136

3.982889722747621

3.995717846477207

3.998929174952731

3.999732275819123

3.999933067834802

3.999983266888701

3.999995816717800

經(jīng)過12次開方，得到：

=0.000000261455223,

=0.000000065363807.

再經(jīng)過2次開方，得到：

=0.000511326923797,

=0.000255663464343.

利用前面推導(dǎo)出計算圓周率π的不等式：

Xn·n/2<π

取n=24576，將X12288及X24576的值代入，可求得圓周率π的范圍：

3.141592649846784<π<3.141592679884800.

通過上述14次開方，計算出圓周率π介于3.1415926和3.1415927之間.上述計算結(jié)果都是我用16位精度的計算器算出的，計算過程保留小數(shù)點后15位，如果減少1位精度，那么計算出的圓周率上限值將超過3.1415927，因此祖沖之開方的時候最少要保留16位精度.

這個開方方法有點像二分法，只是引入步長λ評估從而一次增加1個小數(shù)位，同時充分利用了前面運算的結(jié)果，每一步都確保開方值的平方精確等于被開方值的逼近值.與傳統(tǒng)手工開方法相比，它不需要多次試乘而引入多次乘法(通過試乘規(guī)避了解二次方程，但反復(fù)試驗多位數(shù)乘1位數(shù)的運算，會增加許多運算量).這個開方過程雖然有較多次數(shù)的少位數(shù)加減法，但由于僅有2位數(shù)除以1位數(shù)的除法及多位數(shù)乘以1位數(shù)的乘法(乘除10因太簡單不計入2位數(shù)乘除)，因此很適合原始的籌算(當然也適合現(xiàn)代人用阿拉伯數(shù)字進行手工開方).現(xiàn)代計算機的數(shù)字運算，CPU實現(xiàn)多位乘除法(例如100位乘除100位)仍然具有困難，但實現(xiàn)多位的加減法卻可以較容易地變通實現(xiàn)，利用上述開方方法，可用CPU快速實現(xiàn)多位數(shù)的高精度開方，在某些領(lǐng)域可能會有幫助，有興趣的讀者可以自行編寫這個高精度多位開方函數(shù).

為了說明籌算圓周率π不太困難，我拍了一個簡單的視頻影象放在互聯(lián)網(wǎng)上(www.iqiyi.com和www.youku.com)，視頻的名稱為“兒子和我一起用小棍子精確計算圓周率”，內(nèi)含我和剛上完小學一年級的兒子用算籌(小棍子)進行開方運算，我假設(shè)祖沖之開方時精確到小數(shù)點后17位，讓小孩籌算一個18位有效數(shù)字(1個整數(shù)位，17個小數(shù)位)的開方值，精確到小數(shù)點后17位.