2011年7月15日 星期五

艱苦環境裡的大整數精確計算

艱苦環境裡的大整數精確計算


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


相當的年紀才能有相當豐富的回憶,過往的時日還不能醉生夢死,回憶才有意義。這幾天,這裡天氣變得很冷了,冷到坐下來就能大量的回想到過去的事情,因為,不少老友都找上門來了,自然就勾起了我的許多回憶。

先是彰化的林伯龍醫師,透過此處的網文,捎來了訊息,彼此都感到非常意外。我也立刻回信,告訴他,不會忘記將近 30 年前,中華民國 Forth 語言協會草創初期,他曾經給過我們的金錢資助,並請他問候也住彰化的陳信良先生。當年,他們都捐助過協會不少錢,協會卻是一個地下組織,我拿了錢就能辦事,搞的也不是什麼偉大革命,只是推廣 Forth 。

從小,我也經歷過許多的欺壓與霸凌,但對方後來的下場都不太好。打過我的人,有的還都死掉了。我也常聽說,欺負過我的人,後來的生活都很悽涼、潦倒。可是,曾經幫助過我的人,再度與我聯絡上時,看來都有了滿意的人生,林醫師就是其中之一。我不是鬼神,也不會詛咒別人,蒙受任何打擊時,內心只會要求自己更堅強的站起來,更要充實自己,也有著想當世界天皇的偉大志向,其他甚麼天皇我都不會放在眼裡。

30 幾年前,職業軍人敢在社會上搞地下組織,是屬於違反軍紀的重大犯罪行為,嚴重的,可以撤職查辦,我卻胆大妄為到在社會上私募款項。曾經被也是同行的退伍軍人,一狀告回當時的服務單位,他們把持著國營公司主管的權位,不捐錢也就罷了,還陷害忠良,差一點就令我落得與郭冠英同樣的下場。幸好,承辦告案的魏元勳副所長深明大義,放我一馬,對我算是有著不殺之恩。後來,他出題考我,能不能用 Forth 把圓周率 π 算到小數點後面 1000 位數?我算不出來,於是 Forth 被判定成沒有用。我終身記得此事,便立定志向要自己設法解決,雖然多年以前,我就能解這個問題,但我要等到電腦軟硬體的發展都足夠成熟之後,才願意實際使用 Forth 解決問題。能夠解決,才值得在網文中公佈,算是了了一股心願,從此之後,也才可以不必再談 π 的事情。

我有終身反省的習慣, 20090515 刊出的:『FORTH解大數問題』談到 π ,那時,我用一般 Forth 程式解問題,可以獲得正確答案,但是沒有特別意義。前一篇網文公佈了我想將 π 算到幾位數都可以了,就是為了回報魏副所長的不殺之恩,用的不僅只是實際程式,還是我親自設計出來的程式語言系統。以前,我若堅持已見,想做這種發展,會被視同為不務正業,結果也可能會是撤職查辦。我的本業被上級視同為掃放射性的工作,而且全單位都沒人愛做,我就只好忍氣吞聲的一直幹著吃放射性的工作,難怪七早八早就被國軍淘汰了,也不夠資格領終身俸,但沒有怨言,現在只能搞 Forth 。

我還有另外兩個必須反省之後,透過網文向讀者交待的問題,但不必計較會有多少讀者來點閱。 20100116 刊出的網文:『高級數學函數』,文內提到,當年一月八日出刊的『科學美國人』(Scientific American)裡面的一篇文章,論及一個十進制為232位數的密碼性數字,可以因式分解成兩個均為 116 位數的質數。該文中,我寫下誓言,『它日我必會用到這份資料驗證程式』,現在,實現這個諾言,程式要用大整數除法算對結果才有價值,因為,大數的加、減、乘法都很容易設計,除法不容易。

列示這個程式,代表我對所有研究問題認真的態度,不能設計出算對結果的系統前,我長期被大數無法正確輸出與輸入系統的問題所困擾,請教國際網站的負責人,也得不到答覆,他們反而要我幫忙做另外一件事情。最後,只好自行解決問題,完成這件工作後,令我憶起,從小都是這樣子自立自強長大的,我經常鼓勵中國同胞努力站起來,我自己更應該努力站起來。

0 BIGVARIABLE A1 1000 ALLOT
0 BIGVARIABLE B1 1000 ALLOT
0 BIGVARIABLE C1 1000 ALLOT

: MAIN BASIC
10 LET B{ A1 =
S" 1230186684,5301177551,3049495838,4962720772,8535695953" FIRST-SECTION
S" 3479219732,2452151726,4005072636,5751874520,2199786469" NEXT-SECTION
S" 3899564749,4277406384,5925192557,3263034537,3154826850" NEXT-SECTION
S" 7917026122,1429134616,7042921431,1602221240,4792747377" NEXT-SECTION
S" 9408066535,1419597459,8569021434,13"
NEXT-SECTION
S>BIG }B
20 LET B{ B1 =
S" 3347807169,8956898786,0441698482,1269081770,4794983713" FIRST-SECTION
S" 7685689124,3138898288,3793878002,2876147116,5253174308" NEXT-SECTION
S" 7737814467,999489"
NEXT-SECTION
S>BIG }B
30 LET B{ C1 = A1 / B1 }B
40 PRINT " A1 = : "
50 RUN A1 BIG.
60 PRINT " B1 = : "
70 RUN B1 BIG.
80 PRINT " C1=A1/B1= : "
90 RUN C1 BIG.
100 END
;

執行結果如下:

MAIN
A1 = :
232 digits
12301866845301177551304949583849627207728535695953
34792197322452151726400507263657518745202199786469
38995647494277406384592519255732630345373154826850
79170261221429134616704292143116022212404792747377
94080665351419597459856902143413
B1 = :
116 digits
33478071698956898786044169848212690817704794983713
76856891243138898288379387800228761471165253174308
7737814467999489
C1=A1/B1= :
116 digits
36746043666799590428244633799627952632279158164343
08764267603228381573966651127923337341714339681027
0092798736308917 ok

驗證無誤。

另一件必須向讀者交待的事情,是有關於算出超大質數(Prime number)之事,我在 20090210 『質數』網文中,提供過算出質數之一般 Forth 程式,它並不是什麼偉大的創作,只適合用來當作小型範例程式,當時(2008-09-18)梅森質數大搜索組織(GIMPS, Great Internet Mersenne Prime Search)活動的結果,已經算出了 2 的 43112609 次方減 1 為新發現的大質數,此數大於 1200 萬位數(10^12978188),現在該組織可能已經發現更大的數字了。

我重提此事的目的,是想告訴大家,尤其是全體中國人,找梅森質數的數學計算方法非常簡單,但一般程式語言系統不易設計出那樣的程式,必須是擁有無窮位數大整數計算功能的程式語言系統,才能用來輕易的設計出程式。我不想批評別種程式語言在這方面的能力,只問自己設計出來的 ABC Forth 數學計算系統,又能做到什麼程度?現在則是展示本領的時候了。

因為我已設計出了輕巧的工具,可以實現這套尋找大質數方法的程式,並快速執行後,得到測試結果。也因此能通順、詳細的解釋這套技術。我們先行詳細解釋一些相關術語。

任何數只要可以表示為 M(P) = 2 ^ P – 1 ,就被稱為梅森數(Mersenne number,梅森是一位法國數學家,1588~1648),如果 P 不是質數,則 M(P) 也不會是質數,若 M(P) 為質數,則 P 必定為質數,但是, P 為質數, M(P) 卻不一定為質數。

梅森當初提出上述研究時,宣稱:當 P = 2,3,5,7,13,17,19,31,67,127,257 等 11 個質數時, M(P) 都是質數,而除了這 11 個質數以外,比 257 小的其他質數 P ,計算出來的 M(P) 都不是質數。

可是,後人發現 M(67) 不是質數, M(257) 也不是質數,但 M(61)、M(89)、M(107) 反而都是質數。

當 P 為較大的質數時,能夠計算出來的 M(P) 就會非常的大,因此,就必須經過更大數字的大量計算測試,才能確定該 M(P) 是否為質數?早期沒有電腦可用的時代,這種工作,實在難以實現。後來,法國數學家盧卡斯(Edouard Lucas),於 1856 年與 1878 年,分別提出過發現了用來測試 M(P) 是否為質數的方法, 1930 年代,萊默(D. N. Lehmer)又研究出改良的測試方法,這種測試方法,後來就被稱為盧卡斯萊默測試法(Lucas Lehmer test, LLT),可以高效率測試出梅森數是否為一質數?目前,此法是被普遍用來發現梅森數中更大質數的最佳方法。國際網路梅森質數大搜索組織,就是專門從事於這種活動的一個組織。

盧卡斯萊默測試法,可以用下列程式語言式的敘述表達清楚,意義能夠直讀,我就不用詳細解釋。

\ Determine if M ( P ) = 2 ^ P – 1 is a prime?
\ Lucas-Lehmer test
Lucas-Lehmer( P )
var S = 4
var M = 2 ^ P – 1
repeat P – 2 times:
S = ( ( S * S ) – 2 ) MOD M
if S = 0 return PRIME else return COMPOSITE

但是,讀者請注意,敘述中的 S 及 M 都必須能夠被處理成無限位數的大整數,上列程式才有效, P 則暫時還可以被處理成 21 億以下的一般單整數。不是所有現成的程式語言都能辦到此事,能辦到的大概也不會是免費的軟體。我改寫出一個精簡交談式的 ABC Forth 程式,簡單到就如下列所示,重點顯現在可以快速測試出能夠告訴大家特性的實際數字結果。

INTEGER I
INTEGER P

0 BIGVARIABLE S 10008 ALLOT
0 BIGVARIABLE M 10008 ALLOT

0 BIGVARIABLE BIG0
1 BIGVARIABLE BIG1
2 BIGVARIABLE BIG2

: LUCAS BASIC
10 PRINT " Input P for (2^P-1), but P must be > 2 !!! "
20 INPUTI P
30 RUN CR CR
40 LET B{ S = I>BIG ( 4 ) }B
50 LET B{ M = BIG2 ^ P - BIG1 }B
70 FOR I = 1 TO P - 2
80 LET B{ S = ( ( S * S ) - BIG2 ) MOD M }B
90 NEXT I
100 IF B{ S = BIG0 }B THEN 130
110 PRINT " For P = " ; P ; " : (2^P-1) is a composite. "
120 GOTO 150
130 PRINT " For P = " ; P ; " : (2^P-1) is a prime. It is: "
140 RUN M BIG.
150 RUN CR CR
160 END
;

程式再經簡單的改寫、調整,就能自動的執行從所有梅森數中找出梅森質數的工作,以下為只計算到 2000 位數時的印出結果:

For P = 3 : (2^P-1) is a prime. It is:
1 digits
7
For P = 5 : (2^P-1) is a prime. It is:
2 digits
31
For P = 7 : (2^P-1) is a prime. It is:
3 digits
127
For P = 13 : (2^P-1) is a prime. It is:
4 digits
8191
For P = 17 : (2^P-1) is a prime. It is:
6 digits
131071
For P = 19 : (2^P-1) is a prime. It is:
6 digits
524287
For P = 31 : (2^P-1) is a prime. It is:
10 digits
2147483647
For P = 61 : (2^P-1) is a prime. It is:
19 digits
2305843009213693951
For P = 89 : (2^P-1) is a prime. It is:
27 digits
618970019642690137449562111
For P = 107 : (2^P-1) is a prime. It is:
33 digits
162259276829213363391578010288127
For P = 127 : (2^P-1) is a prime. It is:
39 digits
170141183460469231731687303715884105727
For P = 521 : (2^P-1) is a prime. It is:
157 digits
68647976601306097149819007990813932172694353001433
05409394463459185543183397656052122559640661454554
97729631139148085803712198799971664381257402829111
5057151
For P = 607 : (2^P-1) is a prime. It is:
183 digits
53113799281676709868958820655246862732959311772703
19231994441382004035598608522427391625022652292856
68889329486246501015346579337652707239409519978766
587351943831270835393219031728127
For P = 1279 : (2^P-1) is a prime. It is:
386 digits
10407932194664399081925240327364085538615262247266
70480531911235040360805967336029801223944173232418
48424216139542810077913835662483234649081399066056
77320762924129509389220345773183349661583550472959
42054768981121169367714754847886696250138443826029
17323488853111608285384165850282556046662248318909
18801847068222203140521026698435488732958028878050
869736186900714720710555703168729087

我以不建立質數表的方式,讓上列程式就直接逐次加1的往上計算到 2000 為止,所以可以確定其間沒有其它的質數是梅森數了。我寫出程式的目的,不在與前人比高下,只想用程式證明,除了上列數字外,確實沒有其他的梅森質數了。我沒有用我的電腦繼續算下去,只算出了刊出這篇文章所需要的數字就夠了,讀者大約很難在其他的場合,見到一個名符其實的 386 位數的真正質數了。我將系統暫時設定成可以得到一個萬位數以下的質數,因此,就目前全世界已經公開的 47 個梅森質數而言,我的系統不改設定,可以算得第 32 個以下的梅森質數。修改設定,可以提高記錄,但我只有艱苦環境裡的電腦,不必與有錢人爭長短。我的電腦在這個程式執行期間, CPU 的冷卻風扇不斷的起動運轉,這是紐西蘭冬日裡,環境溫度接近 0 度 C 時的罕見現象,可見它們跑得有多累!

這些數據都是有用的大數據,用來測試系統時可以派得上用場。我從科學美國人期刊上取得 116 位數的質數,就能檢驗出系統的輸出輸入,會受到微軟作業系統內不良設計的嚴重影響。大家都花了那麼多錢向微軟購買軟體,受到這種氣,卻不能告它來索賠,只能透過這種技術網文,將問題披露出來。我也不是光會在社會上取得資源而不回饋的人,上述資料,就回饋給全世界,比科學美國人期刊能夠提供的數字大三倍位數的質數,更大位數的資料,我也能夠產生。

我在 20090210 :『質數』網文中,詳細介紹的程式及其執行結果,對找出超大質數的研究也有幫助。如果將那些質數建成檔案,再利用本文內的程式,兩者配合起來執行,就能更快速的找出更大的梅森質數,不必以每次加1的方式,逐次的往上計算。兩者配合起來的找法,可以輕易的執行到2的二十一億次方的計算,可以超過第47個梅森質數的限量,但我對安排硬體來執行可能可以打破記錄的事情,沒有興趣。

談完了技術,可以話點家常。

前兩天,一對老朋友夫婦來訪,順便找我協修一些修不好的電器,肉眼能夠瞧出來的問題,大概都難不倒我,所以東西很快就都修好了。

第一個問題電器,是被拆到螺絲頭完全毀損、變型而無法受力之深槽內的螺絲,老友確定已經拆不下來了,所以帶來找我協修。我找遍了自己花了不少錢買來的所有工具,確實也沒有適當的螺絲起子,可以派得上用場,因為那根螺絲是一個故意將螺絲頭設計成三線凹槽式的螺絲,目的就在防止一般使用者,任意自行拆開電器來維修,發生觸電意外。

後來,我突然想到,一般人動手工作的方式與我不同,我因經常自己動手做事,知道一些工作技巧,也許那根螺絲換我直接拆,就能拆得下來,因此,就取用了一根1字型的長起子直接拆卸,果然不錯,深槽內螺絲頭全毀成圓洞的螺絲,我直接拆出來了,電路斷掉了的電器,當然也就修好了。

第二個問題電器,是一個新型獨立的液晶顯示幕,不能正常工作,系統告知傳輸線有一條是斷訊狀態,老友不敢擅拆機背板來查線,所以拿來找我協修。我只看了看整組機子的情況,就告訴他,這種東西敢出廠銷售,多股訊號傳輸線發生機內斷線的可能機率相當低,反倒是訊號傳輸接頭又接了轉換接頭,出問題的可能性較高,不如先拆驗接頭。直接拆開一看,又是果然不錯, D 型多孔接頭的一個凸針,根本沒插入對應插孔,全被壓歪了。以小起子併同尖嘴鉗,小心扳回銅針後,看來還能使用,修好了這個大概並不便宜我都還沒買過的液晶顯示幕。

做這些事情,會令我想起自己曾經有過的艱苦人生。

我拆螺絲時,知道扳轉起子的第一時間,必須壓得夠緊,並只施以小角度的衝擊式扭力,感到螺絲確實鬆開時,必須立刻解除壓緊之力,改成單純的只施以扭轉之力,特別是金屬螺絲安裝於塑膠機殼的情況尤然,否則螺絲要不是拆壞螺頭都還拆不下來,就是發生滑牙現象,從此之後都難再拆裝了。英文有句俗語,能幫助經常搞不清楚扭轉方向的人,記住『旋緊』或『扭鬆』的方向,他們說:『Right ~y tight ~y, Left ~y loose』,~y要拖長了發[i]的音,是『右旋變緊,左旋變鬆』的意思。

我的表達能力很好,這一段話是標準起子使用技巧的經典說明,如果您能讀懂這一段話,終身牢記,身體力行,保證您受用不盡。至於話語的意義,您應自己體會,我大學學的是機械工程,是一名具有工程師執照的淘汰老兵,當然能夠詳細解釋出道理,但我會告訴想學的人,就像學 Forth 程式語言一樣,您只宜拿著起子直接動手照著做,多做熟練之後,手腦才能配合得起來,光知道道理沒有用。大家同時拿著起子拆螺絲,同做這麼簡單的一件工作,為什麼會拆出不一樣的結果?因為我出身卑微,有過艱苦的人生。

從小到井邊打水、挑水回家、生火煮飯,其實都還稱不上艱苦,但現代人都不必做了。

我沒上學之前,就陪著家人到海邊防風林撿柴火,當個小樵童,經常要躲避森林巡守員的為難,這也不算什麼,但聽到穿林的松風,環顧四下無人的荒野森林時,心中卻油然生起強烈的意志,立志將來要讓家人不必這樣過日子。

我年紀雖小,找柴迅速,每次都是我陪著媽媽最先回到產業道路旁約好的集合地點,我們要先藏好媽媽的一擔柴,與我能背的一竹簍子柴,等了一陣子還不見鄰居來會合,我媽媽就會要我去路的盡頭把風,她才高喊:

『花~ ~ 草嬸媽哦 ~ ~,好轉去了 ~ ~ 』(要用台語喊,是要花草嬸回去的意思)

這樣的吆喝,一輩子在我腦海裡迴盪。我永遠記得,來到這個世界,媽媽把我養育長大,我負有重任,要改善我們的環境,要改善我們的生活,一生遭受任何打擊,都不得低頭,必須昂首闊步,堅強前行。

從小生活環境不好,經常要倒馬桶,有夠粗鄙吧?那又如何?夜晚穿過黑漆漆的暗巷,端著馬桶,前往公廁後面的糞便池,您以為我心裡能想甚麼?我只想到將來絕不讓我的家人,還要這樣過日子。

我14歲時就上過大夜班,在木器工廠當童工,賺錢貼補家用,清晨四點鐘,下班騎自行車回家,一路上只有我孤獨的一個人,也是黑漆漆的,只有星光照著路面,路過無人看守剛撞死人的鐵路平交道,那時,屍體都只能放在現場好幾天,只蓋上草蓆,您以為我這個童工這麼勇敢,什麼都不怕?但是我從不跟家人提起這些事情,說出來就表示自己將來一定沒出息。

大學畢業時,母親過世了,家人生活還是很苦,我沒能盡任何孝道,我媽媽也從來不曾要求我要回報她。因此,改善環境、改善生活的志向,轉而奉獻社會,這絕對就是我爸爸、媽媽的意思。

有過這些卑下的人生,才會顯示出我拆螺絲與眾不同之處,通常光憑肉眼看不出來。就算成年踏入社會了,我也並非就能過著好日子,一切仍然都要靠著自己,不努力就不可能有好的結果。前面說過,吃公家飯時,我盡是幹著別人不愛幹的掃放射性工作,那是打從進入公家單位開始,就長期被人欺負的結果,後來毅然決然的離職了,還脫離不了不良的影響。

我估計台灣當今也沒幾個能把放射性測量工作做得比我確實的人,十五年前,移民之前,我自購測量放射性用的蓋革─米勒計數器(Geiger-Muller counter),自行測量家中器具有無輻射?以便決定是否搬到紐西蘭繼續使用。我熟知儀器的性能,就算沒有校準儀器的機會,我也還能量得出有意義的數據。只要儀器能夠反應,您將測量時間加長,精確程度就可增加,這是不變的真理,讀數的意義就靠比較,也能凸顯出來。最後,量出了床上長期使用的兩個枕頭,放射性計量讀數,是環境平均背景讀數的三倍,我知道輻射不算太強,等同於大家每次搭飛機,飛上高空後所必須接受的輻射強度,那時,搭飛機時,沒有不得使用任何儀器的管制,我用同一儀器在搭飛機時測量了幾小時,確實獲得這樣的數據,所以也沒必要將枕頭當輻射物質處理,我太太氣得當場就將兩個枕頭當垃圾丟了。長期在原子爐內工作,不可能每天下班前洗完澡才回家,頭髮就像一個標準的空氣過濾器,什麼東西都能吸附,長期工作下來,能不帶點放射性物質回家?

同樣的,長期使用電腦執行數學計算程式,在現今全世界的數學計算軟體環境中,我怎麼不知道誰會欺負誰?我在學習原子爐爐心設計工程時,是拉計算尺過日子的,那時,買一台進口的小型計算器(Calculator)要幾萬台幣,教我這門課程的鄭季康博士,有一台公家買的裝備,我每星期寫好50張單光紙,留下空格,標示程序,就等著鄭老師來教課時借用他的計算器,利用下課時間,一口氣填完表格,晚上再拉計算尺驗證。這是大四下學期的事情,四下的成績不列入畢業後分發工作的參考,全班同學都不想學太囉嗦的東西了,只剩我一個人努力完成原子爐爐心的工程設計,只有我交了報告,裡面全是大量的數學計算。

今天,中國人的環境,可以說電腦已經普及了,可是,有誰還會在乎與電腦相關的數學計算?大家用的是花大錢買來的鴉片式軟體,吃了就上瘾了,誰還想要自行設計?當老師的,選擇軟體之依據只有一個,看看是否有利可圖?管它什麼公益軟體?我聽到台灣開 Forth 月會時,講員吐出來的心聲,學生拿到文憑進入社會工作後,沒人想用老師強教、強用的鴉片式軟體,要不然就盜用。還想正式的用,就得花大錢,而且愈來愈不知道裡面在算些什麼?大家還是希望能夠見到公益性軟體,可是,去那裡找?請問?我怎麼不知道誰在欺負誰!

最近,寫完書後,開始貼一些比較特殊的數學計算網文,有幾天,見到我的網文點閱記錄顯示,德國人的點閱數,比台灣的國人同胞還要高,大家認為我的心裡應該有什麼感想?德國人全國國民平均素質,優於其它國家,有許多技術冠蓋全世界,他們可能比台灣人更了解,在 Forth 這麼艱苦的環境中,一個人單打獨鬥,產生出來的創作,比什麼都可貴,所以願意看台灣人不願意看的東西,是中文也無妨。

不單純只是我用起子與別人用起子的方式有所不同,也不單純只是我用電腦與別人用電腦的方式有所不同,更不單純只是我用 Forth 與別人用 Forth 的方式有所不同。我過過苦日子,從艱苦環境中長大,意志與人不同,心中有點分寸,我寫的全是通順的中文,要想增進全人類的福祉,也是從自己同胞的生活環境作為開始呀!


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

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 重新整理後貼出。