一百個例題 (51 ~ 55)
Ching-Tang Tseng
Hamilton, New Zealand
15 September 2024
FORTH:慶的世界。 本網頁為中華民國(Republic Of China)FORTH語言協會 永久性網頁
2024年9月15日 星期日
一百個例題 (51 ~ 55)
第(51)個範例是一個完整的浮點數輸出指令的設計。
程式的來源是大家公認一直執有 Forth 商售系統牛耳的 Forth Inc. 公司的產品，但經過改寫後才能使用於 Win32Forth 系統。此範例是已經改好的程式，直接載入後便能測試。
開始討論前，不妨先去體會一下它的功能，再來細訴。測試時，我們也找一個有趣的事情體會系統性質，我教大家如何測試？
載入程式成功後，別管好幾條告訴您指令用名重複的警告訊息，只要最後有 ok 出現便可。
試前，我們先將輸出位數的設定值改大，執行下列命令就能把 32 位元浮點數的輸出位數宣告為上限量 18 位數：
18 fdpl !
然後，您從輸入 1.1e0 開始，直接要系統執行 fs. ，以一般標準科學輸出的方式，把這個浮點數印出來。繼續做 1.2e0 fs. 的事情，直到 1.9e0 fs. 。做完後，仔細檢查所有的輸入後立刻輸出的結果，您會發現，有幾個數字尾巴都帶有誤差值，不是絕對的 0 。這種狀況，就是前曾提及的二進制數字表示法之誤差本質，我們把數字送進系統，連算都沒有算，就叫系統重新印出來，誤差卻已經產生了。這種狀況，在十進制的浮點系統中，完全沒有，也就是誤差恆為 0 。
我並不期望大家有能力把印出數字時，轉換方法的整個程式細節，完全搞清楚，除非有一天，您也在幹自己設計浮點系統之事，否則先別去浪費時間，研究冗長又無聊的細節。就算練成這個功夫，將來可能還會用不上，因為浮點系統從未成為標準，也沒必要成為標準。我追蹤過整套程式後，才收集它到這裡當個範例，但也從未用上過它，我自己的浮點設計是十進制的浮點系統，所以用不上。直接引用別人具有浮點運算能力的系統時，這些浮點輸出指令都是現成的，不勞我費心，因此，這套範例我就從未用過。
我之所以仍將程式收集成範例來留參的主因，是因為這個程式設計得很好，因為它把印出浮點數的整套指令，設計成都對應到印出一般整數的規格上去了，所有的 <#. #. #S. #>. SIGN. 浮點數逐個位數的轉換指令，完全對應到標準的整數之逐個位數的轉換指令 <# # #S #> SIGN。這種設計方法，是一種程式設計的典範，不忘本而延伸，可以令後繼者不必看設計細節，就能知道指令可以用來幹什麼？以及應該如何使用？
大家都知道一般 Forth 系統建立的初期，根本就不考慮浮點數的問題，我愛用的 Lina64 就是。原因很多，不必深究，也因此，浮點系統沒有統一標準，若有標準，會更麻煩，每換一套軟硬體，標準就必須跟著重訂，很不實際。早期的 Forth 系統，也因此不太適合使用於數值分析的發展或應用上， FORTH Inc. 公司所有推出的產品，就是很好的說明。這家公司公開宣稱過，他們不擅長於浮點運算系統的技術，主要的成就，大部份偏重於可以只用整數的單板電腦控制上，為太空計劃、天文資料自動收集發展的系統、或為阿拉伯利雅德機場發展的巨型航管系統，都沒有浮點運算的絕對需求，全能只用整數運算完成，這都是事實。
當初 APPLE II 的時代，有一個領先群雄具有浮點運算功能的 MicroMotion FORTH-79 系統，曾被我們在台灣作為推廣標準。這個系統的作者 Martin Tracy 就是個浮點系統的設計好手，後來，他日正當中時賣掉了系統，轉而效力於 Forth Inc. 公司，專門負責浮點系統的開發。可惜，那時的環境只有單板電腦的微處理機應用方式，完全沒有多媒體應用的影子，浮點運算功能的技術發展就式微了， Martin Tracy 也離開了 Forth Inc. 改去效勞微軟公司，主持多媒體套裝軟體程式的開發，也中止了他曾為 20 幾種 CPU 設計過 Forth 系統的繼續研發，非常可惜。從事於影像處裡的工作時絕對需要使用浮點運算， Martin Tracy 的專長就更能發揮，所以他效勞於微軟是正確的選擇。
我們回過頭來討論浮點系統，每個人想自行設計浮點系統時，最先要做的事情，並不是馬上開始設計四則運算，而是應該先行設計出浮點數字的輸出、輸入指令，否則，設計者如何能知道自己的後續設計能執行得出正確的結果？ 而且，輸出比輸入該先設計，因為還沒有輸入指令前，可以先只使用直接放置位元花樣在記憶體上的方式，達到等效的輸入目的。因此，浮點數的輸出指令，是設計浮點系統前的最重要設計，有了它，後續事情就好辦多了。請注意！這套範例的最大缺點，是用浮點運算中的四則運算來處理浮點數，所以，變成了您必須先有運算功能後，才能有這個範例提供的那種漂亮輸出。因此，此範例也不能用於系統設計之初。
浮點輸出的常用指令為 f. fs. Fe. ，相對於它們的核心指令是 (f.) (fs.) (fe.)，它們的意義及用法，是使用 Forth 時很基本的概念，我不擬浪費時間講得這麼細。
我在介紹百例超過一半後，大家應該可以體會到，這些範例，實際上有許多個，都是我在設計系統時，所必須牽扯到的技術。所以，若有人也想設計系統，就不能不搞通這些技術。想要善用新增添出來的系統功能，那就更應該仔細了解這些分項獨立出來的範例討論。我們還有一半的範例可以探討。
:
程式的來源是大家公認一直執有 Forth 商售系統牛耳的 Forth Inc. 公司的產品，但經過改寫後才能使用於 Win32Forth 系統。此範例是已經改好的程式，直接載入後便能測試。
開始討論前，不妨先去體會一下它的功能，再來細訴。測試時，我們也找一個有趣的事情體會系統性質，我教大家如何測試？
載入程式成功後，別管好幾條告訴您指令用名重複的警告訊息，只要最後有 ok 出現便可。
試前，我們先將輸出位數的設定值改大，執行下列命令就能把 32 位元浮點數的輸出位數宣告為上限量 18 位數：
18 fdpl !
然後，您從輸入 1.1e0 開始，直接要系統執行 fs. ，以一般標準科學輸出的方式，把這個浮點數印出來。繼續做 1.2e0 fs. 的事情，直到 1.9e0 fs. 。做完後，仔細檢查所有的輸入後立刻輸出的結果，您會發現，有幾個數字尾巴都帶有誤差值，不是絕對的 0 。這種狀況，就是前曾提及的二進制數字表示法之誤差本質，我們把數字送進系統，連算都沒有算，就叫系統重新印出來，誤差卻已經產生了。這種狀況，在十進制的浮點系統中，完全沒有，也就是誤差恆為 0 。
我並不期望大家有能力把印出數字時，轉換方法的整個程式細節，完全搞清楚，除非有一天，您也在幹自己設計浮點系統之事，否則先別去浪費時間，研究冗長又無聊的細節。就算練成這個功夫，將來可能還會用不上，因為浮點系統從未成為標準，也沒必要成為標準。我追蹤過整套程式後，才收集它到這裡當個範例，但也從未用上過它，我自己的浮點設計是十進制的浮點系統，所以用不上。直接引用別人具有浮點運算能力的系統時，這些浮點輸出指令都是現成的，不勞我費心，因此，這套範例我就從未用過。
我之所以仍將程式收集成範例來留參的主因，是因為這個程式設計得很好，因為它把印出浮點數的整套指令，設計成都對應到印出一般整數的規格上去了，所有的 <#. #. #S. #>. SIGN. 浮點數逐個位數的轉換指令，完全對應到標準的整數之逐個位數的轉換指令 <# # #S #> SIGN。這種設計方法，是一種程式設計的典範，不忘本而延伸，可以令後繼者不必看設計細節，就能知道指令可以用來幹什麼？以及應該如何使用？
大家都知道一般 Forth 系統建立的初期，根本就不考慮浮點數的問題，我愛用的 Lina64 就是。原因很多，不必深究，也因此，浮點系統沒有統一標準，若有標準，會更麻煩，每換一套軟硬體，標準就必須跟著重訂，很不實際。早期的 Forth 系統，也因此不太適合使用於數值分析的發展或應用上， FORTH Inc. 公司所有推出的產品，就是很好的說明。這家公司公開宣稱過，他們不擅長於浮點運算系統的技術，主要的成就，大部份偏重於可以只用整數的單板電腦控制上，為太空計劃、天文資料自動收集發展的系統、或為阿拉伯利雅德機場發展的巨型航管系統，都沒有浮點運算的絕對需求，全能只用整數運算完成，這都是事實。
當初 APPLE II 的時代，有一個領先群雄具有浮點運算功能的 MicroMotion FORTH-79 系統，曾被我們在台灣作為推廣標準。這個系統的作者 Martin Tracy 就是個浮點系統的設計好手，後來，他日正當中時賣掉了系統，轉而效力於 Forth Inc. 公司，專門負責浮點系統的開發。可惜，那時的環境只有單板電腦的微處理機應用方式，完全沒有多媒體應用的影子，浮點運算功能的技術發展就式微了， Martin Tracy 也離開了 Forth Inc. 改去效勞微軟公司，主持多媒體套裝軟體程式的開發，也中止了他曾為 20 幾種 CPU 設計過 Forth 系統的繼續研發，非常可惜。從事於影像處裡的工作時絕對需要使用浮點運算， Martin Tracy 的專長就更能發揮，所以他效勞於微軟是正確的選擇。
我們回過頭來討論浮點系統，每個人想自行設計浮點系統時，最先要做的事情，並不是馬上開始設計四則運算，而是應該先行設計出浮點數字的輸出、輸入指令，否則，設計者如何能知道自己的後續設計能執行得出正確的結果？ 而且，輸出比輸入該先設計，因為還沒有輸入指令前，可以先只使用直接放置位元花樣在記憶體上的方式，達到等效的輸入目的。因此，浮點數的輸出指令，是設計浮點系統前的最重要設計，有了它，後續事情就好辦多了。請注意！這套範例的最大缺點，是用浮點運算中的四則運算來處理浮點數，所以，變成了您必須先有運算功能後，才能有這個範例提供的那種漂亮輸出。因此，此範例也不能用於系統設計之初。
浮點輸出的常用指令為 f. fs. Fe. ，相對於它們的核心指令是 (f.) (fs.) (fe.)，它們的意義及用法，是使用 Forth 時很基本的概念，我不擬浪費時間講得這麼細。
我在介紹百例超過一半後，大家應該可以體會到，這些範例，實際上有許多個，都是我在設計系統時，所必須牽扯到的技術。所以，若有人也想設計系統，就不能不搞通這些技術。想要善用新增添出來的系統功能，那就更應該仔細了解這些分項獨立出來的範例討論。我們還有一半的範例可以探討。
:
\ Simple Floating-Point Output
\ Revision 2013-10-29
\ This simple floating-point output package has features found in more comprehensive implementations yet
\ is remarkably compact and portable. Based on code and algorithm from Forth Inc.
\ SFPOUT.F
\
\ Simple Floating Point Output
\
\ Main words:
\
\ (F.) (FS.) (FE.) F.R FS.R FE.R F. FS. FE.
\ FDP PLACES
\
\ This package should function correctly on any Forth
\ system with the following limitations:
\
\ - Don't attempt to output non-numbers such as NANs
\ or INFs as it will enter an infinite loop.
\ - Floating-point strings are limited to the size of
\ the system's pictured numeric output buffer.
\
\ History:
\ 131029 Fix (F.) to use FDP. Add F. FS. FE. PLACES
FORTH DEFINITIONS DECIMAL
\ Floating-point pictured numeric output operators
: <#. ( F: r1 -- r2 ) FROUND <# ;
: #. ( F: r1 -- r2 ) 10.E F/ FDUP FLOOR FSWAP FOVER F-
10.E F* FROUND F>D D>S [CHAR] 0 + HOLD ;
: #S. ( F: r1 -- r2 ) BEGIN #. FDUP F0= UNTIL ;
: #>. ( F: r -- ) ( c-addr u ) FDROP 0 0 #> ;
: SIGN. ( flag -- ) IF [CHAR] - HOLD THEN ;
\ Variable controlling trailing decimal point display.
\ Default (ON) is to always display decimal point.
VARIABLE FDP 1 FDP !
: 10^n ( r1 n -- r2 ) 0 ?DO 10.E F* LOOP ;
: #.n ( r1 n -- r2 ) 0 ?DO #. LOOP ;
VARIABLE rscale 1 rscale !
FVARIABLE rstep 10.E rstep F!
VARIABLE fdpl 4 fdpl !
\ Normalize to range 1.0 <= r < STEPSIZE
: fnorm ( r1 -- |r2| sign exp )
FDUP F0< 0 2>R FABS
FDUP F0= 0= IF
BEGIN FDUP rstep F@ F< 0=
WHILE rstep F@ F/ R> rscale @ + >R REPEAT
BEGIN FDUP 1.0E F<
WHILE rstep F@ F* R> rscale @ - >R REPEAT
THEN 2R> ;
\ Convert fixed-point
: fcvt ( r n -- )
>R FDUP F0< ( sign) R> 2>R
FABS FDP @ IF ( always output decimal point )
R> #.n [CHAR] . HOLD
ELSE ( conditionally output decimal point )
R@ #.n R> IF [CHAR] . HOLD THEN
THEN #S. R> SIGN. #>. ;
\ Convert real number r to string c-addr u in exponential
\ notation with n places right of the decimal point.
: (e.) ( r n scale step -- c-addr u )
rstep F! rscale ! 0 MAX >R fnorm
R> 2>R IF FNEGATE THEN 1.E R@ 10^n
FSWAP FOVER F* FROUND ( make integer)
FDUP FABS FROT F/ rstep F@ F< 0= IF ( overflow)
rstep F@ F/ R> R> rscale @ + >R >R THEN
<#. R> R> S>D TUCK DABS # #S 2DROP
0< IF [CHAR] - ELSE [CHAR] + THEN HOLD
[CHAR] E HOLD fcvt ;
\ Convert real number r to string c-addr u in scientific
\ notation with n places right of the decimal point.
: (FS.) ( r n -- c-addr u ) 1 10.E (e.) ;
\ Display real number r in scientific notation right-
\ justified in a field width u with n places right of
\ the decimal point.
: FS.R ( r n u -- ) >R (FS.) R> OVER - SPACES TYPE ;
\ Convert real number r to string c-addr u in engineering
\ notation with n places right of the decimal point.
: (FE.) ( r n -- c-addr u ) 3 1000.E (e.) ;
\ Display real number r in engineering notation right-
\ justified in a field width u with n places right of
\ the decimal point.
: FE.R ( r n u -- ) >R (FE.) R> OVER - SPACES TYPE ;
\ Convert real number r to string c-addr u in fixed-point
\ notation with n places right of the decimal point.
: (F.) ( r n -- c-addr u )
0 MAX DUP >R 10^n <#. ( round) R> fcvt ;
\ Display real number r in fixed-point notation right-
\ justified in a field width u with n places right of
\ the decimal point.
: F.R ( r n u -- ) >R (F.) R> OVER - SPACES TYPE ;
\ Set decimal places control for F. FS. FE.
: PLACES ( n -- ) fdpl ! ;
: F. ( r -- ) fdpl @ 0 F.R SPACE ;
: FS. ( r -- ) fdpl @ 0 FS.R SPACE ;
: FE. ( r -- ) fdpl @ 0 FE.R SPACE ;
[DEFINED] DXFORTH [IF] behead 10^n (e.) [THEN]
\ end
第(52)個範例是個很簡單的 FORTRAN 轉寫成 ABC Forth 格式的精簡展示。
Kincaid-Cheney 所著的數值分析課本相當普及，這個範例是該教科書中的第一個例題。
源程式我就收集在此例後面，以便對照著看出改寫方法，顯示其簡單性。
改寫 FORTRAN 程式時，有幾個要點，掌握住後，改寫起來，就沒有任何困難。改寫者甚至於可以完全不管程式技巧，前人能寫得多漂亮，改寫出來的成品，就可以也有多漂亮。這一個範例只是簡單到可以用來解釋要點，才會被收集在此處。
前面曾經提及，使用 ABC Forth 時，首要重點就是先宣告出資料結構，用到多少個，就宣告出多少個。後續的程式，都僅只是在對宣告出來的資料結構進行數學運算操作，如此而已。
改寫 FORTRAN 程式時，最重要的地方，就是完全不要理會所有的輸入與輸出指令，因為 FORTRAN 執行 I/O 的指令受到格式必須先經過宣告的限制，而 BASIC 程式語言的簡易性，就凸顯在 I/O 格式比 FORTRAN 簡單太多了。這種差異，也是我寫出 ABC Forth 的主因。
仔細看看 FORTRAN 的程式格式，它又何嘗不是先宣告資料結構才寫運算程式？只不過在用到極少數的幾個指標時， FORTRAN 設計成是現成可用，而 ABC Forth 則不打算節省，如此而已。搞清楚這種很基本之哲理上的差異，就能明白系統的設計方法了。
本程式沒有邏輯分支，所以程式簡單。而唯一的一個定量迴路，如果硬要用傳統 Forth 來寫，對應性差異太大，但若以 ABC Forth 改寫，就幾乎是完全等效了。這也是為什麼我在寫數值分析程式時，堅持只用 ABC Forth 而拒絕再用傳統 Forth 的主要原因之一。取其改寫時絕不傷神的優點，改寫就不易出錯。他們之間的對應關係為：
DO 2 n=1,1000
..........
2. CONTINUE
對應到：
FOR n = 1 to 1000
..........
NEXT n
至於程式本身，是用來強調一個 1.幾 很小的基底數字，在遭遇到很大方次的狀況下，系統還能不會爆掉的一種收斂式級數現象。本來，就算只用 2 為基底，求他的很高次方時，很快，系統就會爆掉。這個範例採用一個計算 1.幾 ，只比 1 多一點點的微小增量數字來算，算到 1000 次方了，系統也不會爆掉。這樣，作者就稱其為是足夠收斂的級數了。
資料結構宣告了一個具有 1000 個單元的陣列，範例中也有對應的宣告方式可以參考。
:
Kincaid-Cheney 所著的數值分析課本相當普及，這個範例是該教科書中的第一個例題。
源程式我就收集在此例後面，以便對照著看出改寫方法，顯示其簡單性。
改寫 FORTRAN 程式時，有幾個要點，掌握住後，改寫起來，就沒有任何困難。改寫者甚至於可以完全不管程式技巧，前人能寫得多漂亮，改寫出來的成品，就可以也有多漂亮。這一個範例只是簡單到可以用來解釋要點，才會被收集在此處。
前面曾經提及，使用 ABC Forth 時，首要重點就是先宣告出資料結構，用到多少個，就宣告出多少個。後續的程式，都僅只是在對宣告出來的資料結構進行數學運算操作，如此而已。
改寫 FORTRAN 程式時，最重要的地方，就是完全不要理會所有的輸入與輸出指令，因為 FORTRAN 執行 I/O 的指令受到格式必須先經過宣告的限制，而 BASIC 程式語言的簡易性，就凸顯在 I/O 格式比 FORTRAN 簡單太多了。這種差異，也是我寫出 ABC Forth 的主因。
仔細看看 FORTRAN 的程式格式，它又何嘗不是先宣告資料結構才寫運算程式？只不過在用到極少數的幾個指標時， FORTRAN 設計成是現成可用，而 ABC Forth 則不打算節省，如此而已。搞清楚這種很基本之哲理上的差異，就能明白系統的設計方法了。
本程式沒有邏輯分支，所以程式簡單。而唯一的一個定量迴路，如果硬要用傳統 Forth 來寫，對應性差異太大，但若以 ABC Forth 改寫，就幾乎是完全等效了。這也是為什麼我在寫數值分析程式時，堅持只用 ABC Forth 而拒絕再用傳統 Forth 的主要原因之一。取其改寫時絕不傷神的優點，改寫就不易出錯。他們之間的對應關係為：
DO 2 n=1,1000
..........
2. CONTINUE
對應到：
FOR n = 1 to 1000
..........
NEXT n
至於程式本身，是用來強調一個 1.幾 很小的基底數字，在遭遇到很大方次的狀況下，系統還能不會爆掉的一種收斂式級數現象。本來，就算只用 2 為基底，求他的很高次方時，很快，系統就會爆掉。這個範例採用一個計算 1.幾 ，只比 1 多一點點的微小增量數字來算，算到 1000 次方了，系統也不會爆掉。這樣，作者就稱其為是足夠收斂的級數了。
資料結構宣告了一個具有 1000 個單元的陣列，範例中也有對應的宣告方式可以參考。
:
\ elimit.f
integer i
2 reals x n
1000 array x
: x(n)
BASIC
10 for i = 1 to 1000
20 let { n = i>r ( i ) }
20 let { x ( i ) = ( 1 + 1 / ln ( n ) ) ^ n }
30 next i
40 end
;
: main
BASIC
10 run x(n)
20 print " x( 1) = " ; { x ( 1 ) }
30 print " x( 10) = " ; { x ( 10 ) }
40 print " x( 100) = " ; { x ( 100 ) }
50 print " x(1000) = " ; { x ( 1000 ) }
60 end
;
\s
c file: elimit.f
c
double precision x
dimension x(1000)
c
print *
print *,' Slowly converging sequence for irrational number e'
print *,' Section 1.2, Kincaid-Cheney'
print *
c
do 2 n=1,1000
x(n) = (1.0d0 + 1.0d0/dble(n))**n
2 continue
c
print *,' 1, x(1) =',x(1)
print *,' 10, x(10) =',x(10)
print *,' 30, x(30) =',x(30)
print *,' 50, x(50) =',x(50)
print *,' 1000, x(1000) =',x(1000)
print *,' exp(1.0) =',exp(1.0d0)
print *
print *,' error =',abs(x(1000) - exp(1.0d0))
c
stop
end
第(53)個範例，算出浮點數輸入時，系統以累積之 2 的負幾次方之分數量來換算出最接近輸入值的對應數值時，所需要用到的實際分數量，逐個印出，直到負的 99 次方為止。
程式的來源，還是 Kincaid-Cheney 所著的數值分析課本第二章的例題。此題除了有『定量迴路』，還有『邏輯分支』，所以我刻意採用來當改寫範例。
我改寫後，發現了這份例題有點問題，若直接改寫，照原來的執行步驟，邏輯判斷是多餘的安排。
只看 FORTRAN 程式很不容易看得出來，因為 (t .le. 1.0) 指令成立時，才會停止執行印出計算結果。但是前一列程式卻是 t = s + 1.0 ，t 怎麼可能小於 1 ? 除非前一列的 0.5*s 得到了大於 1.0 的負值，但永遠不可能是負值，這就很不合理，邏輯判斷後分支，等於是永遠沒被執行。
最後，我不理會這個問題了，只考慮這個範例的主旨是在印出所有 2 的負幾次方之值，乾脆自己重寫一個 main2 解決就算了，獲得的答案，除了負的 100 次方也印出來外，完全一樣。
我在抄改寫舊程式時，經常會碰到這類問題，有時覺得問題好像就是原作者故意安排的陷阱，學習者若不仔細實作，根本不可能知道問題發生在那裡？
我在發展這些技術時，為了能讓自己永久留參，牢記有問題的地方，忠實留下發展記錄，絕不盲從。錯誤最能考驗系統性能的正常與否？這是這個範例被留下來的主要特色。
:
程式的來源，還是 Kincaid-Cheney 所著的數值分析課本第二章的例題。此題除了有『定量迴路』，還有『邏輯分支』，所以我刻意採用來當改寫範例。
我改寫後，發現了這份例題有點問題，若直接改寫，照原來的執行步驟，邏輯判斷是多餘的安排。
只看 FORTRAN 程式很不容易看得出來，因為 (t .le. 1.0) 指令成立時，才會停止執行印出計算結果。但是前一列程式卻是 t = s + 1.0 ，t 怎麼可能小於 1 ? 除非前一列的 0.5*s 得到了大於 1.0 的負值，但永遠不可能是負值，這就很不合理，邏輯判斷後分支，等於是永遠沒被執行。
最後，我不理會這個問題了，只考慮這個範例的主旨是在印出所有 2 的負幾次方之值，乾脆自己重寫一個 main2 解決就算了，獲得的答案，除了負的 100 次方也印出來外，完全一樣。
我在抄改寫舊程式時，經常會碰到這類問題，有時覺得問題好像就是原作者故意安排的陷阱，學習者若不仔細實作，根本不可能知道問題發生在那裡？
我在發展這些技術時，為了能讓自己永久留參，牢記有問題的地方，忠實留下發展記錄，絕不盲從。錯誤最能考驗系統性能的正常與否？這是這個範例被留下來的主要特色。
:
\ 2**(-n)
((
c
c Second Edition
c Numerical Analysis:
c The Mathematics of Scientific Computing
c D.R. Kincaid and E.W. Cheney
c Brooks/Cole Publ., 1996
c ISBN 0-534-33892-5
c COPYRIGHT (c) 1996
c
c Section 2.1
c
c Computing an approximate value of machine precision
c
c
c file: epsi.f
c
print *
print *,' Approximate value of machine precision'
print *,' Section 2.1, Kincaid-Cheney'
print *
print *,' n computed 2**(-n)'
c
s = 1.0
c
do 2 k=1,100
s = 0.5*s
t = s + 1.0
if (t .le. 1.0) then
s = 2.0*s
t = 1.0/2.0**(k-1)
print 3,k-1,s,t
stop
endif
2 continue
c
3 format (i3,2x,2(e13.6,2x))
stop
end
))
2 integers k k-1
2 reals s t
20 sigdigits !
: main
BASIC
10 let { s = 1 }
20 for k = 1 to 100
30 let { s = 0.5 * s }
40 let { t = s + 1 }
50 if { t <= 1 } then 100
60 let { s = 2 * s }
70 let k-1 = k - 1
80 let { t = 1 / ( 2 ^ i>r ( k-1 ) ) }
90 print " " ; k-1 , { s , t }
100 next k
110 end
;
: main2
basic
10 for k = 1 to 100
20 let { t = 1 / ( 2 ^ i>r ( k ) ) }
30 print " " ; k , { t }
40 next k
50 end
;
第(54)個範例設計出排列與組合的萬用程式。
舉例而言:從 n 個不同數字取出 r 個來，可以有多少種花樣？ 假設 n = 10 , r = 3 。
如果不論取出之 3 個數字的出現順序，如(1,4,7)與(4,7,1)...等，都只能算是同一花樣時，就叫做組合(combination)，我們以 nCr 符號表示所有可能花樣的總數，計算公式為：nCr = n! / ( (n-r)! * r! ) 。
如果考慮出現的順序不同時，均需各別當作為一種不同的花樣，就叫做排列(permutation)，我們以 nPr 符號表示所有可能花樣的總數，計算公式為：nPr = n! / (n-r)!。 顯然，同一條件下，排列的總可能數比組合的總可能數要多。
以上為排列與組合在統計學上基本定義的簡單回顧，所有的數學計算系統都必須能夠用來算出排列與組合數字，它們相當於也是一種函數，只是像取最大值或取最小值的情況一樣，必須有兩個輸入參數，只輸出單一個結果數字。
單純只求排列與組合答案的程式不難設計，但因有上列公式可以直接演算，所以系統通常並不提供這種函數。演算出階層函數(factorial function)的功能，才是此類計算的根本。這是一個只對整數操作的基本函數，可是，就算 64 位元的系統，可用來表示的整數範圍，仍然相當有限，階層函數就不太被系統設計者重視，我的系統則提供。
另外還有一個因素令系統設計者不太願意提供階層函數，如果您仔細思考各種慣用的函數表示方式，您會發現，類似 sin , cos 之類的函數，是前算符式的函數，寫成 sin x 。其實，四則運算在系統中也是一種函數，但為中算符式函數，寫成 x + y 。而階層函數則為後算符式函數，寫成 x ! 。這樣的差異現象，會令系統設計者難以安排三種函數之正確出現位置能被輕易的設計出來。我辦得到，所以我的系統除了它人已有的四則與常用三角函數外，也有階層函數，而且就可以寫成類似 10 ! 的格式來設計程式。
這個範例就展示排列與組合兩個函數的設計總成，其隱含能力卻不是一般程式語言能辦得到的事情。因為，我的系統能處理無限位數的整數計算，才能設計得出很大數字的排列與組合函數來。有人想刁難系統時，這個範例足夠用來對付。超出慣常 64 位元整數範圍的問題時，這個範例也當然能夠解決問題。對付許多在理論數學中研究數字理論的問題，這個範例也很有用。網上早就有許多古怪的問題，專門用來刁難系統，例如硬要您求取非常高階之二項方程式展開式中某一項的係數為何？這種問題，使用我的系統，就很容易解決。
範例的主旨在展示計算無限位數整數運算的程式該怎麼寫？請注意！除了必須使用 b{ . . . . … }b 把演算式子包裹起來外，沒有太多的差異。四則運算的加減乘除仍可照用，只是必須牢記住，大數只能與大數運算，不是大數，格式結構當然不同，想用，就得先進行格式轉換函數的操作，才能使用。還有，想印出結果時，不能隨便印了，只能用 Big. 來印。注意到 BASIC 式的指令都沒有受到影響了嗎？這就是我設計系統的技巧，可以還有更多的數學體系被加進來，照樣可用 BASIC 格式的程式操作。
能算很大的數字，當然也得包括能算很小的數字，紐西蘭的樂透(lottory)是用 40 個數字選出 6 個數字的方式開獎，兩塊錢買一張票，中頭獎一百萬的機率是三百八十三萬八千三百八十分之一( 1 / 383,8380 )。花兩塊錢去賭這一百萬，每一塊錢的作用效率，也就是中獎機率，實在是太低了，我從不買。若算成排列，或然率就更低，是： 1 / 27,6363,3600，約27億分之一。規矩若改成排列的結果，就更不會有人願意賭了。請輸入下列操作，自己試一試：
40 6 nCr big.
40 6 nPr big.
1000 123 nCr big.
:
舉例而言:從 n 個不同數字取出 r 個來，可以有多少種花樣？ 假設 n = 10 , r = 3 。
如果不論取出之 3 個數字的出現順序，如(1,4,7)與(4,7,1)...等，都只能算是同一花樣時，就叫做組合(combination)，我們以 nCr 符號表示所有可能花樣的總數，計算公式為：nCr = n! / ( (n-r)! * r! ) 。
如果考慮出現的順序不同時，均需各別當作為一種不同的花樣，就叫做排列(permutation)，我們以 nPr 符號表示所有可能花樣的總數，計算公式為：nPr = n! / (n-r)!。 顯然，同一條件下，排列的總可能數比組合的總可能數要多。
以上為排列與組合在統計學上基本定義的簡單回顧，所有的數學計算系統都必須能夠用來算出排列與組合數字，它們相當於也是一種函數，只是像取最大值或取最小值的情況一樣，必須有兩個輸入參數，只輸出單一個結果數字。
單純只求排列與組合答案的程式不難設計，但因有上列公式可以直接演算，所以系統通常並不提供這種函數。演算出階層函數(factorial function)的功能，才是此類計算的根本。這是一個只對整數操作的基本函數，可是，就算 64 位元的系統，可用來表示的整數範圍，仍然相當有限，階層函數就不太被系統設計者重視，我的系統則提供。
另外還有一個因素令系統設計者不太願意提供階層函數，如果您仔細思考各種慣用的函數表示方式，您會發現，類似 sin , cos 之類的函數，是前算符式的函數，寫成 sin x 。其實，四則運算在系統中也是一種函數，但為中算符式函數，寫成 x + y 。而階層函數則為後算符式函數，寫成 x ! 。這樣的差異現象，會令系統設計者難以安排三種函數之正確出現位置能被輕易的設計出來。我辦得到，所以我的系統除了它人已有的四則與常用三角函數外，也有階層函數，而且就可以寫成類似 10 ! 的格式來設計程式。
這個範例就展示排列與組合兩個函數的設計總成，其隱含能力卻不是一般程式語言能辦得到的事情。因為，我的系統能處理無限位數的整數計算，才能設計得出很大數字的排列與組合函數來。有人想刁難系統時，這個範例足夠用來對付。超出慣常 64 位元整數範圍的問題時，這個範例也當然能夠解決問題。對付許多在理論數學中研究數字理論的問題，這個範例也很有用。網上早就有許多古怪的問題，專門用來刁難系統，例如硬要您求取非常高階之二項方程式展開式中某一項的係數為何？這種問題，使用我的系統，就很容易解決。
範例的主旨在展示計算無限位數整數運算的程式該怎麼寫？請注意！除了必須使用 b{ . . . . … }b 把演算式子包裹起來外，沒有太多的差異。四則運算的加減乘除仍可照用，只是必須牢記住，大數只能與大數運算，不是大數，格式結構當然不同，想用，就得先進行格式轉換函數的操作，才能使用。還有，想印出結果時，不能隨便印了，只能用 Big. 來印。注意到 BASIC 式的指令都沒有受到影響了嗎？這就是我設計系統的技巧，可以還有更多的數學體系被加進來，照樣可用 BASIC 格式的程式操作。
能算很大的數字，當然也得包括能算很小的數字，紐西蘭的樂透(lottory)是用 40 個數字選出 6 個數字的方式開獎，兩塊錢買一張票，中頭獎一百萬的機率是三百八十三萬八千三百八十分之一( 1 / 383,8380 )。花兩塊錢去賭這一百萬，每一塊錢的作用效率，也就是中獎機率，實在是太低了，我從不買。若算成排列，或然率就更低，是： 1 / 27,6363,3600，約27億分之一。規矩若改成排列的結果，就更不會有人願意賭了。請輸入下列操作，自己試一試：
40 6 nCr big.
40 6 nPr big.
1000 123 nCr big.
:
\ (54)nCrAndnPr.f
\ 二項式的係數可以應用 nCr 求得
\ 用法：執行 10 5 nCr big. 可得 (1+x)^10 之第 5 項係數為 252
\ 注意！使用 n r nCr 或 n r nPr 時 r 必須 >0 且 r 不得為 0
\ 使用 n r nCrOrg 則 r 可以 = 0
\ 定義 :: nCr = n! / ( (n-r)! * r! ) :: nPr = n! / (n-r)!
1 bigvariable n! 20000 allot
1 bigvariable (n-r)! 20000 allot
1 bigvariable r! 20000 allot
1 bigvariable d 20000 allot
3 integers k n r
: nCrOrg ( n r -- addr ) \ get addr of big d, r>=0
[[ r ]] ! [[ n ]] !
basic
10 LET b{ n! = big1 }b :: b{ (n-r)! = big1 }b
:: b{ r! = big1 }b :: b{ d = big1 }b
20 FOR k = 1 TO n
30 LET b{ n! = n! * i>big ( k ) }b
40 NEXT k
50 FOR k = 1 TO n - r
60 let b{ (n-r)! = (n-r)! * i>big ( k ) }b
70 next k
80 for k = 1 to r
90 let b{ r! = r! * i>big ( k ) }b
100 next k
110 let b{ d = n! / ( (n-r)! * r! ) }b
\ 120 run d big.
130 end
d \ big.
;
: nCr ( n r -- addr ) \ get addr of big d, r>0, r=/=0
[[ r ]] ! [[ n ]] !
basic
10 LET b{ n! = big1 }b :: b{ r! = big1 }b :: b{ d = big1 }b
20 FOR k = n - r + 1 TO n
30 LET b{ n! = n! * i>big ( k ) }b
40 NEXT k
50 FOR k = 1 TO r
60 let b{ r! = r! * i>big ( k ) }b
70 next k
80 let b{ d = n! / r! }b
\ 90 run d big.
100 end
d \ big.
;
: nPr ( n r -- addr ) \ get addr of big d, r>0, r=/=0
[[ r ]] ! [[ n ]] !
basic
10 let b{ d = big1 }b
20 for k = n - r + 1 to n
30 let b{ d = d * i>big ( k ) }b
40 next k
\ 50 run d big.
60 end
d \ big.
;
variable m
: test ( n -- )
m !
page
m @ 1+ 0
do
cr
." Term ( " m @ . ." , " I . ." ):"
m @ I nCrOrg big.
cr
loop
;
\ Usage: 100 test
\S
ok
10 5 nCr big.
3 digits
252 ok
49 6 nCr big.
8 digits
13983816 ok
200 100 nCr big.
59 digits
90548514656103281165404177077484163874504589675413
336841320 ok
10000 5000 nCr big.
3009 digits
15917902635324389483375972736415211886530058374576
14550428319103517772637120095798663262853944222217
74335859829932262055804632908708020739850879872195
95848962041757866458580184099587512068914331597813
53174051453473199670521394502538477277336008312053
78448823951274321755502883180927364644281795459349
36890023546288054736628292721322091972680306215783
97698552486834508478688949946112620233602352989894
58928488427591110374321646235202929095545845304023
49292778714312397841036259290830007542173305536549
24253683062815307296533408892556506908751506476159
44622376204326852230062678211259375951657115342848
24533318106868409528400428469950435925781799643074
13894226494475866262818621837575412803625468813885
44759125956185871468454381861463662350728468211441
65546574399328400579417002212869168618937974722788
62022397883728976020496710189761906178593058261688
08117556117796960379809282174855477301204105813490
54627159851188661377744154110563694305682072524481
94310502564874945796288376042950798729141780053010
24149340722579759834860211640098545723183096418633
68889831214559707246945445665178908193538606256602
93683165225062715958242340375627937873328871136143
52737971292965638066368798136853809235306441396478
97981427998980441958797431047888940127197101544121
68400963446529395285243067100038066963076992572201
04426311836533049067512198270012436774453339363870
02281179253561881400957197317504497933395227608620
35738939329776832343771264615030169561499601195082
06705891127875644018328002477885570580594271739655
61724727970366569861808080196554123575656465556543
39707955136421179968234829408914932867170470389361
58996297545140449708716896119990505242038078767450
45086398524630406716702026949125606462058300176130
06222847575106625661061937714355872185378096200269
13816305961756296827876710659465040754767228071475
82168701916632425820168589328145494184963321901025
03263315943618316059553444266801897513519884512933
06946591872301020473208721181284611163964165765568
93394074097665692587872816840693520731443017872513
61780157927471147290158311709071711945782984829441
64359840658473384707719418659651955333974514346503
81761619761261615704035455946677454877741276547147
86635414188001119626029573352659456865843697213096
86983612640564990207924247805354140963069566603071
19593156917262680235151820878651554693737963876050
46437155479530978766508167970001769266592869187570
94175117347665748132703540903393455409827319346571
30920200412827961158882797284732350179796997256267
19728263470177566063313040160755515205233184045927
56797612244679324194846919392918520452394577675953
32686906744319279375609565885643212422852240351665
84543197040090546963296363638177915596412050056857
02690372838060388519713403611629040056633420468941
76159382456860877054526939045603888375597321562922
27663423267910309912054892793591354641456968021307
92488795541350742383065293811197486421347908348956
55794152699977683783414705903919974789150191636363
96775919453875351801524980522104507017055088380935
44209022455222930021060372371375638589078163387440
553649120 ok
第(55)個範例，是前一(54)範例之延伸應用舉例，但只用到統計學上的 nCr 算出高階二項式展開後的巴斯卡三角塔式係數。
這個範例可以令您首先感到震撼的現象，應該是怎麼那麼一點點的程式就設計出了 nCr ？是的， Forth 就是可以寫得這麼精簡。
如果您還想知道 nPr 又是可以怎樣設計得出來？也是一樣，一列程式就能解決。這個話題，我就不在此增加列示了，反而應該告訴大家，最可靠資源的來處，日後您才可以大用。
Forth 的科學程式庫第 #59 號公益程式，就提供排列與組合函數的源程式，請自行網上拜訪取得。該網頁常換服務位置，是一個叫三合星(taygeta)的公司對全世界服務的網頁，我比較熟，所以只是偶爾去拜訪一下。
即使如此，提醒大家，這個新設計的 nCr 有其適用上限，算得的數據，只能適用於 64 位元整數範圍以下，前一範例沒有上限，這是兩者的差異。還有， #59 號科學程式庫提供的公益程式，排列與組合函數的輸出，都被安排成雙整數，適用範圍就能更為寬廣。我們使用這些資源時，首該注意這些應用條件，用起來才不會搞錯。
您載入程式後，就試昨天所給過的直接操作的例子，看看結果。
40 6 nCr .
還是可以得到正確的解答： 3838380 ，對吧？試大的輸入數字就不行了。
至於計算出三角塔式的係數，是死算，是把固定算法寫成程式的範例練習。
:
這個範例可以令您首先感到震撼的現象，應該是怎麼那麼一點點的程式就設計出了 nCr ？是的， Forth 就是可以寫得這麼精簡。
如果您還想知道 nPr 又是可以怎樣設計得出來？也是一樣，一列程式就能解決。這個話題，我就不在此增加列示了，反而應該告訴大家，最可靠資源的來處，日後您才可以大用。
Forth 的科學程式庫第 #59 號公益程式，就提供排列與組合函數的源程式，請自行網上拜訪取得。該網頁常換服務位置，是一個叫三合星(taygeta)的公司對全世界服務的網頁，我比較熟，所以只是偶爾去拜訪一下。
即使如此，提醒大家，這個新設計的 nCr 有其適用上限，算得的數據，只能適用於 64 位元整數範圍以下，前一範例沒有上限，這是兩者的差異。還有， #59 號科學程式庫提供的公益程式，排列與組合函數的輸出，都被安排成雙整數，適用範圍就能更為寬廣。我們使用這些資源時，首該注意這些應用條件，用起來才不會搞錯。
您載入程式後，就試昨天所給過的直接操作的例子，看看結果。
40 6 nCr .
還是可以得到正確的解答： 3838380 ，對吧？試大的輸入數字就不行了。
至於計算出三角塔式的係數，是死算，是把固定算法寫成程式的範例練習。
:
\ (55)BinomialTable.f
: nCr ( n r -- nCr )
1 swap 0 ?do over i - i 1+ */ loop nip ;
20 value m
: onerows ( row -- )
CR
DUP 4 .R
m 1+ 0
DO dup I nCr
dup 0=
if drop else 7 .R then
LOOP
DROP ;
: Table ( -- )
CR CR 4 spaces
m 1+ 0
DO I 7 .R LOOP ( display column numbers )
m 1+ 1
DO I onerows
LOOP
;
table
訂閱：
張貼留言 (Atom)
沒有留言:
張貼留言