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

1 則留言:

Unknown 提到...

拜讀了網文,感佩! holi