Scratch & Math: 歡樂派(三)

Peter Wei
12 min readMay 22, 2017

--

1991年2月波灣戰爭時,美國布署在沙烏地阿拉伯的愛國者飛彈防禦系統沒有成功地攔劫到伊拉克的飛毛腿飛彈,造成美軍28人死亡和超過100人受傷。事後調查報告顯示,這是愛國者飛彈系統內軟體累積了0.34秒時間誤差所造成的後果。在這一章裡我們實作圓周率的精準計算,然後來瞭解什麼是電腦軟體裡的浮點數精確度和捨入誤差。

數學關鍵字: 圓周率。

數學家:萊布尼茲、高斯、勒讓德、余智恆。

程式關鍵字: 浮點數精確度(floating point precision)、捨入誤差(round-off error)。

Scratch程式:

圓周率π的精確度

在網路上可以利用「整數數列線上大全」(On-Line Encyclopedia of Integer Sequences,縮寫:OEIS)來查詢圓周率π的數值和相關資料。「整數數列線上大全」是美國AT&T實驗室的研究員尼爾·斯洛恩(Neil Sloane)自1956年開始蒐集的各種數列,到了2016年底已經累積了大約28萬條各式各樣的數列。登錄在編號A000796的數列就是圓周率π,下面的數字顯示圓周率 π 精確度到小數點後兩百位,

圓周率π是一個無理數,也就是說如果將圓周率π寫成小數形式,小數點之後會有無限多個不會循環的小數。我們可以想像如果持續計算,上面這個圓周率π可以列出無窮無盡的精確位數,直到天荒地老也沒有個盡頭。

阿基米德用多邊形逼近法得到圓周率π的範圍介於,

在這種情況下圓周率π可能等於3.1410,也可能等於3.1420,或是任何介於上下兩個邊界值的數。由這個式子只能推斷圓周率π的精確度到小數點後第二位,也就是3.14,而小數點後第三位就無法由這個式子確定是那一個絶對的數字了。

阿基米德的時代距離今天大約2300年前,他以手算出96邊形的比例然後推導圓周率π值,這樣的計算量就已經足以讓我們頭昏眼花。現代電腦的強大運算能力成為計算圓周率π的利器,我們在動手做5–1中用scratch程式求3600邊形的比例推導出圓周率π值,彈指間就可以算出3.141592,精確度到小數點後六位,真的是又快又準。

余智恆(Alexander Yee)是在美國出生的華裔第二代,和一般的年輕人一樣喜歡打電動、聽音樂、迷日本動漫,在他22歲的時候他寫的程式 y-cruncher創下了計算圓周率π到小數點後五兆位數的世界紀錄。五兆位數是個什麼樣的概念?基本上這個數有5*1012位數,如果用A4的紙每頁500字的規格來列印要印100億張紙才印得完。

一般工程或天文運算其實並不需要用到成千上萬位精確度的圓周率 π,計算圓周率π值到小數點後39位的精確度已經足以計算出誤差小於一個原子大小的可觀測宇宙周長。像 y-cruncher 這類計算高精確度圓周率π值的程式多半應用於電腦軟硬體對於精確度和運算速度的測試。當然如果能夠打破圓周率π值精確度計算的世界紀錄當然也是挺讓人開心的一件事。

要計算這麼大量的數值資料,y-cruncher程式必須搭配多核心的運算處理器,以多執行緒(multithreading)的運算方法,在同一個時間執行多個執行緒來提升計算效能。除了強大的運算處理器外,硬體方面也會要求要有龐大的記憶體和足夠的磁碟容量來計算和儲存結果。即便如此,如此大量的資料也會須要上百天的時間進行運算。

動手做5–4 公式求解圓周率π

[程式設計需求規格]

演算法的選擇在大量運算時會是影響程式執行時間的關鍵。現代電腦程式會利用快速收斂的遞迴公式求解、快速的乘法演算法和運用多執行緒的程式架構等等,來加速計算,以求得大量而精確圓周率π值。

萊布尼茲,十七世紀德國數學家、哲學家,發展出微積分、拓樸學,並對二進位的發展做出貢獻。(取自維基百科Gottfried Wilhelm Leibniz

動手做5–4我們要利用三種不同的公式來解出圓周率π值,並且比較結果的精確度。

[公式1]萊布尼茲圓周率公式

十六、七世紀的歐洲數學家開始用無窮級數(infinite series)來計算圓周率π。無窮級數是無限個數列(series)數字的總和,用這種方法計算出的圓周率π值會比阿基米德的多邊形逼近法要來得精確。著名的德國數學家萊布尼茲(Leibniz)在1674年用下面這個無窮級數計算圓周率π值,

[公式2]高斯-勒讓德公式

基於十八世紀德國數學家高斯(Gauss)和法國數學家勒讓德(Legendre)的研究,1975年史丹福大學的教授布蘭特(Richard Brent)和另一位數學家沙勒明(Eugene Salamin)各自推導出了下面的演算法來計算圓周率π值。這個演算法分為二位步驟,

  1. 初始值設定,

2. 重複下面計算式直到an和bn的差值小於想要的範圍,

3. 計算圓周率π值,

[公式3]連分數

圓周率π是一個無理數,因此無法表示成一個整數的比例,但是可以表示成如下的連分數(continued fractions),

接著我們就要將上述三種公式用scratch實做出來,然後比較圓周率π的精確度。

[程式設計角色和舞台]

在這個程式裡只會使用顯示變數而隱藏所有角色,可以使用任意造型的角色。舞台上以文字輸入三種公式選項。

放大的舞台背景1畫面如下,

[程式設計解決方案]

[程式檢視]

這個程式主要是將不同的數學公式轉換為程式語言,我們來審視一下每一段程式的細節和結果。

[公式1]萊布尼茲圓周率公式

萊布尼茲公式是一個無窮級數,我們依據它的特徵轉化為程式,

  1. 無窮級數中每一個數字前的正負號呈現規律地一正一負變化。設定變數n為由1開始,表示第n項數字 ,偵測n除以2的餘數就可以判斷變數正負號。

2. 無窮級數是由一連串的分數組成的數列和,數列裡的分數分母是逐次遞增的奇數。

3. 程式只能計算有限多項數字和,由於分數的分母越來越大,代表越後面的分數值會越小終至可以忽略,所以在這裡我們總共計算100,0000項數字和來逼近無窮數列的值。最後要記得乘以4來得到圓周率 π 值。

執行程式可以求得圓周率π=3.1415936535…,

也就是說程式計算了一百萬個數列數字和,可以計算出圓周率π值精確度到小數點後五位,可以看到這並不是一個可以很快收斂圓周率π值的方法。事實上根據數學推算,用萊布尼茲圓周率公式要計算五十億個數列數字和,才可以精確到小數點後第十位。

[公式2]高斯-勒讓德公式

高斯-勒讓德公式裡我們建立三個清單:

a清單存放an和bn的算數平均數,

b清單存放an和bn的幾何平均數,

t清單存放相對映的tn計算值,

值得注意的是scratch裡沒有2的次方指令方塊,我們必須利用10的次方來轉換,

高斯-勒讓德公式相對是一個收斂很快的公式,運算3次就可以得到圓周率π=3.14159265358979,精確度達到小數點後第十四位。

[公式3]連分數

圓周率π的連分數可以拆成第一層整數3加上一個連分數的部分,

從第二層開始都是整數6加上一個連分數,且數字分子為漸增的奇數平方。建立客製方塊做遞迴運算。

同樣地程式只能計算有限多個連分數值來估算圓周率,在這裡我們總共計了200,0000層連分數。圓周率π=3.141592653589…,精確度到小數點後第十二位。

以上三種公式計算圓周率 π,高斯-勒讓德公式的收斂速度最快,只要計算3次就可以使精確度到小數點後第十四位。但是與像y-cruncher這樣的程式相比,這樣的精確度實在不夠看。

我們計算出的圓周率π值和實際圓周率 π 間的差值被稱為計算誤差,基本上在我們的scratch程式有兩種誤差來源,

1. 以有限數列來逼近無窮級數,省略沒有計算的數列數字會造成誤差。

2. 捨入誤差(round-off error),電腦是以有限大小的位元數來儲存變數,所以當電腦要儲存一個有無限多位小數點後數字的數時,就會有捨入誤差。

像y-cruncher可以計算精確度到小數點22兆位數的程式,通常要應用特別的「任意精確度計算法」來做加減乘除四則運算來避免任何可能的誤差。

觀念平台:浮點數精確度(floating point precision)與捨入誤差(round-off error)電腦以一種稱為浮點數(floating point)的形式儲存小數的數值。電腦的浮點數看起來有點像科學記號表示法,例如說十進位的數字-0.875可以表示成:
64位元的電腦會將-0.875分為三個部分做儲存,1. 正負號(sign):在這裡會存下數字1,使得 (-1)¹=-1。如果是正數則會存下數字0,使得 (-1)⁰=1。2. 有效數字(fraction):存下數字7。3. 冪次(exponent):存下數字1020。這三個部分的資料會用一個64位元的空間來儲存,
(取自維基百科Double-precision_floating-point_format, GFDL
儲存的空間總共有64位元,每一個位元都是0或1的二進位數字。當電腦用有限位數的的浮點數來進行實數的四則運算,與相應的數學運算相比很可能會出現精確度上的誤差。舉例來說,數學裡可以出現很小很小的正實數,這個數還是大於0。但是電腦以64個位元來儲存一個正實數卻有最小的下限,超過這個下限就會被視為0,這種狀況就稱為算數下溢(underflow)。同樣的道理,電腦以64個位元來儲存一個正實數也有最大的上限,超過這個上限就會被視為無窮大,這種狀況稱為算數上溢(overflow)。我們可以用一個簡單的小程式來測試一下scratch變數可以顯示的最小正實數值,超過這個下限,scratch變數就會變為0。
在這個程式裡我們設定變數n初始值為1,再反覆將變數n除以2。當變數n的值超過儲存數值下限時就會被視為0。
從執行結果看當變數n反覆除以2,超過1023次之後變數n就被視為0,我們得到的最小變數n的值為:
這個值和理論上64位元電腦可以儲存的最小正實數值(1*2^–1023)值很接近但略有出入,精確度到小數點後12位。這個微小的差值就是scratch在做四則運算時累積產生的捨入誤差(round-off error)。
可別小看這個微小的捨入誤差,當它們累積起來的時候就會對數值計算產生顯著的影響,在真實的世界裡這樣的影響可能造成的巨大的財務、甚至是生命的損傷。一個著名的案例是在1991年2月波灣戰爭時,美國布署在沙烏地阿拉伯的愛國者飛彈防禦系統沒有成功地攔劫到伊拉克的飛毛腿飛彈,造成美軍28人死亡和超過100人受傷。事後調查報告顯示,這是愛國者飛彈系統內軟體累積了微小的時間誤差所造成的後果。愛國者飛彈系統內的軟體以24位元儲存數字資料,系統內計算時間的單位是1/10秒,但由於捨入誤差的關係存在0.000000095秒的誤差。在系統連續運算100個小時後,對於時間計算累積的捨入誤差達到0.34秒。而伊拉克的飛毛腿飛彈每秒飛行速度為1676公尺,因此這致命的0.34秒誤差就讓愛國者飛彈產生500公尺以上的距離誤差而攔劫失敗。

延伸閱讀:

Scratch & Math: 天花板上的蜘蛛(一) Scratch & Math: 天花板上的蜘蛛(二)

Scratch & Math: 直線裡的宇宙觀(一) Scratch & Math: 直線裡的宇宙觀(二)

Scratch & Math: 不能說的祕密(一) Scratch & Math: 不能說的祕密(二)

Scratch & Math: 無理的道理(一) Scratch & Math: 無理的道理(二) Scratch & Math: 無理的道理(三)

Scratch & Math: 歡樂派(一) Scratch & Math: 歡樂派(二)

Scratch & Math: 歡樂派(三) Scratch & Math: 歡樂派(四)

--

--

Peter Wei

DIY, Scratch, Math, Essay, Book, Travel, Movie, Music