2011年7月1日 星期五

ABC FORTH的新發展

ABC FORTH的新發展


曾慶潭 Ching-Tang Tseng
ilikeforth@gmail.com
Hamilton, New Zealand
2 July 2011


1. 前言

這可能又是一篇沒人愛看的文章了,其中卻蘊含了不少難以發展的技術,您愛不愛看就悉聽尊便,我不會很在乎,您若還想繼續讀下去,千萬別受我個人意念的影響,否則得不到好結果。

我的各種研究與發展,長久以來就不被重視,沒人理睬,更不可能得到任何支持,我早已習慣,處之泰然,唯一覺得欣慰的事情,就是發佈網文之後,還真有不少人來拜訪,全世界都有,但未必是因我的研究與發展而來。

事實上,我個人在 ABC Forth 數學計算系統上的發展,已經與大眾漸行漸遠。當初為了顧及系統的實用性,曾將 ABC Forth 系統發展之事,暫停了一段時日,改為努力地寫出一本能讓大眾接受的書本。結束寫書之長期工作後,我便再度開始大肆發展這個系統,第一件挑選的重任,就是實現我曾經在『玄數』網文中介紹過的技術,將超大整數的數學計算能力加進系統,實現了!耗時約兩星期,還有後續發展,準備再收集幾個超級大整數的特殊函數功能程式,將它們全加進系統。

這些新發展,與我已經寫成的書本不相關;這些新功能,也與一般大眾的使用需求不相關。我長期思考過這些問題,所以發展再難,也應該可以一一突破;我也做好心理上的準備,所以完成一項,就還會冒出新的一項。每個新發展,都需要讀不少書,也要研究許多先聖先賢設計出來的大型程式,全面融會貫通之後,才能組織出自己的東西。能有具體結果的發展,都不是輕鬆的工作,如果您親眼目睹我這個60歲的老頭,這麼老了還這樣子幹活,您適合在年輕時及時享樂而懈怠學習嗎?希望我的努力能夠成為年輕人的鼓勵。

首先,讓大家見一見超大整數數學計算系統的基本功能,以網文容量容得下,真實資料刊得出來為原則,利用將圓周率 Pi(π,為打字便利計,本文一律簡稱 Pi)算到 2500 位數的計算技術,展示程式與執行結果。

我對大數計算系統性能上的設計基本要求是:程式必須讓高中生也能讀懂,除了清楚還必須精確。清楚,就能讓大家容易利用這個系統來設計大數的數學計算應用程式,精確,才能作為大數問題的分析工具。

本文除了展示,免不了必須談一點理論根據與比較分析,這些內容可以告訴大家,新聞報導將Pi計算到五兆(5T digits)位數,只是技術不錯,卻並不是甚麼天大的突破,所以,技術值得鼓勵,卻值不上因此而應頒給實現者任何學術學位。使用 Forth 程式語言,能夠輕而易舉的估算出計算這種問題的記憶體耗用量,我在發展出 ABC Forth 系統的大數計算能力後,便有資格告訴大家,現行個人電腦的能力,大約就是能將Pi處理到達十兆位數,以後呢?以後再說。

2. 程式及其執行結果

\ 求 Pi 至 2500 位數的程式,以 2501 位整數來表示

INTEGER I
INTEGER KK

0 BIGVARIABLE X
0 BIGVARIABLE DD
0 BIGVARIABLE PI#1 2500 ALLOT
0 BIGVARIABLE W# 2500 ALLOT
0 BIGVARIABLE ATAN# 2500 ALLOT
0 BIGVARIABLE PI#2 2500 ALLOT

0 BIGVARIABLE BIG0
10 BIGVARIABLE BIG10

: RESET BASIC
10 LET B{ W# = BIG0 }B
20 LET B{ ATAN# = BIG0 }B
30 LET B{ DD = BIG0 }B
40 LET KK = 0
50 END
;

: SUBATAN BASIC
10 REM Subroutine ATAN(X), in fact, X=1/N
20 LET B{ DD = X * X }B
30 LET B{ W# = X }B
40 FOR I = 1 TO 2500
50 LET B{ W# = W# * BIG10 }B
60 NEXT I
70 LET B{ ATAN# = BIG0 }B
80 LET KK = 1
90 LET B{ W# = W# / DD }B
100 LET B{ ATAN# = ATAN# + ( W# / I>BIG ( KK ) ) }B
110 LET KK = KK + 2
120 LET B{ W# = W# / DD }B
130 LET B{ ATAN# = ATAN# - ( W# / I>BIG ( KK ) ) }B
140 LET KK = KK + 2
150 IF B{ W# <> BIG0 }B THEN -90
160 END
;

: PAI1 BASIC
10 LET B{ X = BIG0 }B
20 LET B{ PI#1 = BIG0 }B
30 RUN RESET
40 LET B{ X = I>BIG ( 5 ) }B
50 RUN SUBATAN
60 LET B{ PI#1 = PI#1 + ( I>BIG ( 16 ) * ATAN# ) }B
70 RUN RESET
80 LET B{ X = I>BIG ( 239 ) }B
90 RUN SUBATAN
100 LET B{ PI#1 = PI#1 - ( I>BIG ( 4 ) * ATAN# ) }B
110 RUN CR CR
120 PRINT " 2500 digits pi = : "
130 RUN PI#1 BIG.
140 END
;

: PAI2 BASIC
10 LET B{ X = BIG0 }B
20 LET B{ PI#2 = BIG0 }B
30 RUN RESET
40 LET B{ X = I>BIG ( 57 ) }B
50 RUN SUBATAN
60 LET B{ PI#2 = PI#2 + ( I>BIG ( 176 ) * ATAN# ) }B
70 RUN RESET
80 LET B{ X = I>BIG ( 239 ) }B
90 RUN SUBATAN
100 LET B{ PI#2 = PI#2 + ( I>BIG ( 28 ) * ATAN# ) }B
110 RUN RESET
120 LET B{ X = I>BIG ( 682 ) }B
130 RUN SUBATAN
140 LET B{ PI#2 = PI#2 - ( I>BIG ( 48 ) * ATAN# ) }B
150 RUN RESET
160 LET B{ X = I>BIG ( 12943 ) }B
170 RUN SUBATAN
180 LET B{ PI#2 = PI#2 + ( I>BIG ( 96 ) * ATAN# ) }B
190 RUN CR CR
200 PRINT " 2500 digits pi = : "
210 RUN PI#2 BIG.
220 END
;

ok

PAI1

2500 digits pi = :
2501 digits
31415926535897932384626433832795028841971693993751
05820974944592307816406286208998628034825342117067
98214808651328230664709384460955058223172535940812
84811174502841027019385211055596446229489549303819
64428810975665933446128475648233786783165271201909
14564856692346034861045432664821339360726024914127
37245870066063155881748815209209628292540917153643
67892590360011330530548820466521384146951941511609
43305727036575959195309218611738193261179310511854
80744623799627495673518857527248912279381830119491
29833673362440656643086021394946395224737190702179
86094370277053921717629317675238467481846766940513
20005681271452635608277857713427577896091736371787
21468440901224953430146549585371050792279689258923
54201995611212902196086403441815981362977477130996
05187072113499999983729780499510597317328160963185
95024459455346908302642522308253344685035261931188
17101000313783875288658753320838142061717766914730
35982534904287554687311595628638823537875937519577
81857780532171226806613001927876611195909216420198
93809525720106548586327886593615338182796823030195
20353018529689957736225994138912497217752834791315
15574857242454150695950829533116861727855889075098
38175463746493931925506040092770167113900984882401
28583616035637076601047101819429555961989467678374
49448255379774726847104047534646208046684259069491
29331367702898915210475216205696602405803815019351
12533824300355876402474964732639141992726042699227
96782354781636009341721641219924586315030286182974
55570674983850549458858692699569092721079750930295
53211653449872027559602364806654991198818347977535
66369807426542527862551818417574672890977772793800
08164706001614524919217321721477235014144197356854
81613611573525521334757418494684385233239073941433
34547762416862518983569485562099219222184272550254
25688767179049460165346680498862723279178608578438
38279679766814541009538837863609506800642251252051
17392984896084128488626945604241965285022210661186
30674427862203919494504712371378696095636437191728
74677646575739624138908658326459958133904780275900
99465764078951269468398352595709825822620522489407
72671947826848260147699090264013639443745530506820
34962524517493996514314298091906592509372216964615
15709858387410597885959772975498930161753928468138
26868386894277415599185592524595395943104997252468
08459872736446958486538367362226260991246080512438
84390451244136549762780797715691435997700129616089
44169486855584840635342207222582848864815845602850
60168427394522674676788952521385225499546667278239
86456596116354886230577456498035593634568174324095
2 ok


PAI2


2500 digits pi = :
2501 digits
31415926535897932384626433832795028841971693993751
05820974944592307816406286208998628034825342117067
98214808651328230664709384460955058223172535940812
84811174502841027019385211055596446229489549303819
64428810975665933446128475648233786783165271201909
14564856692346034861045432664821339360726024914127
37245870066063155881748815209209628292540917153643
67892590360011330530548820466521384146951941511609
43305727036575959195309218611738193261179310511854
80744623799627495673518857527248912279381830119491
29833673362440656643086021394946395224737190702179
86094370277053921717629317675238467481846766940513
20005681271452635608277857713427577896091736371787
21468440901224953430146549585371050792279689258923
54201995611212902196086403441815981362977477130996
05187072113499999983729780499510597317328160963185
95024459455346908302642522308253344685035261931188
17101000313783875288658753320838142061717766914730
35982534904287554687311595628638823537875937519577
81857780532171226806613001927876611195909216420198
93809525720106548586327886593615338182796823030195
20353018529689957736225994138912497217752834791315
15574857242454150695950829533116861727855889075098
38175463746493931925506040092770167113900984882401
28583616035637076601047101819429555961989467678374
49448255379774726847104047534646208046684259069491
29331367702898915210475216205696602405803815019351
12533824300355876402474964732639141992726042699227
96782354781636009341721641219924586315030286182974
55570674983850549458858692699569092721079750930295
53211653449872027559602364806654991198818347977535
66369807426542527862551818417574672890977772793800
08164706001614524919217321721477235014144197356854
81613611573525521334757418494684385233239073941433
34547762416862518983569485562099219222184272550254
25688767179049460165346680498862723279178608578438
38279679766814541009538837863609506800642251252051
17392984896084128488626945604241965285022210661186
30674427862203919494504712371378696095636437191728
74677646575739624138908658326459958133904780275900
99465764078951269468398352595709825822620522489407
72671947826848260147699090264013639443745530506820
34962524517493996514314298091906592509372216964615
15709858387410597885959772975498930161753928468138
26868386894277415599185592524595395943104997252468
08459872736446958486538367362226260991246080512438
84390451244136549762780797715691435997700129616089
44169486855584840635342207222582848864815845602850
60168427394522674676788952521385225499546667278239
86456596116354886230577456498035593634568174324130
4 ok



上列程式中包括了兩個同樣可以求得 Pi 的設計, PAI1 為根據標準論文發佈的式子計算而得。 PAI2 引用的式子,則為網路上德國人曾經公開過的設計,號稱也能獲得同樣的結果。經上列二程式執行、測試後的結果發現,當計算到 2500 位數時,兩者僅最後四位數有所不同, PAI1 為 4,0952 而 PAI2 為 4,1304 ,經與 uBASIC 執行後的結果比較,有更多的 5 位數不同。因此,這種尾部幾位數的誤差,純粹起因於電腦執行多次數字計算,最後引致的捨位誤差(Round off error),調整上列設定為 2510 位數後,可得正確的結果應為4,1125.............。

目前,我將 ABC Forth 數學計算系統設定成最大輸出為 10008 位數,要增加位數,更改相關設定便能辦到。我用別人廢棄的電腦來設計系統, Pentium II 0.266GHz 的 CPU ,執行出一萬位數,耗時約五分鐘。夠了,再多的位數,不是我的目的所在,所以不去浪費時間提升位數。

我用我自己設計的大數計算系統,把求解 Pi 到幾位數的計算原理澈底搞清楚,是實際執行後的了解。也把現行系統能夠處理 Pi 到幾位數的極限能力澈底搞清楚,也是實際設計出程式後的了解。為什麼選 2500 位數刊入網文?因為它是 uBASIC 系統的最高上限能力, 2500 位數就能夠拿來進行有根據的實際比較。現行的套裝軟體 Maple 或 Mathematica 也都能獲得比較數據,但是它們很貴,我買不起,它們也不提供像我這樣的源程式,只能得到不知所云、沒有意義的答案,無補於增進知識或技術,我一直當它們是腐化人心的套裝軟體,全國只須買一套來進行數據比對就夠了。

別的程式語言,早就能計算 Pi 到許多位數了,傳統的 Forth 程式語言也行,網文上可以輕鬆得到許多這方面的資料。但是,說實在的,一般程式都很難讀懂或比對,因為程式與實際的數學表示式相去太遠,除了 BASIC 或 FORTRAN 式的數學專用程式語言外,大概都是如此。所以我才設計 ABC Forth 系統的大數計算能力。實現了設計,需要夠份量的測試程式,求解 Pi 到許多位數的問題,就是最佳範例。

3. 理論依據

如果您最近光顧過 Forth 程式語言發明人─ Charles H. Moore 先生的個人網頁,您會發現他也留了一篇新近討論 Pi 的簡短網文,標題是『Thoughts about Pi』。他當寫詩一般的方式寫這篇短文,還有哲學推理般的結論,他說:根據推論,一個四度空間內的超級球形體,圓球面積應為直徑 d 的 3 次方乘以 Pi 。

我用會講話的程式語言 ABC Forth 系統,讀 Charles H. Moore 先生網頁中他所留下來的愛詩,韻味十足,還刻意去調整詩中押韻部份的腔調,簡直就像是一種特殊的享受。分享他的喜好之外,也响應他的討論,本文就討論 Pi ,用的也全是他留下來的技術所設計出來的 ABC Forth 數學計算系統。他網文中也提到,最近有人將 Pi 算到一兆(A Trillion)位數,還建議了可參考之『Pi計算技術』的網頁。我因很久以前便已見過許多這種網頁,雖未再度按址前往一覽究竟,但對計算 Pi 的技術,已有一點認識。

能寫這篇網文之前,缺乏的就是實際使用程式進行各種快速的測試。以前雖寫過許多次這種程式了,但都因為發展困難,程式寫成,得到結果之後,便作罷了,不敢再度多寫另外的程式。這一次,我設計出了極為方便的系統,可以用來快速設計出相關程式,以便深入了解問題,現在,就培養出了許多的認識。

計算 Pi 的最原始根據,可以是只需要一個非常簡單的數學式子:

π = 16 * ATAN ( 1 / 5 ) – 4 * ATAN ( 1 / 239 ) ........(1)

以及 ATAN(X) 的冪級數展開式(Power series expansion)

ATAN ( X ) = X – X ^ 3 / 3 + X ^ 5 / 5 – ........(2)

現在,我深刻的了解到,如果能夠設計出一個具有大整數計算能力的系統,使用上述最原始的計算依據就夠了,並不需要其它的理論式子。後來被發現出來的許多可用計算式子,是為了讓計算過程,能因收斂較快而得到結果,才被提出來的,在位數的精確度上,並沒有另外實質上的貢獻。換句話說,不管您用那一套數學式子來計算許多位數的Pi,都能算得出來,差別只在算得結果的快慢而已。

因為,上列的(1)與(2)式,都是恆等式,不是近似式,恆等式代表理論上絕對相等,只要您有本事一路精確的計算下去,計算結果就絕對正確,執行 PAI1 的程式,用的就是(1)與(2)式, PAI2 則使用(2)式與下列的(3)式。

π = 176 * ATAN ( 1 / 57 ) + 28 * ATAN ( 1 / 239 ) – 48 * ATAN ( 1 / 682 ) + 96 * ATAN ( 1 / 12943 ) ........(3)

那麼,本文實在沒有必要再寬列各種其它的計算式子了,甚至於最近算出五兆位數Pi者,所依據的數學式子,也沒有必要去探尋。因為,如果我將大數計算系統設計得很好,再大的數字都能計算,給我一個十兆( 10 的 14 次方)位數的數字,送得進系統,就使用本文設計出來的程式,也一樣能得到十兆位數的 Pi 。不過,若僅用上述程式處理這樣的問題,單將10自乘十兆次以得到十兆位數的超級大整數,就已經是一個相當驚人的乘法計算量了,隨後還得進行更多的計算,因此,不是輕鬆的計算。

我所設計出來的大數計算系統,並沒有計算位數的上限,現行的上限只受限於硬體,指的是 Forth 系統以上所有未被用到的記憶體容量,全部都能用來放置計算之值。算法中只需設計三個位數容量相同的大數,而每 8 個位元組(Bytes)能表示大約 21 億以上,也就是大約 10 位數。 Forth 系統的選址能力可以到達多少空間?計算過程就能存取到相同的程度。那麼十兆位數就需要三十兆乘以 0.8 的位元組量,能裝 24T Bytes 的記憶體,執行 64 位元為一單元選址能力的 Forth 系統就能辦到。 Forth 是能將類似硬碟之大批儲存媒體的實體記憶體,納入選址範圍的程式語言,這樣的 Forth 系統可以設計得出來。我設計的 ABC Forth 系統,全用 Forth 高階指令設計,可隨基礎 Forth 系統而轉殖,上列程式也就能繼續被執行。

系統的選址能力而不是計算能力,成了解決這個問題的絕對關鍵。這也就是為什麼能算五兆位數Pi的人,沒有能夠得到學術地位的原因了。以學術方式,首先提出(1)或(3)式的人,才能獲得學術名聲,還有許多人提出過各種冪級數來算 Pi ,都寫成了論文,並證明提出的式子是一個恆等式,學術界以這種方式接受研究結果,也很合情合理,重點就在『恆等』而不僅只是『幾位數』。

本文列出了實際程式與理論依據的數學式子,兩者間的對應,是直接的關係,不太需要解釋。計算Pi到小數點後面幾位數的問題,必須使用固定點算術來計算,而不可以使用浮點算術,因為浮點數表示的數字,離小數點愈遠的位數,就有愈多的棄用數字(De-normal number),它的不精確數字表示方式,算不出精確的 Pi 。固定點數字還必須當作放大了等同位數的大整數來處理,才能得到運算後的正確絕對值,因此,設計程式時,想算幾位數?就得先行將那麼大的一個超級大整數先行準備好。例如:我們想算2500位數,就得準備一個 2500 位數的大整數,它代表一個 10 的 2500 次方的超大數字,再以它為根據來進行 ATAN 函數的計算,算得出來,就能得到正確的結果。

上述程式, ATAN 的累積演算結果,就存放在 ATAN# 內,演算過程中的暫存數字,放在 W# 中,要計算冪級數多少項?才會令暫存數字 W# 內已無數可算,是一個算了才知道結果的問題,因此,把 W# 算到 0 了,就表示計算位數到達了設定的位數了。算得的 ATAN 再經過(1)或(3)式的逐項累積計算,就能存進 PAI1 或 PAI2 ,就能得到最終結果。
隨著系統的不同,將 W# 算到 0 的計算次數會有所不同,我的系統統計數字如下:

1/N........KK
5........3581
239......1057
57.......1429
682......885
12943....613

換句話說,如果 1/N = 5 ,程式便計算了 ( 3581 - 1 ) / 2 = 1790項,才讓 W# 等於 0 ,也就是算到了 X^3851 / 3851 這一項, W# 的內容才消失成 0 。

解釋到這裡,讀者應該明白,我們已經有了(1)式,為什麼還去發展出(3)式來的原因了。因為

(1)式要計算:

( 3851 - 1 + 1057 -1 ) / 2 = 2453(次)

(3)式要計算:

( 1429 - 1 + 1057 -1 + 885 - 1 + 613 - 1 ) / 2 = 1990(次)

如果要求位數愈多,計算次數的差別就會愈來愈大,所以不要小看從德國人那裡得來的資料,(3)式有它存在的價值。

還有,演算過程中使用了大量的除法,除到最後,免不了會有除不盡的捨位問題,計算次數愈多,捨位的可能次數就愈多,結果就造成了 2500 位數最後的幾位數不夠精確的問題。 PAI1 與 PAI2 的計算結果都有 173 與 179 的誤差,是算了幾千次累積出來的結果。 uBASIC 的誤差則為 17013 , uBASIC 內的計算方式不得而知,顯然計算次數至少也是幾萬次,否則不會累積出這麼大的誤差。

耗用記憶體與執行速度,均非現行電腦程式設計上必須嚴格考慮的問題,因此,我在設計程式時, 2500 位數就直接宣告 2500 ALLOT ,要求系統控留出 2500 個位元組(Bytes)的記憶體空間來,不去考慮可以乘以 0.8 以節省記憶體用量的問題。要等幾分鐘才能得到一萬位數的答案,大概是我個人可以容忍的時間耗用量,因此,我也把系統處理位數的設定,暫時固定在一萬位數,系統內就固定規劃了 10008 個位元組的印出數字緩衝區,不過,這些規矩純粹只是我的個人喜好設定,可以隨時修改。

4. 討論

上述簡單的說明,都只是針對求 Pi 到幾位數的問題而進行實測後的敘述,若非基於簡明的程式與輕易的操作性能,想體會出這些結果,並不容易。指令 SUBATAN 中,在 150 列後面新插入一列

155 PRINT KK

重新載入程式再執行,您就可以得到第 3 節中的分析數據。調整新的設定數據,也是輕而易舉、簡單、明白。這是我設計出這個大數計算系統的基本目標。

我還有後續工作必須進行,就是為這樣的系統添加可有的常用函數,例如:大數的整數方次,或大數的開平方,其它類似一般整數運算所需要的函數已經建好了。待函數建立齊全之後,就能落定 ABC Forth 數學計算系統的大數計算全套功能,算是系統性能的一大躍升。

發展這個系統最大的難題,不在數學計算自身,因為源程式是先聖先賢的創作,為了完成我自己的要求,源程式不得不進行大肆修改,發展時,主要的困難出現在如何安排出恰當的大數輸入、輸出指令。

上列程式中尚未顯示輸入功能,只顯示了一個名為 I>BIG 的小整數轉換成大整數的指令。我下了很大的功夫,以自己設計的 ABC Forth 系統追蹤出一些問題,發現如果僅藉助於S” ..... ”標準指令來安排大數的輸入,系統文字輸入緩衝區,規劃了最多只能為 260 個字數的限量,限制了輸入大數的可用數字量。直接輸入的數字,如果超過限制量,便無法正確轉換。因此,必須另行設計新的指令或新的格式來解決問題。大數的後面,必須配合另一個名為 S>BIG 的指令,它可以根據任意指定位址作為起始位址,轉換出該位址以後任意位數的數字。目前,僅暫時使用S” ..... ”標準指令設計一般程式。

如果您嫌直接將大數寫在程式內太麻煩,而想把一個超大數字明明白白的儲存在純文字檔案中,例如就儲存了 10K Bytes 大約為一萬位數的數字,想要轉換進我的大數計算系統,一點問題都沒有,依靠 S>BIG 指令就能完成。轉換的位數沒有上限,我估計,只要您單憑打字打得出來的數字,它都能辦得到。上列程式顯示了 BIG.(Big dot,大數拓印)為數字的輸出指令,它的輸出上限,基本上暫定為 10008 位數,是為了避免輸出您我都看不完的位數而故意設定的,它也可以沒有上限。

幾個典型的測試程式,就如下列,它們都是執行無誤的標準程式。

\ 大數測試程式

INTEGER I

0 BIGVARIABLE A1 2504 ALLOT
0 BIGVARIABLE B1 100 ALLOT
0 BIGVARIABLE C1 100 ALLOT

10 BIGVARIABLE BIG10
0 BIGVARIABLE DD

: TEST01 BASIC
10 LET B{ A1 = S" 12,3456,7890,1234,5678,9012,3456,7890 " S>BIG }B
20 LET B{ B1 = S" 98,7654,3210,9876,5432,1098,7654,3210 " S>BIG }B
30 LET B{ C1 = A1 + B1 }B
40 RUN CR
50 PRINT " C1=A1+B1 and C1 = : "
60 RUN C1 BIG.
70 RUN CR CR
80 END
;

: TEST02
B{{ A1 = S" 3 " S>BIG }}B
B{{ B1 = S" 5 " S>BIG }}B
B{{ C1 = A1 * B1 }}B
CR ." A1=3 , B1=5 , C1=A1*B1 , and C1 = : "
C1 BIG.
;

: TEST03 BASIC
10 LET B{ A1 = BIG0 }B
20 IF B{ A1 = BIG0 }B THEN 60
30 PRINT " A1<>0 "
40 GOTO 999
60 PRINT " A1=0 "
999 END
;

: TEST04 BASIC
10 FOR I = 1 TO 10
20 LET B{ DD = BIG0 + I>BIG ( I ) }B
30 RUN DD BIG.
40 NEXT I
50 END
;

系統還有其它許多功能,舉例而言,本文列示的程式,實際上不單只是求 Pi 到幾位數的功能而已,別忘了,程式中設計了 ATAN 函數,它就代表一個非常標準,能將 ATAN 函數精確值求到非常多位數的程式。凡是可以寫成冪級數方式表示的任何函數,都可以依照本文中設計的基本方式設計出來,這種函數有無限多個。

我曾經參與 fig 總會科學程式庫中一個公益程式的審核工作,程式希望將特殊函數之值,精確的計算到小數點後面十幾位數,該函數也能用冪級數展開式表示。那時,專注於原作者建議的多精度浮點計算系統的運用,自己尚未發展這一套大整數計算,所以是糊里糊塗的接受別人提供,以 C 程式語言寫成的的多精度浮點功能,根本不知道計算結果到底精不精確?或絕對精確到幾位數?現在自己設計出了這個大數計算系統,再度從事那樣的工作時,會有最好的比較依據了。

大整數計算的精確度,是可以辦得到絕對精確的,多精度浮點數系統則辦不到,這是兩者最大的差別所在。

多年以前,我就曾經在 eForth 系統上設計過固定點數字計算系統,也實現了全套必要基本函數的設計,有過充足的設計經驗。現在,我如果還想硬幹,再設計一套小數點後面有一千位數精確度,小數點以上也有一千位數大數容納能力的固定點數字計算系統,絕非難事。那些斤斤計較於小數點後面十幾位數還算不準的浮點計算問題,相形之下,就根本算不了什麼了。從執行這個算出 2500 位數 Pi 值的程式顯示,我的舊電腦所須耗用的時間,還在我可以容忍的範圍之內,因此,必要時,我會用這樣的固定點計算系統,去解釋別人的浮點計算結果精不精確?

大學裡教我數值分析課程的王九龍老師,教過我們,二的一百次方(2^100)在程式中要改寫成下式:

(((((2^2)^2)^5)^5)

系統才會算得更準,為什麼?

如果您手頭有科學用小型計算器,直接算算看,我所得到的答案如下:

2^100 = 1.267650599E30
(((((2^2)^2)^5)^5) = 1.2676506E30

兩個答案確實有一點不一樣,那一個比較準?據說是後面算出來的這一個。現在,我有了我自己設計出來的大數計算系統,何不自己驗證?

\ 精算2的100次方

INTEGER I
0 BIGVARIABLE A1
1 BIGVARIABLE BIG1
2 BIGVARIABLE BIG2

: 2^100 BASIC
10 LET B{ A1 = BIG1 }B
20 FOR I = 1 TO 100
30 LET B{ A1 = A1 * BIG2 }B
40 NEXT I
50 RUN CR CR
60 PRINT " A1=: "
70 RUN A1 BIG.
80 END
;

執行2^100後顯示:

A1=:
31 digits
1267650600228229401496703205376 ok

我的系統可以明確的告訴您,那一個才比較準了,確實是後面那一個。不僅如此,我還能驗證系統內軟體或硬體的精確度,到達什麼程度?如果您勤快一點,就分別去開啟任何具有浮點數功能的 Forth 系統,例如: Win32Forth 的前後兩種版本 V4.2 與 V6.14 ,然後執行下列指令

30 SIGDIGITS !
2E0 100E0 F** FE.

您將得到浮點系統告訴您之 2 的 100 次方的答案,它們都來自於硬體的數學處理器(Co-processor),再經軟體設計的輸出指令告訴您結果, V6.14 版使用的是 8 Bytes/Float ,而 V4.2 版使用的是 10 Bytes/Float ,實際輸出結果如下:

V6.14版:1.26765060022822944000000000000 E30
V4.2版: 1.26765060022822940000000000000 E30
絕對值 1267650600228229401496703205376

您請自己看看上面列示的評鑑結果,不管您未來是幾位元的系統,或者是您很自以為是的使用現成邏輯閘(Logic gate),燒製出了一顆能算出許多位元的數學處理器(Co-processor),我現在就準備好了工具,能用來明確說出系統計算數字後的絕對精確程度。

我也曾經公佈過將階乘( Factorial, ! )計算到一億階乘(100000000 !)的網文,用的是浮點數計算方式。現在,可以用大整數來算出非常大的精確答案了,這些指出來的功能,都能讓許多高等數學教科書的內容進行改寫。我不是數學家,但我設計出了研究這種問題所需要的最佳工具系統。

完成的發展證明了我在『玄數』網文中的數字資料格式觀念是正確的,所以設計程式時,才能腦筋清楚的進行規劃。同樣的模式,可以運用在基因圖譜的運算系統發展上了。生物學我尚外行,基因圖譜該如何運算?得再讀不少書才能了解。但我知道,世上還沒有能用 BASIC 程式語言來設計基因圖譜比對運算程式的系統。我想設計一個輸入 900 個 AGTC 基因圖譜後,能自動到圖譜資料庫中,比對出資料的系統來。或者是進行兩套基因圖譜的比較後,獲得有百分之幾相同的數據。甚至於生物學家能夠告訴我的運算需求,我應該都設計得出來。

我的研究發展確實是離一般群眾愈來愈遠,所以這些網文的讀者也會愈來愈少,但不管怎麼說,這些技術都還是一些有用的東西,已經寫出來公佈了,就放著吧。表面上,我的生活習慣看起來似乎很窮,使用 32 位元老掉牙的電腦發展系統,太落後了。大家都期盼著64位元的裝備普及起來,就有 64 位元的 Forth 系統可用了,也許 64 位元的裝備確實能夠提供許多難以想像的電腦性能,但就數學計算而言,今天我所完成的 ABC Forth 系統大數計算功能,已經突破了不僅只是 64 位元的整數範圍而已,沒有上限!而且可以『永遠』沒有上限!


附註 : 20241021 重新整理後貼出。

沒有留言: