一百個例題 (46 ~ 50)
Ching-Tang Tseng
Hamilton, New Zealand
12 September 2024
FORTH:慶的世界。 本網頁為中華民國(Republic Of China)FORTH語言協會 永久性網頁
2024年9月12日 星期四
一百個例題 (46 ~ 50)
第(46)個範例於 20150412，以英文寫成文章，貼出於此個人網頁上，文章的標題為 “Curve fitting demo in ABC FORTH”。
源程式的出處為 Lawrence J. Mayes ELECTRONIC ENGINEERING, 57, No 698; FEBRUARY 1985 p 38，原文標題為 “Lagrange interpolating polynomial curve fitting“，我在英文貼文及這個範例程式中均標明來源，絕不掠人之美，只用它來考驗我所設計的系統，找出問題，詳列差異，仔細解釋，永久留參。
我在學習數值分析課程時，採用的課本為 1969 年出版的 Applied Numerical Methods ，作者有三人，兩位是密西根大學的 Brice Carnahan, James O. Wilkes 及德州大學的 H. A. Luther ，書本很大很厚。第一章就教插值與近似技術(interpolation and Approximation)。主要內容就是講解這個範例程式中所用到的除差分法(divided-differences method)求解出多項式係數的技術。
上列電子工程期刊是中山科學研究院圖書館的館藏書籍，我影印該文自用。早年期刊上的論文，品質都非常好，條理分明、言簡意賅，使用當年廣被群眾愛用的 BASIC 程式語言發展而成。文末附圖，但不提供繪圖程式，只強調插值技術的所得結果，近似度令人滿意、可用。我設計的 ABC Forth 創作，就是 BASIC 式語法，但結構更為簡單，文內所實現的程式，就是證明 BASIC 能怎麼用 ? 我創作的 ABC Forth 也能照樣使用。
由於範例程式中與我的個人貼文均有詳細的解釋，技術方面，我不擬多作強調，有興趣，就請詳讀資料。刊文所貼附圖，也是我使用 ABC Forth 系統所繪製，這部份技術，超出範例能夠涵蓋的討論範圍，是另一項龐大的主題，只能日後視情況討論。
:
源程式的出處為 Lawrence J. Mayes ELECTRONIC ENGINEERING, 57, No 698; FEBRUARY 1985 p 38，原文標題為 “Lagrange interpolating polynomial curve fitting“，我在英文貼文及這個範例程式中均標明來源，絕不掠人之美，只用它來考驗我所設計的系統，找出問題，詳列差異，仔細解釋，永久留參。
我在學習數值分析課程時，採用的課本為 1969 年出版的 Applied Numerical Methods ，作者有三人，兩位是密西根大學的 Brice Carnahan, James O. Wilkes 及德州大學的 H. A. Luther ，書本很大很厚。第一章就教插值與近似技術(interpolation and Approximation)。主要內容就是講解這個範例程式中所用到的除差分法(divided-differences method)求解出多項式係數的技術。
上列電子工程期刊是中山科學研究院圖書館的館藏書籍，我影印該文自用。早年期刊上的論文，品質都非常好，條理分明、言簡意賅，使用當年廣被群眾愛用的 BASIC 程式語言發展而成。文末附圖，但不提供繪圖程式，只強調插值技術的所得結果，近似度令人滿意、可用。我設計的 ABC Forth 創作，就是 BASIC 式語法，但結構更為簡單，文內所實現的程式，就是證明 BASIC 能怎麼用 ? 我創作的 ABC Forth 也能照樣使用。
由於範例程式中與我的個人貼文均有詳細的解釋，技術方面，我不擬多作強調，有興趣，就請詳讀資料。刊文所貼附圖，也是我使用 ABC Forth 系統所繪製，這部份技術，超出範例能夠涵蓋的討論範圍，是另一項龐大的主題，只能日後視情況討論。
:
\ (46)Lagrange插值多項式曲線調適.f
\ 原始設計日期：2012-08-06
\ 原本設計ABC FORTH系統的陣列與矩陣時，並未規劃指標為0時的儲存方式。
\ 規劃完成後，下列程式方得以執行。
\ ABC Forth demo code: Lagrange interpolating polynomial curve fitting.
\ Author: Ching-Tang Tseng, Hamilton, New Zealand
\ Date: 2015-04-09
\ Contact: ilikeforth@gmail.com
\ Reference:
\ Lawrence J. Mayes ELECTRONIC ENGINEERING, 57, No 698; FEBRUARY 1985 p 38.
\ Description:
\ This program uses a method based on the Lagrange interpolating polynomial.
\ ZZ(i) is the input raw data array.
\ Matrix AAA(i,j) is the difference table manipulation matrix.
\ HH(i) is the output coefficients array of the fitted polynomial.
5 INTEGERS I J N N1 P
: INPUT-ORDER BASIC
10 PRINT " Order of fit: "
20 INPUTI P
30 IF P > 56 THEN 50
40 GOTO 70
50 PRINT " Order P can not > 56, Re-enter please. "
60 GOTO -10
70 END ;
56 ARRAY ZZ
: SUB245 BASIC
10 RUN INPUT-ORDER
20 PRINT " Input data: "
30 PRINT " =========== "
40 FOR I = 0 TO P \ 指標從0開始
50 PRINT " Y value for X abscissa value: " ; I
60 INPUTR ZZ ( I )
70 PRINT " Y( " ; I ; " )=" ; { ZZ ( I ) } \ 注意前一I為整數
80 NEXT I
90 END ;
57 ARRAY CC
: SUB235 BASIC
10 FOR I = 0 TO P + 1 \ 指標可為算式
20 LET { CC ( I ) = 0 }
30 NEXT I
40 END ;
56 ARRAY HH
57 56 MATRIX AAA
: SUB275 BASIC
10 PRINT " Polynomial fit: "
20 PRINT " =============== "
30 FOR I = 0 TO P
40 LET { HH ( I ) = 0 }
50 NEXT I
60 FOR J = 0 TO P
70 FOR I = 0 TO P
80 LET { HH ( I ) = ZZ ( J ) * AAA ( I J ) + HH ( I ) }
90 NEXT I
100 NEXT J
110 FOR I = 0 TO P
120 PRINT " HH( " ; I ; " )= " ; { HH ( I ) }
130 NEXT I
140 END ;
INTEGER I-1
: SUB220 BASIC
10 FOR I = P TO 1 STEP -1
20 LET I-1 = I - 1
30 LET { CC ( I-1 ) = CC ( I ) * I>R ( J ) + CC ( I-1 ) }
40 NEXT I
50 END ;
INTEGER Q
REAL B0
: MAIN BASIC
10 RUN SUB245
20 FOR J = 0 TO P
30 RUN SUB235
40 LET { CC ( 0 ) = 1 }
50 FOR N = 0 TO P
60 LET N1 = N + 1
70 IF N = J THEN 90
80 GOTO 100
90 LET N = N + 1
100 IF N > P THEN 200
105 LET Q = NEGATE N
110 LET { B0 = I>R ( Q ) }
120 FOR I = 1 TO N1
130 LET I-1 = I - 1
\ 下式實際上就是差值 CC(i-1)-CC(i)
140 LET { AAA ( I J ) = B0 * CC ( I ) + CC ( I-1 ) }
150 NEXT I
160 LET { AAA ( 0 J ) = B0 * CC ( 0 ) }
170 FOR I = 0 TO N1
180 LET { CC ( I ) = AAA ( I J ) }
190 NEXT I
200 NEXT N
210 RUN SUB220
220 FOR I = 0 TO P
230 LET { AAA ( I J ) = AAA ( I J ) / CC ( 0 ) }
240 NEXT I
250 NEXT J
260 RUN SUB275
270 PRINT " Finished."
280 END ;
\ 使用例：
\ 先執行下列指令 GIVEN 以便獲得可用之輸入數據
\ 再執行 MAIN 然後於對答時輸入由 GIVEN所獲得的數據
: GIVEN BASIC
10 RUN CR CR
20 FOR I = 0 TO 5
30 PRINT " Y( " ; I ; " )= " ; { EXP ( I>R ( I ) ) }
40 NEXT I
50 RUN CR CR
60 END ;
\S
參考文獻及執行結果:
***********************************************************
100 COM A(57,56), C(57), H(56), Z(56)
105 INTEGER I, J, N, N1, P
110 CLEAR
115 LINPUT "Order of fit:",P$ : P = VAL(P$) : IF P > 56 THEN BEEP : GOTO 115
120 GOSUB 245
125 FOR J = 0 TO P
130 GOSUB 235
135 C(0) = 1
140 FOR N = 0 TO P
145 N1 = N + 1
150 IF N = J THEN N = N + 1
155 IF N > P THEN 190
160 B0 = -N
165 FOR I = 1 TO N1
170 A(I,J) = B0 * C(I) + C(I-1)
175 NEXT I
180 A(0,J) = B0 * C(0)
185 FOR I = 0 TO N1 : C(I) = A(I,J) : NEXT I
190 NEXT N
195 GOSUB 220
200 FOR I = 0 TO P : A(I,J) = A(I,J)/C(0) : NEXT I
205 NEXT J
210 GOSUB 275
215 BEEP 100,190 : DISP "finished" : STOP
220 FOR I = P TO 1 STEP -1
225 C(I-1) = C(I) * J+C(I-1) : NEXT I
230 RETURN
235 FOR I = 0 TO P + 1 : C(I) = 0 : NEXT I
240 RETURN
245 PRINT : PRINT "INPUT DATA" : PRINT "==========" : PRINT
250 FOR I = 0 TO P
255 DISP : DISP "Y value for abscissa value:"; I : INPUT Z(I)
260 PRINT "Z("; I; ")= "; Z(I)
265 NEXT I
270 RETURN
275 PRINT : PRINT "POLYNOMIAL FIT" : PRINT "==============" : PRINT
280 FOR I = 0 TO P : H(I) = 0 : NEXT I
285 FOR J = 0 TO P
290 FOR I = 0 TO P
295 H(I) = Z(J) * A(I,J) + H(I)
300 NEXT I
305 NEXT J
310 FOR I = 0 TO P : PRINT "H("; I; ")= "; H(I) : NEXT I
315 RETURN
320 END
**************************************************
我的執行結果如下
GIVEN
Y( 0 )= 1.00000000000
Y( 1 )= 2.71828182846
Y( 2 )= 7.38905609893
Y( 3 )= 20.0855369232
Y( 4 )= 54.5981500331
Y( 5 )= 148.413159103
ok
MAIN
Order of fit: ? 5
Input data:
===========
Y value for X abscissa value: 0 ? 1.00000000000
Y( 0 )= 1.00000000000
Y value for X abscissa value: 1 ? 2.71828182846
Y( 1 )= 2.71828182846
Y value for X abscissa value: 2 ? 7.38905609893
Y( 2 )= 7.38905609893
Y value for X abscissa value: 3 ? 20.0855369232
Y( 3 )= 20.0855369232
Y value for X abscissa value: 4 ? 54.5981500331
Y( 4 )= 54.5981500331
Y value for X abscissa value: 5 ? 148.413159103
Y( 5 )= 148.413159103
Polynomial fit:
===============
HH( 0 )= 1.00000000000
HH( 1 )= 2.74952933754
HH( 2 )= -3.30606647675
HH( 3 )= 3.03499879102
HH( 4 )= -.885001709373
HH( 5 )= .124821886021
Finished. ok
驗證無誤
**************************************************
原始資料輸出結果如下
H0 = 1
H1 = 2.74952
H2 = -3.30606
H3 = 3.03500
H4 = -.885002
H5 = .124822
**************************************************
(1)輸入之 X(i) 值固定為整數值 0 ... n，Y(i)值則需人工輸入，但可以為浮點數。
(2)程式中不存 X(i) 之值。只存 Y(i)之值。存於 ZZ(i) 陣列中。
(3)本程式中使用之 Y(i) 輸入範例為 exp(X(i)) 之值，為方便起見，加建了GIVEN指令。
(4)GIVEN執行後，可直接提供 exp(0) ... exp(5) 六個 Y(i) 值，做為對答時輸入之用。
(5)執行 main 指令後得到的 6 個 HH(0) ... HH(5) 係數，即為所求多項式的係數答案。
(6)本程式的缺點是輸入數據只能為整數的固定差值 0, 1, 2, 3 .... n ，不能為其他值。
(7)程式應用範圍被限制於最大為 56，因為陣列 ZZ(56) 及矩陣 AAA(57,56)
被設定成只能操作於 56 個單元的上限量。
(8)就建立差分表(difference table)應用於求取內插值(interpolation)問題的數值分析
技術研究而言，這是一個很直接便能了解之在列表內(在矩陣 AAA(i,j)內)操作計算的
良好範例。
(9)本程式設計成只處理等距數據(equally spaced data)，若欲處理不均勻差量數據
(nonuniformly spaced data)，可修改原始資料內的程式160列內容來達到目的。
(10)2015年時，這方面的技術已發展得相當完善，網上已有許多公益程式可以參考使用，
故不對此程式進行任何修改，僅就現況使用或研究便可。
(11)插值技術的應用，主要是強調內插(interpolation)而非外插(extrapolation)，因此，
插值範圍以外的數值不應以內插結果來求取。
(12本程式的輸出結果，可以利用函數繪圖程式予以圖形表示效果。
(13)精確度設定為 6 位數時的執行結果如下：
given
Y( 0 )= 1.00000
Y( 1 )= 2.71828
Y( 2 )= 7.38906
Y( 3 )= 20.0855
Y( 4 )= 54.5982
Y( 5 )= 148.413
ok
main
Order of fit: ? 5
Input data:
===========
Y value for X abscissa value: 0 ? 1
Y( 0 )= 1.00000
Y value for X abscissa value: 1 ? 2.71828
Y( 1 )= 2.71828
Y value for X abscissa value: 2 ? 7.38906
Y( 2 )= 7.38906
Y value for X abscissa value: 3 ? 20.0855
Y( 3 )= 20.0855
Y value for X abscissa value: 4 ? 54.5982
Y( 4 )= 54.5982
Y value for X abscissa value: 5 ? 148.413
Y( 5 )= 148.413
Polynomial fit:
===============
HH( 0 )= 1.00000
HH( 1 )= 2.74928
HH( 2 )= -3.30559
HH( 3 )= 3.03469
HH( 4 )= -.884923
HH( 5 )= .124815
Finished. ok
第(47)個範例展示一個 Forth 與眾不同的性能，能夠直接操控硬體的能力。
繼前一個相當複雜的數學計算應用範例之後，今天介紹的範例，就顯得比較容易了解。
Forth 的能力，是能簡、能繁。繁的部分，改善出 BASIC 語法格式的程式設計能力，都能辦到。簡的部分，是 Forth 程式語言的特性，長久用、時常用、用久了就能體會，不用就難體會。
這些百例，非一天的收集而成，我在多年收集之後，每次回顧內容，都會產生新的感受，手癢到很想再做點變化，重新再試一試。也覺得百例仍能擴充出千例，只是願不願繼續做下去的問題，因為我粗略估算過，幾十年來，我寫成過的程式，能有幾千個了，有的時候，單個之內，就能再細分變化出幾十個，總量已經很難估算。
這個範例可能是電腦中最簡單的硬體應用範例，它直接操控 CPU 旁邊的一個計時器(timer)，以最快速的低階指令取出計時器的內容，拿來應用，我以它舉例，來觀測 Forth 系統中指令的執行速度，應用起來，可做出多種變化，可用來體會系統的執行速度。
這套應用技術的原始細節，是來自於我詳細研究 Forth 程式語言發明人 Charles H. Moore 在他向全世界公開之系統軟體 Color Forth 源程式後的延伸性應用。他設計 Color Forth 供自己用，用來設計 Forth CPU ，系統顯示的字體特別大，圖更是，顏色鮮明，彩色疊影照樣清晰，輸入一個指令就立刻被執行掉，系統緊密，能獨立運轉，輸入時可以只使用三個鍵，在 Windows‘ OS 中也有他親手寫成的版本。他手頭上有設計 CPU 時所需要的硬體模組圖檔，他老了，眼力不好，就用這樣的工具自己在家裡畫出生產 CPU 時所需要的罩刻圖樣，幾十層的疊影也照樣能夠處理。他有設計 CPU 硬體時最在乎的 Timing 數據庫，那一個零件模組傳訊號到另一個模組需要多少時延？都可以用 Color Forth 模擬出來，也就是他也自己測試時脈，讓設計成果能夠直接交付給像台積電那樣的公司，直接生產小奈米(nm)高頻率(>GHz)的 GA144 百蕊以上的 CPU ，產品很多年前就上市了，但我沒用過 GA144 。Color Forth 的介紹網頁如下，提供有興趣者參考：
https://colorforth.github.io/cf.htm
上訴提到用 Forth 模擬 Timing 測試的事情，這就需要動用到 CPU 旁邊的 Timer 來幫忙實現了，因此， Charles Moore 設計了它的驅動程式，組合語言寫的，非常簡單、快速，我就照抄來用，使用久了，也才發現，所有現行 64 位元系統，換成了新的 CPU 後，系統中的 Timer 並沒有跟著換成也是 64 位元的 Timer ，仍是 32 位元的 Timer。您若想改寫出 64 位元的驅動程式時，要學 Charles Moore 那樣，把取得的 32 位元值，轉供 64 位元的組合語言指令使用才行。
我所舉的這個範例很有趣，想測試任何程式的執行耗時多久時，您先執行一個 mark0 指定，然後輸入被測指令，只要不會令系統當機，給幾個都無所謂，最後面只需緊接著跟隨一個 mark1 指令，就行了。想看結果？則再使用 HowLong 指令，讓測試的數據印出來，就獲得了。
如果您想看我設計之取得時脈數的指令是怎麼設計的，執行
See GetTicks
便可看到全由低階組合語言設計的內容。
如果您對測試的動態能力沒有感覺，我再給您一列精簡程式，可以讓您直接看到 Timer 快速動態變化的內容，以 GHz 快速變化的數字，就顯示在所開小視窗的最上面一列內，想終止顯示，壓任何鍵就能恢復正常。
: test getxy begin 20 0 gotoxy get ticks 20 d.r key? Until gotoxy ;
現行電腦已經是以 GHz 的時脈在運轉，早在 APPLE II 的時代，我們就有過以 MHz 時脈來測量原子爐停爐棒跌降時間的技術，測量要求只須以毫秒(10^3分之一秒)等級顯示 1.8 秒內的跌降曲線。這個範例執行出來的結果，時間的精確度已達 10^9 分之一秒，能用來測量指令執行所耗用的時脈量。 GPS 或 BDS 衛星定位系統則以原子鐘等級的時脈供電波傳送耗用時間來算出定位，精確度等級到達了 10^15 分之一秒。這幾天科技期刊刊出了核子鐘(Nuclear clock)的技術，精確度能比原子鐘高出百萬倍，也就是到達了 10^21 分之一秒的程度了。無論時基精確度到達什麼程度，能精確測量，必須以電腦取得時基計量的數據來顯示或換算，這個範例程式在展示這種技術。這是核子鐘的貼文圖片：
:
繼前一個相當複雜的數學計算應用範例之後，今天介紹的範例，就顯得比較容易了解。
Forth 的能力，是能簡、能繁。繁的部分，改善出 BASIC 語法格式的程式設計能力，都能辦到。簡的部分，是 Forth 程式語言的特性，長久用、時常用、用久了就能體會，不用就難體會。
這些百例，非一天的收集而成，我在多年收集之後，每次回顧內容，都會產生新的感受，手癢到很想再做點變化，重新再試一試。也覺得百例仍能擴充出千例，只是願不願繼續做下去的問題，因為我粗略估算過，幾十年來，我寫成過的程式，能有幾千個了，有的時候，單個之內，就能再細分變化出幾十個，總量已經很難估算。
這個範例可能是電腦中最簡單的硬體應用範例，它直接操控 CPU 旁邊的一個計時器(timer)，以最快速的低階指令取出計時器的內容，拿來應用，我以它舉例，來觀測 Forth 系統中指令的執行速度，應用起來，可做出多種變化，可用來體會系統的執行速度。
這套應用技術的原始細節，是來自於我詳細研究 Forth 程式語言發明人 Charles H. Moore 在他向全世界公開之系統軟體 Color Forth 源程式後的延伸性應用。他設計 Color Forth 供自己用，用來設計 Forth CPU ，系統顯示的字體特別大，圖更是，顏色鮮明，彩色疊影照樣清晰，輸入一個指令就立刻被執行掉，系統緊密，能獨立運轉，輸入時可以只使用三個鍵，在 Windows‘ OS 中也有他親手寫成的版本。他手頭上有設計 CPU 時所需要的硬體模組圖檔，他老了，眼力不好，就用這樣的工具自己在家裡畫出生產 CPU 時所需要的罩刻圖樣，幾十層的疊影也照樣能夠處理。他有設計 CPU 硬體時最在乎的 Timing 數據庫，那一個零件模組傳訊號到另一個模組需要多少時延？都可以用 Color Forth 模擬出來，也就是他也自己測試時脈，讓設計成果能夠直接交付給像台積電那樣的公司，直接生產小奈米(nm)高頻率(>GHz)的 GA144 百蕊以上的 CPU ，產品很多年前就上市了，但我沒用過 GA144 。Color Forth 的介紹網頁如下，提供有興趣者參考：
https://colorforth.github.io/cf.htm
上訴提到用 Forth 模擬 Timing 測試的事情，這就需要動用到 CPU 旁邊的 Timer 來幫忙實現了，因此， Charles Moore 設計了它的驅動程式，組合語言寫的，非常簡單、快速，我就照抄來用，使用久了，也才發現，所有現行 64 位元系統，換成了新的 CPU 後，系統中的 Timer 並沒有跟著換成也是 64 位元的 Timer ，仍是 32 位元的 Timer。您若想改寫出 64 位元的驅動程式時，要學 Charles Moore 那樣，把取得的 32 位元值，轉供 64 位元的組合語言指令使用才行。
我所舉的這個範例很有趣，想測試任何程式的執行耗時多久時，您先執行一個 mark0 指定，然後輸入被測指令，只要不會令系統當機，給幾個都無所謂，最後面只需緊接著跟隨一個 mark1 指令，就行了。想看結果？則再使用 HowLong 指令，讓測試的數據印出來，就獲得了。
如果您想看我設計之取得時脈數的指令是怎麼設計的，執行
See GetTicks
便可看到全由低階組合語言設計的內容。
如果您對測試的動態能力沒有感覺，我再給您一列精簡程式，可以讓您直接看到 Timer 快速動態變化的內容，以 GHz 快速變化的數字，就顯示在所開小視窗的最上面一列內，想終止顯示，壓任何鍵就能恢復正常。
: test getxy begin 20 0 gotoxy get ticks 20 d.r key? Until gotoxy ;
現行電腦已經是以 GHz 的時脈在運轉，早在 APPLE II 的時代，我們就有過以 MHz 時脈來測量原子爐停爐棒跌降時間的技術，測量要求只須以毫秒(10^3分之一秒)等級顯示 1.8 秒內的跌降曲線。這個範例執行出來的結果，時間的精確度已達 10^9 分之一秒，能用來測量指令執行所耗用的時脈量。 GPS 或 BDS 衛星定位系統則以原子鐘等級的時脈供電波傳送耗用時間來算出定位，精確度等級到達了 10^15 分之一秒。這幾天科技期刊刊出了核子鐘(Nuclear clock)的技術，精確度能比原子鐘高出百萬倍，也就是到達了 10^21 分之一秒的程度了。無論時基精確度到達什麼程度，能精確測量，必須以電腦取得時基計量的數據來顯示或換算，這個範例程式在展示這種技術。這是核子鐘的貼文圖片：
:
\ (47)Howlong.f \ 應用 GerTicks 實測執行指令時，所耗用之CPU時脈(clock pulses)數量 \ get the ticks as soon as possible without any tedious code VARIABLE CLK0H VARIABLE CLK0L VARIABLE CLK1H VARIABLE CLK1L : Mark0 ( -- ) GetTicks CLK0H ! CLK0L ! ; : Mark1 ( -- ) GetTicks CLK1H ! CLK1L ! ; \ ignore the most significant bit by doing ABS to the high byte \ this bit is impossible changed in a short period of time while testing \ it can proved that ticks number is always to be a positive integer : HowLong ( -- ) CLK1L @ CLK1H @ ABS CLK0L @ CLK0H @ ABS DNEGATE D+ D. ; \ USAGE: \ (1) Mark0 NOOP Mark1 or \ (2) Mark0 3 5 + Mark1 \ (3) HowLong
第(48)個範例是一個亂數產生程式，精簡而扼要。
亂數函數是數學計算系統中的一個必要函數，我在設計系統時一定會建立，它有許多用途，所以是設計系統者不能不懂的東西。
＊亂數可以用來模擬系統狀況。
＊可以用來做事，例如供作抽籤程式之用。
＊可以用來快速提供測試系統時所需要的大量數據。
＊幾乎所有的電玩程式都要用到它。
＊很多求解數學方程式的問題，常用它提供的數據做為起始值。
＊人工智慧發展中，龐大矩陣的內容，就用亂數來產生極小的數值，做為代表尚未學習前的數學模式之狀態...等等等等應用。
上列用法，我幾乎都做過了。
亂數產生程式的設計方法，我用過許多種、也能講出十幾種，而且是 40 幾年前就已深入接觸過。
會設計亂數產生程式者，還得懂得評鑑亂數的品質，做過分析，才能確認方法夠好、可用。這個範例的方法，也經歷過這些發展過程，最後才落定於此，並經常被使用到。
此範例的方法叫做位元炒作法(bits cooking method)是炒作位元許多次後的輸出結果。炒作的方法很簡單，就是拿前一個存放於種子變數中的數字來炒作，炒作的原理就是先左移或右移質數量個次數，每次都將移位前與後的兩個數值進行一是一非則得真的 XOR 邏輯運算作為新結果，如此炒來炒去幾次，讓最高與最低效的位元都被變化到你很難想像它應該會是多少了，就炒成了，也得到了新的亂數，這就是這一套亂數產生函數的設計方法。
產生整數亂數是一切其他亂數的根本，浮點亂數，複數亂數也都可以憑藉整數亂數而產生。亂數在應用時，通常都另需指定使用範圍，處理出這樣的亂數，方法也有許多種，用時，不得不小心。
實際的應用上，另需注意一個原則性的問題，也就是無論您怎麼用亂數產生函數，由電腦產生的亂數，恆為一種虛亂之數(pseudo random number)，就只因為它是被電腦計算出來的，不是真正的自然地亂。後來，我增加了一個使用前先行亂化一下種子變數的技巧來改善狀況，亂化的方法，就是取用前一範例中取得之高速亂變時基的低效數值，32位元時直接放，64位元時放兩次。
網路上，英國學術界有人建立網頁提供亂數，介紹時號稱亂數的起源來自於大自然，他們使用英國北海多變氣象的數據亂組合而成，有興趣者可以自行搜索。我為了增廣見識，拜訪使用過該網頁，並介紹給本地華人組織，於春節聯歡會中臨時取得摸彩號碼之需。
2019 年， Google 公司發佈過他們新搞出了 53 個量子態的電腦，首個應用，不是計算，就是用來產生亂數，是量子態熱擾動後幾乎不耗時間的量子亂數產生函數應用。亂數被送進模擬系統使用後，據報導說是效果很好，瞬間就能得模擬系統的答案，詳情我並不熟知，但我知道用這種方法產生亂數的正面意義。
亂數產生函數都是計算產生的，所以很耗時間，因為一旦想要使用，用量必定都很龐大，產生得不快就不好。廣義地講，炒作位元的方法，還是用到邏輯運算，也是計算，只是相對於加減乘除的計算而言，比較不耗時間而已。所以，我在 64 位元環境中設計亂數產生程式時，都已改採位元炒作法產生亂數，我已固定了設計方法。
:
亂數函數是數學計算系統中的一個必要函數，我在設計系統時一定會建立，它有許多用途，所以是設計系統者不能不懂的東西。
＊亂數可以用來模擬系統狀況。
＊可以用來做事，例如供作抽籤程式之用。
＊可以用來快速提供測試系統時所需要的大量數據。
＊幾乎所有的電玩程式都要用到它。
＊很多求解數學方程式的問題，常用它提供的數據做為起始值。
＊人工智慧發展中，龐大矩陣的內容，就用亂數來產生極小的數值，做為代表尚未學習前的數學模式之狀態...等等等等應用。
上列用法，我幾乎都做過了。
亂數產生程式的設計方法，我用過許多種、也能講出十幾種，而且是 40 幾年前就已深入接觸過。
會設計亂數產生程式者，還得懂得評鑑亂數的品質，做過分析，才能確認方法夠好、可用。這個範例的方法，也經歷過這些發展過程，最後才落定於此，並經常被使用到。
此範例的方法叫做位元炒作法(bits cooking method)是炒作位元許多次後的輸出結果。炒作的方法很簡單，就是拿前一個存放於種子變數中的數字來炒作，炒作的原理就是先左移或右移質數量個次數，每次都將移位前與後的兩個數值進行一是一非則得真的 XOR 邏輯運算作為新結果，如此炒來炒去幾次，讓最高與最低效的位元都被變化到你很難想像它應該會是多少了，就炒成了，也得到了新的亂數，這就是這一套亂數產生函數的設計方法。
產生整數亂數是一切其他亂數的根本，浮點亂數，複數亂數也都可以憑藉整數亂數而產生。亂數在應用時，通常都另需指定使用範圍，處理出這樣的亂數，方法也有許多種，用時，不得不小心。
實際的應用上，另需注意一個原則性的問題，也就是無論您怎麼用亂數產生函數，由電腦產生的亂數，恆為一種虛亂之數(pseudo random number)，就只因為它是被電腦計算出來的，不是真正的自然地亂。後來，我增加了一個使用前先行亂化一下種子變數的技巧來改善狀況，亂化的方法，就是取用前一範例中取得之高速亂變時基的低效數值，32位元時直接放，64位元時放兩次。
網路上，英國學術界有人建立網頁提供亂數，介紹時號稱亂數的起源來自於大自然，他們使用英國北海多變氣象的數據亂組合而成，有興趣者可以自行搜索。我為了增廣見識，拜訪使用過該網頁，並介紹給本地華人組織，於春節聯歡會中臨時取得摸彩號碼之需。
2019 年， Google 公司發佈過他們新搞出了 53 個量子態的電腦，首個應用，不是計算，就是用來產生亂數，是量子態熱擾動後幾乎不耗時間的量子亂數產生函數應用。亂數被送進模擬系統使用後，據報導說是效果很好，瞬間就能得模擬系統的答案，詳情我並不熟知，但我知道用這種方法產生亂數的正面意義。
亂數產生函數都是計算產生的，所以很耗時間，因為一旦想要使用，用量必定都很龐大，產生得不快就不好。廣義地講，炒作位元的方法，還是用到邏輯運算，也是計算，只是相對於加減乘除的計算而言，比較不耗時間而已。所以，我在 64 位元環境中設計亂數產生程式時，都已改採位元炒作法產生亂數，我已固定了設計方法。
:
\ (48)RANDOM.F VARIABLE SEED 987654321 SEED ! : RENEWRANDOM ( -- ) GetTicks DROP SEED ! ; : COOKING ( n -- n ) DUP 13 LSHIFT XOR DUP 5 RSHIFT XOR DUP 7 LSHIFT XOR ; : RND ( -- ) SEED @ COOKING DUP SEED ! ABS 1- ; : RANDOM ( n -- n ) RND SWAP MOD ;
第(49)個範例就是亂數函數的實際應用。
應用的結果，不必去論其好壞，理論上，它就是確實能這樣用。
普天之下，求解圓周率的方法無數，很多方法也都被拿來評鑑電腦計算系統的好壞。求 pi 的技術，已經是一門古老的科技功夫，網上處處可見。我國古時候的名人祖沖之（西元 429 年 ～ 500 年），字文遠，范陽郡逎縣（今河北省保定市淶水縣）人，劉宋時代的數學家，就已經開始研究這個問題了。文獻記載過他提出之『密率』用法，為 355 除以 113 可得相當精確的 pi 值，直到今天，科技已經這麼先進，密率仍然很有用。我們 Forth 界的名著：『Starting Forth』，書中照樣引用密率。順便提示：相對的疏率，則是 22 除以 7 。
我在幾十年前學過 Monte Carlo 法的數值分析技術，於是設計了這個範例，可以用來解釋亂數產生程式的一種用法，就是使用亂數產生函數來求 pi 。
怎麼求？我學核子科學之父 Rutherford ，建立一種科學模型，比如以獵槍射擊乾草堆，以定量、定性的高能粒子轟擊固定的原子結構，例如完美的晶體，然後研究散射結果，就能解釋被轟擊物體的內容。
我的 Monte Carlo 法模型，則是以亂數產生一大群可以均勻分佈在一個正方形平面區域內所有座標點的數學模型，然後根據四分之一正方形內切圓的幾何條件，把圓內的座標點累積計數出來，圓外的座標點則放棄不要，就能得到結果。由於函數只取第一象限的關係，只得四分之一個圓的結果，而全圓又與四分之一圓上下左右對稱，因此，最後結果乘上四倍，就能得到 Monte Carlo 法算出來的 pi 。
話語講得通後，要寫成標準的 Forth 程式很不容易，根據幾何原理設計的判斷程式為 in-circle? 指令，根據亂數產生程式獲得座標點的程式為 point 指令。剩下其他的程式，就都是配合性的設計。
我自己模擬過使用結果，就列印在範例的最後面，數值僅能參考，因為您執行此範例程式後的所得結果，一定跟我的不一樣，這就是亂數的本質，對比之後，所得數據不是太離譜，就應該可以接受。
理論上，如果您能用量子電腦快速產生無限個亂數座標點後，您的所得結果必定就是圓周率 pi 。我才用了三個數字：一千、一萬跟十萬個座標點，當然就會只有那樣的不精確程度了。可是所得結果都不離譜，用十萬點時，已經具有小數點後面兩位數的精確度了： 3.14.... ，這就表示答案完全合理。
範例的目的，在介紹應用技巧，與實測結果，沒有人真願用這種方式來求得 pi 的。可是，有一天，若您也想用亂槍打鳥的方式模擬自然現象，研究某種問題時，學會這個範例，您就能知道如何運用 Monte Carlo 法分析問題了，而且，用的點數越多，結果就越精確。
您沒有風洞時，您可以有 Monte Carlo 法產生出來的點，照樣能進行氣動力學研究。
您沒有超臨界核子反應器時，您可以有 Monte Carlo 法產生之無數個起始中子狀態點，照樣能進行核子試爆研究。
您沒有經費對整個社會進行普查、取得數據時，您可以有 Monte Carlo 法產生出來的社會性狀分佈數據，照樣能進行國家政策想賴以決定之大事研究。
:
應用的結果，不必去論其好壞，理論上，它就是確實能這樣用。
普天之下，求解圓周率的方法無數，很多方法也都被拿來評鑑電腦計算系統的好壞。求 pi 的技術，已經是一門古老的科技功夫，網上處處可見。我國古時候的名人祖沖之（西元 429 年 ～ 500 年），字文遠，范陽郡逎縣（今河北省保定市淶水縣）人，劉宋時代的數學家，就已經開始研究這個問題了。文獻記載過他提出之『密率』用法，為 355 除以 113 可得相當精確的 pi 值，直到今天，科技已經這麼先進，密率仍然很有用。我們 Forth 界的名著：『Starting Forth』，書中照樣引用密率。順便提示：相對的疏率，則是 22 除以 7 。
我在幾十年前學過 Monte Carlo 法的數值分析技術，於是設計了這個範例，可以用來解釋亂數產生程式的一種用法，就是使用亂數產生函數來求 pi 。
怎麼求？我學核子科學之父 Rutherford ，建立一種科學模型，比如以獵槍射擊乾草堆，以定量、定性的高能粒子轟擊固定的原子結構，例如完美的晶體，然後研究散射結果，就能解釋被轟擊物體的內容。
我的 Monte Carlo 法模型，則是以亂數產生一大群可以均勻分佈在一個正方形平面區域內所有座標點的數學模型，然後根據四分之一正方形內切圓的幾何條件，把圓內的座標點累積計數出來，圓外的座標點則放棄不要，就能得到結果。由於函數只取第一象限的關係，只得四分之一個圓的結果，而全圓又與四分之一圓上下左右對稱，因此，最後結果乘上四倍，就能得到 Monte Carlo 法算出來的 pi 。
話語講得通後，要寫成標準的 Forth 程式很不容易，根據幾何原理設計的判斷程式為 in-circle? 指令，根據亂數產生程式獲得座標點的程式為 point 指令。剩下其他的程式，就都是配合性的設計。
我自己模擬過使用結果，就列印在範例的最後面，數值僅能參考，因為您執行此範例程式後的所得結果，一定跟我的不一樣，這就是亂數的本質，對比之後，所得數據不是太離譜，就應該可以接受。
理論上，如果您能用量子電腦快速產生無限個亂數座標點後，您的所得結果必定就是圓周率 pi 。我才用了三個數字：一千、一萬跟十萬個座標點，當然就會只有那樣的不精確程度了。可是所得結果都不離譜，用十萬點時，已經具有小數點後面兩位數的精確度了： 3.14.... ，這就表示答案完全合理。
範例的目的，在介紹應用技巧，與實測結果，沒有人真願用這種方式來求得 pi 的。可是，有一天，若您也想用亂槍打鳥的方式模擬自然現象，研究某種問題時，學會這個範例，您就能知道如何運用 Monte Carlo 法分析問題了，而且，用的點數越多，結果就越精確。
您沒有風洞時，您可以有 Monte Carlo 法產生出來的點，照樣能進行氣動力學研究。
您沒有超臨界核子反應器時，您可以有 Monte Carlo 法產生之無數個起始中子狀態點，照樣能進行核子試爆研究。
您沒有經費對整個社會進行普查、取得數據時，您可以有 Monte Carlo 法產生出來的社會性狀分佈數據，照樣能進行國家政策想賴以決定之大事研究。
:
\ Monte Carlo on the Pi.f
\ a place to keep a count of points that lie within the circle
\ (all points lie within the square, so no need to count them)
variable in-circle
\ square a number ( n -- n^2 )
: square dup * ;
\ generate a 2D point ( -- x y )
: point 1000 random 1000 random ;
\ determine if point lies within the circle ( x y -- flg )
: in-circle? square swap square + 1000000 < ;
\ run a simulation of n points ( n -- )
\ reset within circle counter
\ generate a point, test for within circle
\ count points within circle
\ do next point
: sim
0 in-circle !
0 do point in-circle?
if 1 in-circle +! then
loop ;
: pi sim 4 in-circle @ * cr . ;
((
run a simulation, calculate pi, and print result
(don't divide by area of square because we want all digits)
Here are some test runs:
1000 pi 3232
10000 pi 31472
100000 pi 314096
))
第(50)個範例是完全下載得來的佳例，用來把記憶體的內容傾印出來。
收集這個範例的主因有幾個。一是傾印方式跟 F83 系統的完全一樣，很好用也很恰到好處。二是創作者完全不以特殊指令來綁規格，從 16 位元到 32 到 64 位元的任何傳統 Forth 系統，都可以直接載入後便能使用。三是只有 GNU 版權，宣告後就能公開嵌入自己的創作，公諸於世。四是夠精簡，使用宣告而成的獨立緩衝區，放置處理中的暫態數據，使用時完全不受系統中其他指令或程式的影響。
為便於測試，我只將 dump 指令改成 ddump 後控存，以防載入時出現『指令用名重複』的警告訊息。
在發展 Forth 程式的過程中，免不了都會需要除錯工具。實際上， dump 指令早就是除錯工具之一，幾乎所有的 Forth 系統都直接提供。 ddump 它的用法，就是給一個起始位址與打算印出多少個位元組，就能執行 ddump 指令了。我給一個例子展示如何使用：
‘ ddump 40 - 100 ddump
就能看到這個 ddump 指令被編譯進系統後，存在於記憶體中整個的指令內容。程式介紹，上訴說明就夠了，改談一點這位貢獻公益程式的作者。
Krishna Myneni是一位印度裔的加拿大高能物理博士，是純以 c 設計出 kForth 系統的獨立作者。他在捐獻誤差函數 (erf) 程式給 Forth Scientifical Library 組織時，計畫主持人 Montgomery 請我審核程式，三位共同作者中唯一一位願意提供參考文獻給我的好心作者就是他。他也在讀過我的網文後，極力推薦，確認我的審核資格，我們因此成為網上的好朋友，可惜我不跑任何用 c 產生的 Forth 系統，沒能幫他。他也因我堅持傳統 Forth 的可貴，而開始修改 c 式 Forth 辦不到的設計，一改就熬了 20 幾年，還是沒有全面改成，也耽誤了他從 32 轉進 64 位元環境的開發進度，一直拖到去年，他才也搞出 64 位元 kForth 的胚胎版本。
但是，他很願意接受我的批評，聽我的建議後，設法把所有可取得的公益 Forth 程式，全面改寫成傳統 Forth 可以執行得起來的新版程式，外放給全世界免費使用。因此卻得罪了歐洲所有 Forth 系統的作者及該生態群體中的成員，常見他在國際論壇上被那一群惡徒攻擊，很多人都嫌他多事，把特殊用法又劣化回古老不良的設計格式，甚至於連個指令名稱妥不妥當之事都要攻擊他，德、英兩國的人最惡毒，可能是他的貢獻影響了英國的 MVP Forth 及德國的 iForth 兩個系統的銷售，及免費系統 gForth 跟用者的數量。
大家看了上述敘事，說到這裡，應該可以體會西方國家凡事永遠都想當盟主而打擊東方人的惡毒程度了，只因人種不同，西方國家就都聯合起來阻絕東方人，在 Forth 界也是如此，德國人的 gForth 就是用來腐化人心的學校集體創作，很多那個系統才能執行的公開程式，在傳統規格的 Forth 系統上都跑不起來，幾乎都是被故意弄成如此的。只能在 iForth 系統中執行的程式更是， MVP Forth 偷了美國 LMI Forth 的技術後，另行整合出優化程式碼的編譯工具，讓產生的成品什麼環境都不能用。
不知道這些底細的人，才會跟用歐洲的 Forth 系統，但擺明了，跟用者絕對摸不通整套 Forth 的精神，只能跟著用，也只能培養出二級程度的使用能力。台灣就有很多這種後來者，跟他們談 c 的部分還可以，無法深談 Forth 的核心技術。舉例來說，我用這個 ddump 看系統內容時，因歐式系統都強調防止使用者接觸 Forth 系統內的細部結構，想談就談不成。真談起來，還會發現他們的結構有很多但書，會讓執行情況經常轉向，所以談了也是白談。連單一個 COMPILE 指令，都得印出整整幾十頁的獨立檔案源程式後，才可能用來彼此溝通。這樣搞 Forth ，一輩子也搞不通，搞沒幾年，就必定會像搞流行程式語言那樣，退出群體不搞了，我看了四十幾年，屢見不鮮，很了解此前台灣的情況。四十幾年後，我還能在這裡跟大家談 Forth ，就是因為我沒中過 c 式 Forth 的毒，一直永遠講得出 Forth 的核心精神。 c 我不太熟悉，一直沒有什麼進展，我只接觸 Forth 也應必備的組合語言。
:
收集這個範例的主因有幾個。一是傾印方式跟 F83 系統的完全一樣，很好用也很恰到好處。二是創作者完全不以特殊指令來綁規格，從 16 位元到 32 到 64 位元的任何傳統 Forth 系統，都可以直接載入後便能使用。三是只有 GNU 版權，宣告後就能公開嵌入自己的創作，公諸於世。四是夠精簡，使用宣告而成的獨立緩衝區，放置處理中的暫態數據，使用時完全不受系統中其他指令或程式的影響。
為便於測試，我只將 dump 指令改成 ddump 後控存，以防載入時出現『指令用名重複』的警告訊息。
在發展 Forth 程式的過程中，免不了都會需要除錯工具。實際上， dump 指令早就是除錯工具之一，幾乎所有的 Forth 系統都直接提供。 ddump 它的用法，就是給一個起始位址與打算印出多少個位元組，就能執行 ddump 指令了。我給一個例子展示如何使用：
‘ ddump 40 - 100 ddump
就能看到這個 ddump 指令被編譯進系統後，存在於記憶體中整個的指令內容。程式介紹，上訴說明就夠了，改談一點這位貢獻公益程式的作者。
Krishna Myneni是一位印度裔的加拿大高能物理博士，是純以 c 設計出 kForth 系統的獨立作者。他在捐獻誤差函數 (erf) 程式給 Forth Scientifical Library 組織時，計畫主持人 Montgomery 請我審核程式，三位共同作者中唯一一位願意提供參考文獻給我的好心作者就是他。他也在讀過我的網文後，極力推薦，確認我的審核資格，我們因此成為網上的好朋友，可惜我不跑任何用 c 產生的 Forth 系統，沒能幫他。他也因我堅持傳統 Forth 的可貴，而開始修改 c 式 Forth 辦不到的設計，一改就熬了 20 幾年，還是沒有全面改成，也耽誤了他從 32 轉進 64 位元環境的開發進度，一直拖到去年，他才也搞出 64 位元 kForth 的胚胎版本。
但是，他很願意接受我的批評，聽我的建議後，設法把所有可取得的公益 Forth 程式，全面改寫成傳統 Forth 可以執行得起來的新版程式，外放給全世界免費使用。因此卻得罪了歐洲所有 Forth 系統的作者及該生態群體中的成員，常見他在國際論壇上被那一群惡徒攻擊，很多人都嫌他多事，把特殊用法又劣化回古老不良的設計格式，甚至於連個指令名稱妥不妥當之事都要攻擊他，德、英兩國的人最惡毒，可能是他的貢獻影響了英國的 MVP Forth 及德國的 iForth 兩個系統的銷售，及免費系統 gForth 跟用者的數量。
大家看了上述敘事，說到這裡，應該可以體會西方國家凡事永遠都想當盟主而打擊東方人的惡毒程度了，只因人種不同，西方國家就都聯合起來阻絕東方人，在 Forth 界也是如此，德國人的 gForth 就是用來腐化人心的學校集體創作，很多那個系統才能執行的公開程式，在傳統規格的 Forth 系統上都跑不起來，幾乎都是被故意弄成如此的。只能在 iForth 系統中執行的程式更是， MVP Forth 偷了美國 LMI Forth 的技術後，另行整合出優化程式碼的編譯工具，讓產生的成品什麼環境都不能用。
不知道這些底細的人，才會跟用歐洲的 Forth 系統，但擺明了，跟用者絕對摸不通整套 Forth 的精神，只能跟著用，也只能培養出二級程度的使用能力。台灣就有很多這種後來者，跟他們談 c 的部分還可以，無法深談 Forth 的核心技術。舉例來說，我用這個 ddump 看系統內容時，因歐式系統都強調防止使用者接觸 Forth 系統內的細部結構，想談就談不成。真談起來，還會發現他們的結構有很多但書，會讓執行情況經常轉向，所以談了也是白談。連單一個 COMPILE 指令，都得印出整整幾十頁的獨立檔案源程式後，才可能用來彼此溝通。這樣搞 Forth ，一輩子也搞不通，搞沒幾年，就必定會像搞流行程式語言那樣，退出群體不搞了，我看了四十幾年，屢見不鮮，很了解此前台灣的情況。四十幾年後，我還能在這裡跟大家談 Forth ，就是因為我沒中過 c 式 Forth 的毒，一直永遠講得出 Forth 的核心精神。 c 我不太熟悉，一直沒有什麼進展，我只接觸 Forth 也應必備的組合語言。
:
\ dump.4th \ \ memory dump utility \ \ Copyright (c) 1999 Krishna Myneni \ Creative Consulting for Research and Education \ \ This software is provided under the terms of the GNU \ General Public License. \ \ Last Revised: 4-25-1999 \ create dump_display_buf 20 allot : hexchar ( n -- m | return ascii hex char value for n: 0 - 15 ) dup 9 > if 10 - 65 + else 48 + then ; : dump_display_char ( n -- n|46 ) dup 33 < over 126 > or if drop 46 then ; : ddump ( a n -- | display n bytes starting at a ) dup 0> if 0 do i 16 mod 0= if cr dup . 58 emit 2 spaces then dup c@ 16 /mod hexchar emit hexchar emit 2 spaces dup c@ dump_display_char dump_display_buf i 16 mod + c! i 16 mod 15 = if dump_display_buf 16 type then 1+ loop drop else 2drop then ;
訂閱：
張貼留言 (Atom)
沒有留言:
張貼留言