FORTH : 慶的世界: 9月 2024

2024年9月30日星期一

一百個例題 (76 ~ 80)

Ching-Tang Tseng

Hamilton, New Zealand

1 October 2024

第(76)個範例，記錄本人自創系統中浮點數字開平方(FSQRT)的發展歷程。

回顧已介紹過的範例，範例(12)是整數的開平方(SQRT)程式，範例(32)是複數的開平方(ZSQRT)程式。相關問題均已在說明中解釋過。浮點數的開平方(FSQRT)程式又與上述兩者完全不同了，只因這些數字在系統中的資料結構完全不同，程式當然就會不同，除了不同，還有許多特質彼此也不相同，算出結果是一回事，了解它們各自的本質，則是另外一回事，全都不可混為一談。

我自創的浮點系統是十進制式的浮點系統，所有的函數本可僅由整數計算而成，沒錯，確實是這樣，我設計的函數都可以只用整數的四則運算處理出來，可是，當我最後檢視成果時，發現只有開平方不符合要求。因為，誰都知道，十八位數的整數開完平方後，只能得到九位數的整數，要這樣用，緊接其後的另外九位數只好填0了，這怎麼可以？所以，在 64 位元的時代，用這種方法設計浮點數系統中的開平方函數時，您只好用 128 位元的整數來開出 64 位元的結果，程式會非常離譜， 40 位數的整數運算，等同於處理天文數字，光算一個，就能耗掉您難以忍受的時間。

除此之外，還有其他問題，在較少位數環境能用的快速收斂算法，在多位數環境中都會凸顯出問題來，通常在首尾的邊界都不適用了。快速收斂的特殊規則，由於是利用高變率的特點來演算，就納入了易於循環的隱憂，在某些特殊數字點上就會循環而不收斂。

基於上述的困擾，我自創的系統在最後階段才解決浮點數開平方的問題，深知浮點數的開平方，絕對不宜再用雙整數的開平方方式來完成設計。可是，全世界根本沒有可供參考的文獻能提供現成的浮點數開平方程式資源，全都是依靠硬體數學計算處理器直接算出來的。我長期受此困擾，苦思不得其解，後來，為了準備一篇二分法(Bisection method)求解任何方程式實數之根，強調可以絕對都收斂的技術文章時，發現了這套永遠可靠的技術，可以用它來解決浮點數開平方函數的設計問題。於是動手做了大量的研究與發展，最後落定了設計，固定在我自創的系統之中。這套辦法能將任意位數的浮點數，開平方後得出相同位數的結果。

關於二分法的詳細說明，請自行參考 20180416 我的個人網頁英文網文：Bisection method to solve equation 。

發展這套技術的過程，是全用我自創系統完成的，剛開始時，我只寫 BASIC 式程式，解得答案，發展速度才快，確定沒有問題後，我才改寫成類似第(75)個範例中的那種只有中算符的 Forth 程式格式，就是這個範例中起始處的格式。接著，我再精簡使用宣告變數的使用數量，繼續改寫，最後改成只用純粹 Forth 寫成的格式，得到了可以永久有用的結果，才轉進到我的自創系統中使用。至此，工作並未結束，因為我自創系統的浮點數資料格式結構與眾不同，我還能進一步，憑此特殊結構，優化執行速度。例如：驗證兩浮點數是否為同一正負號時，我不執行實際的乘算後判定正負，改採只取正或負數的位元，直接進行邏輯運算後得到結果。我也在許多小地方進行優化處理，例如：除以 2 時，改採用乘以 0.5 ，令運算更快、更準，也可以大量節省執行時間。還有其他的優化細節，不屬於此範例之範疇者，就不多談。

二分法技術，是一種絕對收斂的數值分析技術，所有學電腦數學計算的人都應該了解，我對所有絕對收斂的技術都很有興趣，我也曾經研發成功過一套簡單型法(Simplex method)的曲線調適應用程式，也是個可以萬用無誤的技術，與這個範例的性能相當，非常珍貴。這些東西也都是我個人自創出來的方法，很實用。
:

 
\ (76)BisecFsqrt.f

variable count

\ full forth style
fvariable n
fvariable a
fvariable b
fvariable p
: test  ( f -- fsqrt )
  0 count !
  fdup n f!  a f!  0.0 E 0 b f!
  begin
  a f@  b f@ f- fabs 1. e -14 f>
\ 以10Bytes/Float來算，設定便可以改為1.e-18，就與系統計算結果完全一致，
\ 否則只能用1.e-14。
  while
  a f@ b f@ f+ 2.0 e 0 f/ p f!
  a f@ fdup F* n f@ f-  	\ f(a)=a*a-n
  fdup f0= if fdrop a f@ exit then
  p f@ fdup f* n f@ f-  	\ f(p)=p*p-n
  fdup f0= if fdrop p f@ exit then
  f*  f0>
  if p f@ a f! else p f@ b f! then
  1 count +!
  repeat
  a f@ b f@ f+ 2. e 0 f/
  dup -16 < if fround then 
;

\ 18 sigdigits !
: main
  101 1
  do
    cr i . ." count  = " count @ .
    cr i . ."  fsqrt = " i s>f fsqrt fs.
    cr i . ." bfsqrt = " i s>f test  fs.

  loop ;
main

\s
((
9 reals n a b p fa fb fp eps ft

: test1   ( f -- fsqrt )
  0 count !
  fdup {{ n }} f!  {{ a }} f!  0.E0 {{ b }} f!
  begin
  {{ fa = a * a - n }}
  {{ fb = b * b - n }}
  {{ eps = abs ( a - b ) - 1.E-15 }}
  eps f0>
  while
  {{ p = ( a + b ) / 2 }}
  {{ fp = p * p - n }}
  {{ ft = fa * fp }}
  ft f0>
A  if {{ a = p }} else {{ b = p }} then
  1 count +!
  repeat
  cr ." count = " count @ .
  cr ." fsqrt1 = " {{ ft = ( a + b ) / 2 }} ft fs.
  cr ." fsqrt2 = " n fsqrt fs. ;

\ need only n a b p ft
: test2
  0 count !
  fdup {{ n }} f!  {{ a }} f!  0.E0 {{ b }} f!
  begin
  {{ ft = abs ( a - b ) - 1.E-15 }}
  ft f0>
  while
  {{ p = ( a + b ) / 2 }}
  {{ ft = ( ( a * a ) - n ) * ( ( p * p ) - n ) }}
  ft f0>
  if {{ a = p }} else {{ b = p }} then
  1 count +!
  repeat
  cr ." count = " count @ .
  cr ." fsqrt1 = " {{ ft = ( a + b ) / 2 }} ft fs.
  cr ." fsqrt2 = " n fsqrt fs. ;

\ need only n a b p
: test3
  0 count !
  fdup {{ n }} f!  {{ a }} f!  0.E0 {{ b }} f!
  begin
  a b f- fabs 1.e-15 f>
  while
  {{ p = ( a + b ) / 2 }}
  a a F* n f- p p f* n f- f* f0>
  if {{ a = p }} else {{ b = p }} then
  1 count +!
  repeat
  cr ." count = " count @ .
  cr ." fsqrt1 = "  a b f+ 2.e0 f/ fs.
  cr ." fsqrt2 = " n fsqrt fs. ;
))

第(77)個範例也是我個人的創作，以最省記憶體的方式，教您如何存放暫態數據。

前面第(35)個範例，開始談到可以把暫態數據放到那裡去的問題，那邊講的就是一種方法，跟今天打算介紹的另一套新辦法不同。有什麼不同？自己跑一跑程式，就能體會，都能達到目的，用法都有限制，性質都有差異。

在範例(35)中的設計方法是先規規矩矩的宣告出資料結構，用完就忘掉它們的名稱，程式可以照跑，用過的名稱就可以重新再被使用，系統不會亂掉。今天要講的方法，是根本不給資料結構的名稱，只在用到的指令之前，宣告出必須控留的記憶體位址量，然後單憑我設計的萬用指令，叫作戴上了『帽子』(HAT)，來指揮各個暫態變數，能被定出位址，就能使用了。這樣用，是把記憶體省到不能再省的地步來存放暫態數據。此前，全世界沒見過有人像我這樣設計系統過，所以說這個 HAT 的用法，算是獨創。

我在走進 64 位元時代之後，才深深體會到 Forth 系統每當產生一個指令或資料結構單元時，頭部耗用記憶體的數量，至少都超過幾十個位元組(bytes)，非常驚人，而用來放置資料的體部，浮點數仍然也只需兩個位元組就夠了。所以，我才會想要這樣子存放暫態數據。

前曾提過，局部變數的設計，更耗系統內的記憶體，少則論百個位元組來算，多了也可以是成千上萬，若要論節省，今天講的這套方法可能最省，想用幾個？就只先行宣告出該配置(ALLOT)幾個就行了，用時，就靠執行 n HAT 來取得位址，直接了當。請注意，這個 HAT 的設計，會隨著系統的不同而不同，可能每套系統都不一樣。所以，如果一個人接觸新系統時，不知道或無法知道系統的結構，是不可能運用這種設計的。我對後來推廣或販售系統的公司抱著不讓使用者深究系統結構的觀念來使用 Forth 不以為然。要那樣用 Forth ，不如不用，用了就難有創意，無可發揮，必定會妨礙 Forth 的未來發展。

此範例中有三套系統都用到這項觀念，經我測試，都能成功。我從不為任何商售系統研發、推廣或測試，它們能用，我也不用，這是我的原則，就算我也設計出過那些系統的 ABC Forth 版本，我也絕不推出。我不接收他人賺錢，我卻免費為人服務的現象。所以，此範例中只列三個不要錢的可用程式： Win32Forth , Lina64 , gForth64 。

範例從 HAT 的設計開始，舉例，先行配置出記憶體的方式，再在緊接著的指令中，使用如此產生的記憶體，可以執行暫態存取，完成所設計的程式的工作。執行後，就能看到結果，從頭到尾一氣呵成，不難了解，以後設計程式時，這樣用也不會有問題。

但是，這套用法，也有它先天上的使用限制，就是啟動系統後載入自己的設計時，這樣用不會有問題。

如果是自己在設計系統，而且要把系統固定下來，供後續使用，在自己設計的系統中，就不能這樣運用這套設計。

問題牽扯到系統在作業系統中是絕對選址還是相對選址而被執行的狀況上去了，沒有頭部名稱的暫態數據儲存位址，在以絕對選址方式於作業系統中被執行的狀況，這套辦法不能使用。這部分觀念，解釋起來比較麻煩，就此打住。

如果不求節省，我個人建議，還是規規矩矩宣告出變數名稱來存取資料，才是正途。這些包羅萬象的使用範例，目的不是在教大家搞花俏，反倒是適合用來深刻的了解系統內部的結構。
:

 
\ (77)hat.f

\ Make an exclusive local array beside a word
: hat ( n -- addr )
  current @ swap cells - ;

\ Usage and test demo

\ Integer: compute 3*n1*n1+10*n1-7*n2
2 cells allot
: test ( n1 n2 -- )
\ n2 to 2 hat   n1 to 1 hat
  2 hat !   1 hat !
  3   1 hat  @  dup * *    dup cr .
  10  1 hat  @  *        + dup cr .
  7   2 hat  @  *        -     cr .
;

1 cells allot
: test2 ( n1 -- n2 )
  1 hat !
  1 hat @ dup * 2 *
  1 hat @ -5 * -
  10 + cr . ;

\ Floatiing point: compute torque according to power and rotating speed
\ ª«²zžÑµª p.183 , torque = hp / rpm
\ 2 cells/fp

8 cells allot
: torque ( hp rpm -- )
  6 hat f!   3 hat f!                           \ rpm to 4 hat,  hp to 2 hat
  3 hat f@   550.e0    f*                       \ hp  --> ft-lb/sec
  6 hat f@   2.e0 fpi  f*   f*   60.e0 f/       \ rpm --> rad/sec
  f/
  cr ." torque = " fs. space ." lb-ft " ;

4 cells allot
: torque2 ( hp rpm -- )
  550.e0   3 hat f!                             \ convert factore to 2 hat
  2.e0 fpi F*  f*  60.e0 f/                     \ rpm --> rad/sec
  fswap 3 hat f@ f* fswap f/
  cr ." torque = " fs. space ." lb-ft " ;

\S
\ for wina32 and lina64
\ (6). hat, fpi
\ **************************************
\ Make an exclusive local array beside a word
: hat ( n -- addr )
  current @ >nfa @
  swap 1+ cells - ;

\ testing started from this point
\ 1 cells allot
\ : tt1 ( n -- )
\   1 hat ! 10 30 * 5 + 1 hat @ /mod
\   cr ." This is tt1." 3 spaces . . ;

\ : tt2 ( -- )
\   cr ." This is tt2." ;
\ : tt3 ( -- )
\   cr ." This is tt3." ;

\ 2 cells allot
\ : tt4 ( n1 n2 -- )
\   2 hat ! 1 hat !
\   10 1 hat @ * 2 hat @ +
\   cr ." This is tt4." 3 spaces . ;


\ 2+1-->3, 5+1-->6, 3 cells/fp
\ 6 cells allot
\ : tt5 ( r1 r2 -- )
\   5 hat f!
\   2 hat f!
\   2 hat f@  1.0 E 0 f*
\   5 hat f@  3.0 E 0 f*
\   f+ cr fs. ;

3.14159265358979328 E 0 fconstant fpi

\ 6 cells allot
\ : tt6 ( hp rpm -- )
\   5 hat f!   2 hat f!                        \ rpm to 5 hat,  hp to 2 hat
\   2 hat f@   550. E 0    f*                  \ hp  --> ft-lb/sec
\   5 hat f@   fpi fdup f+   f*   60. E 0 f/   \ rpm --> rad/sec
\   f/
\   cr ." torque = " fs. space ." lb-ft " ;

\ 3 cells allot
\ : tt7 ( r -- )
\   2 hat f!
\   10. E 0 fdup f*
\   2 hat f@
\   f/ cr fs. ;

\ testing finished here
\ **************************************

\ for gforth64
\ : hat ( n -- addr )
\   current @ >body swap 3 + cells - ;

\ 3 cells allot
\ : tt8 ( -- )
\   1 1 hat ! 2 2 hat ! 3 3 hat !
\   1 hat @ 2 hat @ +
\   3 hat @ +
\   cr . ;
\ execute tt8 get 6

第(78)個範例是一個能除出無限位數，並顯示出結果的特殊程式。

這種程式，此前，網上沒有，只有我設計過，並給予妥善的顯示方式，要幾位數就能有幾位數，要看第幾位數就能直接看到該位數，整個程式是下過一番功夫才設計出來的，相關指令，我已全部固定在我設計的 ABC Forth 系統中，可以直接使用。想印出幾位數時，例如 500 位數，就先執行

500 MaxDigits !

想得到結果時，例如想知道祖沖之的密率 355/113＝？時，就執行

355 113 /..

便能得到結果。

程式的複雜度可能令您一時看不懂它的設計，沒關係，留著，改天想用時，直接用，改天有空想學時，再回頭研究我提供的這個源程式。

四則運算是一切數學計算的基礎，仔細研究四則運算，看似非常脆弱，使用過 Forth 的人，都會覺得四則運算簡直不堪一擊。因為，輸入的數字，若超過能表現的範圍時，系統完全不管。四則運算的結果，超出單整數能表現的範圍時，系統也完全不管。甚至於連警告都沒有，系統還是照樣的運行下去。

早年，在低位元的時代，這些問題確實是個嚴重問題，設計程式時，不得不自己小心。等到高位元的時代來臨後，以 32 位元為分界點，上述問題就不那麼嚴重了，因為人類的活動環境中，很少有需要 10 位數以上的四則運算。也就是在這個時候，我也才開始有興趣自己設計數學運算系統，造出來的成果才有實用價值。我這樣子使用自創的系統，已經超過十五年。

如果把能夠輸入的數字範圍，限定成只能為單整數，以此為大前題來討論四則運算。加、減、乘三則運算，都能得到絕對精確的運算結果，單整數不夠用來表現結果時，可以透過雙整數的處理來解決問題。只有除法運算，藉用雙整數，也不足以用來解決問題。例如：我就問您 355 除以 113 是個幾位數的循環小數？只用雙整數，解決不了問題。

1500 年前，中國的古代數學家進行過這方面的研究，他們沒有電腦可用。今天，不一樣了，有了電腦，但還是沒有幾個人知道，能用電腦把這樣的數學問題解決到什麼程度？我研發 Forth ，就在發展這種能力，設計了這個程式，能用來解決常人難以回答的這個問題。

程式的設計方法很簡單，就只是利用雙整數除以單整數，得到結果後，保留餘數，再將餘數補 0，擴充成雙整數，繼續除以固定不變的除數，依次類推，執行到指定的限定位數為止。

最後，才一次把結果漂漂亮亮的印出來。不一位數一位數印出結果的原因，就是電腦系統在執行印出文字時，必須叫用作業系統來輔助完成印出功能，叫用又放掉，再叫用又再放掉，這樣執行程式，速度會很慢，我實際測試過，知其問題。所以，設計成一次印出結果的方式來設計程式。

研究數字體系中關於有理數與無理數之劃分時，使用 /.. 指令可以很清楚的看出所有的循環小數無論是有多少位數才開始循環都是有理數，因為循環小數都能改寫成 P/Q 型式表示的分數，而 P 和 Q 都為整數，無理數則不行。

請使用 /.. 指令體會一下 355/113 仍是個有理數，而 pi 是個無理數，密率只是一個無理數 pi 的有理數之極佳近似值。使用 /.. 可以算出密率到底是幾位數的循環小數。另外給一個趣味數字：以 1 除以 998001 ，會得到一個 3000 位數而且很規則的循環小數，請用 /.. 指令去認識它。
:

 
\ (78) long dot /..

: (u.) ( u -- addr len )
  0 <# #s #> ;

\ Maximum digits great than 253 characters.
variable MaxDigits
400 MaxDigits !
: pad1 ( -- addr ) HERE 1024 + ; 

\ Big counted string at addr1
: BigCount ( addr1 -- addr2 len )
  dup cell + swap @ ;

\ Add string addr1 cnt to the big counted string addr2
: BigAppend  ( addr1 len addr2 -- )
  2DUP 2>R  BigCount CHARS +  SWAP CMOVE  2R>
  DUP >R @ + R> ! ;

\ Add one char to the big counted string addr
: cBigAppend  ( c addr -- )
  dup >r BigCount + c! r@ @ 1+ r> ! ;

: (/..) ( Numerator Denominator -- addr len )
  dup 0= abort" Warning: /0 error! "  
  pad1 256 0 fill
  2dup 0< swap 0< xor
  0< if 45 pad1 cBigAppend then    \ -
  abs >r abs r>
  swap over /mod
  (u.) pad1 BigAppend              \ integer part
  46 pad1 cBigAppend               \ .
  MaxDigits @ 1- 0                 \ fraction part
  do
     over >r 10 um* r> um/mod
     (u.) pad1 BigAppend
     dup 0= if leave then
  loop
  2drop pad1 BigCount ;

variable 50c/l

: RightSideMark ( -- )
  3 SPACES ." :" 
  50 50c/l +! 
  base @ decimal 
  50c/l @ . 
  base ! CR ;

: BigTypeWithRightSideMark ( addr len -- )
  0 50c/l !
  >r
  begin
        r@ 50 >
  while dup 50 type RightSideMark
        r> 50 - >r
        50 +
  repeat
  r> type ;

: 50cRuler ( -- )
  cr ." =========1=========2=========3=========4=========5"
  cr ." 12345678901234567890123456789012345678901234567890"
  cr ." ========(c) 2014 Copyright, 50 chars Ruler========"
  cr ;

: /..  ( Numerator Denominator -- )                    
  (/..) cr
  BigTypeWithRightSideMark           
  50cRuler ;

第(79)個範例是雙整數四則運算中之三則運算( D* , D/MOD , D/)的設計程式。 D+ 指令，通常在所有的 Forth 系統中，都已直接提供。

收集百例的後期，我已逐漸不耐於繼續使用 32 位元的系統，自 2012 年後，我幾乎已改為完全使用 64 位元 Forth 的境地，這個範例就是在這種情況下收集得來的，因此，我首先令其在 Lina64 系統下能夠執行，然後，回頭再修正出在 Win32Forth 系統也能執行的程式。這個程式，在兩個系統中都能使用，達到了技術能夠古今適用、接續傳承的目的。

雙整數的乘法指令，前曾討論過，這裡的收集則更進一步，讓 D* 更適應於例外狀況，如發生溢量(overflow)時的處理狀況。凡是能夠處理例外情況的指令，都得添加許多邏輯判斷後分支處理的設計，程式就會冗長、麻煩。

我在收集程式時，若能同時獲得理論依據，我都會順便收錄進範例說明中，關於雙整數的乘法與除法的理論根據式，此程式中亦有，日後不只是應用，若需轉往別的系統，增加這種同樣的設計時，就可以參考。

雙整數的除法比較麻煩，這個範例中的設計方法，使用中間過度期擴充運算範圍的方式來設計，也就是在運算過程中，採用了比雙整數更高範圍的三倍整數(Triple number)來處理運算。如此一來，就能得到比較精確的運算結果。

我自己設計過一氣呵成的雙整數除法，可以少增加指令名稱，但程式必定冗長。我承認，我的設計不如這一套收集。

研究雙整數四則運算，搞通之後，當然也能憑同樣道理，設計出四倍整數的四則運算指令，我就做過，只是練習，也只是好玩，研究到最後，就會渴望獲得一個能放諸四海而皆準的萬用程式。有的，那就是無限位數的四則運算功能，我已經將其融入我自創的系統中了，變成永遠可以使用的工具，前面已介紹過使用的範例。

雙整數四則運算的測試方法很簡單，您只需牢記所有雙整數的輸入方式，就是在數字後面多添加一個小數點就夠了，想獲得結果時，也要牢記住，別再使用 dot( . ) ，而是該使用 D dot( D. ) 。

如果因為數字位數太大，沒有把握驗證運算結果是否正確，解決方法也很簡單，您只需使用 1 後面全是 0 ，然後最後補上一個小數點的方式輸入數字便可。雖然這樣的數字看起來很整齊，系統轉換後的內容，則是非常不整齊的花樣，顯示很整齊，只是使用者的外觀感覺，指令能處理完再輸出一個也是很整齊的結果時，就表示程式設計得很正確了。

從這個範例的實際應用中，我們可以發現， Forth 對數字的認知，實際上只需兩種就夠了，一個是不帶小數點的整數，另一個就是帶小數點的雙整數。如果能夠體認出這樣的觀念，那麼，您就很容易接受我後來自創的十進制浮點系統了，因為，我的浮點數，其假數(mantissa)部分，就處理成雙整數，其指數(exponents)部分，就處理成單整數。這樣子設計浮點系統，就得用到此範例中雙整數的四則運算。
:

 
\ (79) Double precision integer arithmetic.

\ Revamped for ciforth from Dr. Tim Hendtlass, Real Time Forth, Swinburne University of Technology, 1993, P.62
\ 20170317
\ Multiply two double precision numbers to give a double precision product. 

: d- dnegate d+ ;

\ : PICK 1+ CELLS DSP@ + @ ;

\ : ?dnegate 0< if dnegate then ;

\ With overflow check.

: UD*C ( ud1 ud2 -- ud3 )        		\ all numbers unsigned doubles
dup >r rot dup >r >r over >r     		\ put a c c b on return stack
>r swap dup >r                   		\ put a d onto return stack
um*                             		\ b*d
0 r> r> um* d+ r> r> um* d+      		\ offset 16 bits, add on a*d+b*c
0 r> r> um* d+                   		\ offset another 16 bits, add on a*c
or 0<> abort" D* overflow"       		\ check for overflow
;

\ Without overflow check.

: UD* ( ud1 ud2 -- ud3 )         		\ all numbers unsigned doubles
rot >r over >r >r over >r        		\ put c b a d on return stack
um*                              		\ b*d = part of 32 bit answer
r> r> * r> r> * + +              		\ a*d+b*c= addition to top 16 bits
;

: D* ( d1 d2 -- d3 )             		\ all numbers signed doubles
dup >r dabs 2swap dup >r dabs    		\ #s +ve, keep info to work out final sign
ud*                              		\ get 32 bit answer. Change this to ud*c to get overflow c heck
r> r> xor 0< if dnegate then
\ r> r> xor ?dnegate             		\ work out and apply final sign
;

\ More clear D*
\ : D* ( d1 d2 -- d3 )
\   >r >r
\   over r@ um*
\   rot r> * +
\   rot r> * + ;

\ Division - ( U0 * 216 +U1 ) / (V0 * 216 + V1 ) = (A0 * 216 + A1 )
\ Use fast algorithm, remainder needs an additional 32 bit multiplication and subtraction.

: T* ( ud un -- ut )             		\ Unsigned double * unsigned single = unsigned triple
dup rot um* >r >r                		\ high-part of answer to return stack
um* 0 r> r> d+                   		\ get low-part,offset 16 bits,add high-part
;

: T/ ( ut un -- ud )             		\ Unsigned triple / unsigned single = unsigned double
>r r@ um/mod swap                		\ divisor > r, divide top 16 bits, rem to top
rot 0 r@ um/mod swap             		\ combine with next 16, divide these by divisor
rot r> um/mod swap drop          		\ repeat for last 16 bits, lose final remainder
0 2swap swap d+                  		\ combine parts of answer to for final answer
;

: U*/ ( ud un1 un2 -- ud2 )      		\ ud * un1 / un2, triple intermediate product
>r t* r> t/
;

: UD/ ( U1 U0 V1 V0 -- A1 A0)    		\ Unsigned 32 bit by 32 bit divide. No remainder
dup 0=                           		\ top 16 bits of divisor = 0?
if swap t/                       		\ simple case, make it a triple, do /
else                             		\ more involved case
dup 0 1 rot 1+ um/mod >r         		\ work out scaling factor,copy to r. (I did 65536. ==> 0 1)
drop r@ t* drop >r >r            		\ scale denominator, move to return stack
dup 0 r> r> 2dup >r >r u*/ d-    		\ calculate (U-U0*W1/W0)
r> r> r> -rot nip u*/            		\ multiply by (D/W0)
nip 0                            		\ /2^16, make answer double
then
;

: D/MOD ( dn1 dn2 -- drem dquot ) 	        \ Divide two signed double numbers 
2 pick over xor >r               		\ work out sign of answer
dabs 2swap dabs 2swap            		\ convert numbers to positive
2over 2over ud/ 2dup >r >r       		\ do the division, save copy of quotient
ud* d-                                        \ calculate the remainder
r> r>
r> 0< if dnegate then
\ r> ?dnegate                   		\ retrieve answer,apply final sign
;

: D/ ( dn1 dn2 -- dquot )        		\ Divide two signed doubles, no remainder
2 pick over xor >r               		\ work out sign of answer
dabs 2swap dabs 2swap            		\ convert numbers to positive
ud/                                           \ do the division
r> 0< if dnegate then
\ r> ?dnegate                    		\ retrieve answer, apply final sign
;

第(80)個範例是一個比較複雜、需要叫用動態連結庫存程式(Dynamic-Link Library, DLL)才能完成的聲音輸出程式，同時展示樂譜的產生原理。

在範例(65)中，我們也曾利用控制聲音輸出的現成指令 Beep 來產生 Morse 無線電通訊碼之 Tik-Dak 聲音輸出的效果，那是令 Forth 系統發出聲音之最簡單的方法。這一個範例展示從 Forth 系統外部叫用動態連結程式產生 Beep 功能的方法，並進行了更進一步之設計應用，用它來產生單音音樂。

可作為動態連結程式使用的軟體，是經過特殊編譯器產生的軟體，產生的方式隨著系統不同而不同，甚至於在設計時的程式書寫方式也與一般程式語言略為不同，可供外用的指令名稱都必須按照設計規格另行宣告。編譯出來的可用軟體才能在被叫用時隨便被安放在記憶體中使用。它的產生技術超出討論範圍，這裡僅做簡單說明。

動態連結程式的用法，首先就是得把動態連結程式放到系統指定的記憶體區域內，然後才能根據安放位址找到動態連結程式中原本賦予的可用指令，這個指令的名稱通常是大小寫都有所區分，根據它，就能找到它的起始執行位址。使用時功能指令可以繼續延用同一名稱，也能將其轉換成使用者自己的新命名。下面這三列程式就在執行出這樣的事情：

"kernel32.dll" LOAD-DLL: K32
"Beep" 'k32 DLL-ADDRESS: BBB
: sing ( time frequency -- )
BBB call drop ;

--------------------------------

範例(80)原本包括了五個連續的程式，但都不能直接載入 Win32Forth 系統，因為，它們都是百例收集工作中，晚期的作品。我已長期脫離 32 位元的環境，改用 64 位元所有可以免費取得的 Forth 系統。這幾年來，我大約用過不下 10 種 64 位元的 Forth 系統，最後覺得，還是 ciForth 最好，前後連貫、一脈相承、永久有效。

然而，這幾個範例，只有範例(80-1)才是 ciForth 使用的程式，它叫 Wina64 Forth 。後續的幾個，就都不是了，我看了程式最前面的內容，也想不起來是該給那一套 64 位元 Forth 系統使用的測試程式了。我曾用過裝在 hp 手提電腦上面 64 位元的 W7 OS ，這些程式，就是那時的產品。那時，我發展的軟體，都只放在外加硬碟上。後來，由於硬碟用了五年以上，轉壞了，我就改裝 SSD 取代硬碟，也不再裝 W7 OS 了。這些留下來的程式，就成了古董，但因係辛苦測試而成，捨不得丟掉，就當範例留著。我已不打算再回頭去恢復那些 64 位元的 Forth 系統，如果強要說明，我只能用猜的，可能是一個叫做 FasmForth64 的系統，我對 Fasm 也很熟，知道怎麼用。

若想把這個程式恢復成 Win32Forth 系統可以使用的程式，方法很簡單，只有兩個要點。

首先，請注意在 WIN32Forth 系統中特有的設計，如何將一長段程式，處理成只算是說明或註解的方法？就是使用一對 (( ..... )) 指令，用的方式，只需如下列所示便可：

((
********
*************
(這一區段內，可以是任何說明文字或程式)
*************
********
))

另外，請找到這套範例中的關鍵指令 sing ，把它以上的部份，全用上列雙括弧把它處理掉，改成下列新添加的設計於此處，請特別注意！在這個指令設計中之 CALL Beep ，大小寫必須嚴格區分，否則叫不到靜態連結程式庫中的功能指令。

: sing ( duration freqency -- )
CALL Beep DROP ;

只做這兩步操作就夠了。重新載入程式，您就能試出結果。

範例(80-5)就是曾被修改成在 W10 可以執行的程式，這裡不再加貼，執行程式中的幾條樂曲，都能演唱出單音音樂。

----------------------------

我設計的音樂程式與眾不同，節奏與旋律都用音樂術語，單音就不帶和聲。要設計，我也設計得出來，方法就是用 Forth 控制好幾個分開但樂理上相關的樂譜檔案，同時在一個喇叭上演出單音，就能達到和聲的效果。這是音樂方面的題外話，我不想多談。現代作業系統都使用先佔式多工(temptative mutitask)原理來運作在作業系統中所有被執行的軟體，Forth 系統能同時操作出許多個可執行程式。

上述說明能如此精簡改成的理由，是 Win32Forth 系統把很多叫用動態連結庫存程式才能完成的功能，都以叫用靜態連結庫存程式的方式固定在系統中了，使用者才能這樣使用牽扯到涉及外部影、音、圖、字輸出的硬體驅動指令。動態連結與靜態連結使用的觀念與所處的地位，可以用一張簡圖來表示，我在2019年1月2日的網文『組合語言在電腦中產生機器碼的程序』貼出過，文末再貼一次以加深印象、促進認知。

近幾年才陸續釋出的 64 位元 Forth 公益系統，在發展初期，都不具有這方面比較高階的應用功能，我領先使用，就得自己發展出許多特出的使用方法。單純使用 Win32Forth 時，不能體會出如何叫用動態連結庫存程式中之功能指令來使用的方法，方法也是各個系統都不相同，這些技術也與 Forth 本身無關。古今的做法也完全不相同。這些Forth對外連結軟體以產生功能的技術，我跟用了40年，沒有停過，它們涉及的知識非常廣泛，無法在單一個程式中交代清楚，只能用到時順便解釋一下。

留下這些範例的目的，不只是只能用來教學，也具有協助我記憶的功效，每逢我再度碰到新的 Forth 系統時，我也仍然需要回顧這些程式的內容來依樣畫葫蘆，這等同於教學相長。這種範例程式，別人也無法提供給您，至少，此前網上我未曾見過。

動態連結程式與靜態連結程式對照關係之貼圖：
:

 
\ (80-1)play music in wina forth                                  
                                                  
: make-constant ( n adr -- )                                                
  BODY> >R    
  R@ >DFA  !   
  'BL >CFA @   
  R> >CFA ! ;           
                                                   
: LOAD-DLL: ( sc -- u ) 
  CREATE $, DROP 
  DOES>  DUP >R $@  LOAD-DLL           
         DUP R> make-constant ;                                     
                                               
: DLL-ADDRESS:  ( sc xt -- adr ) 
  CREATE , $, DROP     
  DOES>  DUP >R   CELL+ $@  
         R@ @ EXECUTE DLL-ADDRESS   
         DUP R> make-constant ;          

"kernel32.dll" LOAD-DLL: K32 
"Beep" 'k32 DLL-ADDRESS: BBB
: sing ( time frequency -- )
  BBB call drop ;

\ (64)tone.f
\ 20160402
\ 基本原理只利用給音長及音頻後叫用Beep功能程式便可唱出單音來
\ 試用後效果顯示在XP環境中尚可，在W7環境中尾音被切除得很清楚，效果變差。
\ 選用192此值，係因可作為拍子設定之最小公倍數使用，以便都能除得整數。
\ 192 = 2*3*4*8

192 3 * VALUE Duration

\ n1:Duration   n2:Tone Frequency
\ : sing ( n1 n2 -- ) CALL Beep drop ;

: silence ( n1 -- ) MS ;
: | ; immediate

\ Get duration ( -- n1 )
: 4T   Duration 4 * ;
: 3T   Duration 3 * ;
: 2T   Duration 2 * ;
: 1T   Duration     ;
: T/2  Duration 2 / ;
: T/3  Duration 3 / ;
: T/4  Duration 4 / ;
: T/8  Duration 8 / ;
: T/16 Duration 16 / ;

\ Get tone frequency ( -- n2 )

\ 全面降八度時用，注意！每八度均成整數倍的關係。
\ : MDO 262 ;
\ : MRE 294 ;
\ : MMI 330 ;
\ : MFA 349 ;
\ : MSO 392 ;
\ : MLA 440 ;
\ : MSI 494 ;


: MDO 523 ;
: MRE 587 ;
: MMI 659 ;
: MFA 698 ;
: MSO 784 ;
: MLA 880 ;
: MSI 988 ;

\ Sing it ( n1 -- )
: Dl MDo sing ;   : LDo MDo 2 / sing ;   : HDo MDo 2 * sing ;
: Re MRe sing ;   : LRe MRe 2 / sing ;   : HRe MRe 2 * sing ;
: Mi MMi sing ;   : LMi MMi 2 / sing ;   : HMi MMi 2 * sing ;
: Fa MFa sing ;   : LFa MFa 2 / sing ;   : HFa MFa 2 * sing ;
: So MSo sing ;   : LSo MSo 2 / sing ;   : HSo MSo 2 * sing ;
: La MLa sing ;   : LLa MLa 2 / sing ;   : HLa MLa 2 * sing ;
: Si MSi sing ;   : LSi MSi 2 / sing ;   : HSi MSi 2 * sing ;

: Waltzing-Matilda
  1T SO   T/2 SO  T/2 SO  1T SO   1T MI
  1T HDO  T/2 HDO T/2 SI  1T LA   1T SO
  1T SO   T/2 SO  T/2 SO  1T LA   1T SO
  1T SO   T/2 FA  T/2 MI  1T RE
  T/2 Dl T/2 RE  1T MI   1T MI   1T RE  1T RE
  T/2 Dl T/2 RE  T/2 MI T/2 Dl  T/2 LLA T/2 LSI  1T Dl
  1T LSO  T/2 Dl T/2 MI  1T SO   T/2 FA T/2 MI
  1T RE   T/2 RE  T/2 RE  1T Dl  ;
: Jasmine
  1T MI T/2 MI T/2 SO T/2 LA T/2 HDO T/2 HDO T/2 LA 1T SO T/2 SO T/2 LA 2T SO
  1T MI T/2 MI T/2 SO T/2 LA T/2 HDO T/2 HDO T/2 LA 1T SO T/2 SO T/2 LA 2T SO
  1T SO 1T SO 1T SO T/2 MI T/2 SO 1T LA 1T LA 2T SO
  1T MI T/2 RE T/2 MI 1T SO T/2 MI T/2 RE 1T Dl T/2 Dl T/2 RE 2T Dl
  T/2 MI T/2 RE T/2 Dl T/2 MI 1T T/2 + RE T/2 MI 1T SO T/2 LA T/2 HDO 2T SO
  1T RE T/2 MI T/2 SO T/2 RE T/2 MI T/2 Dl T/2 LLA 2T LSO
  1T LLA 1T Dl 1T T/2 + RE T/2 MI T/2 Dl T/2 RE T/2 Dl T/2 LLA 2T LSO ;
: I-am-a-little-bird
  1T DL 1T DL 1T DL 1T T/2 + MI T/2 RE 1T DL
  1T MI 1T MI 1T MI 1T T/2 + SO T/2 FA 1T MI
  1T SO 1T FA 1T MI 3T RE
  2T RE T/2 DL T/2 SI 1T DL 1T RE 1T MI
  2T FA T/2 MI T/2 RE 1T MI 1T FA 1T SO
  T/2 SO T/2 FA 1T MI 1T RE 2T DL ;
: Pocalicalina
  1T SO T/2 HMI T/2 HMI T/2 HRE T/2 HRE 2T HDO 3T HMI
;

: TEST 1 0
  DO
  Waltzing-Matilda     2T silence
  Pocalicalina         2T silence
  Jasmine              2T silence
  I-am-a-little-bird   2T silence
  LOOP ;

2024年9月27日星期五

一百個例題 (71 ~ 75)

Ching-Tang Tseng

Hamilton, New Zealand

28 September 2024

第(71)個範例，收集了四套以高階定義方式設計出來的雙整數乘法指令 D* 。

百例中的後半段範例，大部份都是與我在發展系統的過程中相關事務的研究，如果想學系統設計，這些事務都得搞清楚。能有四套不同的作法實現程式，這個範例就不可能是一次或一天內落實下來的。有的時候，為了怕忘記程式執行的原理，我還得在收集之後，又回頭去查貢獻者當時所提供的說明文獻。

您若仔細研究這四套的各個設計，可以總結出 Forth 設計程式時的基本精神，就是在搞堆疊上數據的翻來覆去以理順秩序。這就是 Forth 跟所有傳統程式語言完全不一樣的地方。範例介紹進度已經超過 70% 了，更不能回頭重談 Forth 系統中最基本的堆疊操作指令。還有，為了暫時安放數據，這些程式是不是也顯現了開始投機使用可以臨時暫放數據的地方？不放堆疊，以減輕堆疊操作的負荷？這些精神，我在前面的範例中都說明過，所有設計系統發展到後來的人，都會走上同樣的研究路途。

我得提醒大家，許多 Forth 系統，在提供雙整數的四則運算指令時，只有 D+ ，就也有 D- ，因為減法只是將減數取負值後再與被減數 D+ ，就得到 D- 了。但是，大部份 Forth 系統不提供 D* 指令，更不提供 D/ 的指令，為什麼？我長期閱讀國際論壇上的討論，真的發現，只有一個理由可以解釋這種現象。就是 D* 還很容易設計， D/ 就是比較不容易設計，所以才都不提供 D/ 指令。

討論範例時主要使用的 Win32Forth 系統，沒有 D* 也沒有 D/ 指令。

能設計出 D/ 指令，就成了一個寶貝資源，我是會設計，我也將成果留進了庫存程式。方法只是一再的執行移位與相減， 128 位元作 128 次就得到結果了，算法沒有特出處，所以程式也只留供自己庫存。至於平常例行的使用，很少用到 D/ ，若非用不可，我設計的系統已經具有涵蓋範圍更為廣泛的大整數除法，使用起來根本沒有位數限制，所以，我自己設計的 D/ 指令就沒有必要再添入系統了。

Forth 系統中有一個為了提高計算精確度而設計的 */ 指令，裡面涉及的除法，只是用雙整數除以單整數得到單整數，不是雙整數除以雙整數得到雙整數。

四則運算是一切計算的基礎，我也只憑此基礎設計出所有的函數，最後形成了所有的數學體系。雙整數的運用，在低位元的時代，也就是 32 位元以下的時代，因為最多只有 10 位數精確度，可能還具有常用價值。但在進入 64 位元時代之後，一個單整數就具有將近 20 位數的精確度，再搞雙整數去得到 40 位數的精確度，就有點不太實際了。
:

 
\ Dmultiply.f

\ In FSL_UTIL the following is defined:
: d*1  ( d1 d2 -- dprod )   \ double multiply
  dup 3 pick xor >r
  dabs 2swap dabs
  udm*
  r> 0= IF exit ENDIF
  dxor 2swap dxor
  1. UD+c >r                \ ud1 ud2 carry1 -- ud1+ud2 carry2
  2swap r> UD+c drop ;

: d*2   ( d1 d2 -- dprod )
  >r >r
  over r@ um*
  rot r> * +
  rot r> * + ;

: d*3 ( d1 d2 -- d1*d2 )
  >r swap >r 2dup um* 2swap r> * swap r> * + + ;

: d*4 	( multiplicand . multiplier . -- product . )
  3 PICK * >R  TUCK * >R  UM*  R> +  R> + ;

第(72)個範例，是我自創之系統中一個典型的檔案操作源程式。

這個主題，是一個可以完全獨立來討論的範例，您無法在全世界網上的任何角落找到這樣的資源。對能搞檔案操作程式的設計者而言，一方面是知者不說，另方面是做法不同，所以網上沒有相關訊息。我公開展示範例，等於是將自己的技術傾囊相授，能否完全吸收？就只能看您自己了。

一般 Forth 系統，能給您的檔案操作指令，就是 create-file, open-file, read-file, write-file, close-file, copy-file, delete-file ..... 等等等之類的現成指令。可是，您若不照我這樣的設計來使用，您可以自己試試看，可能連一個檔案也開啟不了，還可能把某個檔案搞爛了，最後想殺也殺不掉。

我在發展這個範例程式期間，毀掉過 CentOS7 好幾次，最後才放棄了，此後便堅決改用 Ubuntu 。我家人有從事於軟體設計工作者，不接受我批評 CentOS7 不好的意見，說我這樣玩系統，當然會把系統搞壞。他也不肯幫我另裝 Ubuntu ，他認為我老是在搞病毒與駭客的玩意兒，很危險。只肯跟我寫 email 溝通，最好不要安裝我創作的任何程式，以免令他的公司受害，他就會丟掉飯碗。我也不置可否，反正不能跟他談 Forth 。

這就是為什麼程式設計人，搞到後來，免不了就得自行花很多時間，細讀許多技術文獻，搞清各種所謂的『協定』(protocol)之主因。我也是這麼過來的，否則，您不可能憑猜測，就知道可讀可寫檔案，屬性的位元花樣狀態為何？每個作業系統不同版本的規格也都不一樣，我連書也沒買過，學也沒上過，人也沒問過，就在網上長期搜索，找出、讀通、試成、使用這些資源。歸根結底，就是我自己知道該怎麼做，才做得成檔案操作之事。我有個生產鋁合金錠的大亨朋友，當年他就是這麼說我的，公司的錢不怕讓我隨便用，因為我勇於創新。

大家都已經不能去 CentOS7 跑這個範例程式了，可是，將其重新載入 Lina64 ，照樣能用。換句話說，我已經把它寫成通用格式了，是長久發展後的結果。這個範例在 Win32Forth 系統中不能使用，也不是我在 Ubuntu 環境中的落定程式，這是在 CentOS7 中發展系統時產生的可貴記錄，只適合當虛擬教材使用。

被操作檔案可以放置的地方，有很多可以安排的方式，可以放在 Forth 系統之內，也可以放在系統之外。能放在系統的浮動記憶體內，也可以放在系統的固定記憶體內。範例程式內，都留有這些安排方式的程式技巧。這又牽扯到必須對系統內部結構必須了解的知識要求，要講整個系統的記憶圖譜，有點麻煩，此處不深談。

請注意！我在歸納學問之後，發現，只要設計出這三個指令： GetFile , NewFile , SaveFile 。就夠用來寫任何具有檔案操作需求的程式了。您想執行 delet , copy , move ..... 等等額外要求時，別在程式內搞，請自己用滑鼠去點開資料夾，單用點選操作來達到目的，才是合理的做法，因為『眼見為真』，看著做，就不會出差錯。想要設計出所有的功能指令，我也辦得到，讀讀書，寫程式，就能完成，但沒必要。

搞檔案，是個非常龐大的話題，您難得見到，我只講使用這三個指令就夠了的範例。現在，您只需了解這些東西，按照我在範例程式中旁白所給出的使用方法，搞自己的東西，就能全用程式來解決問題。改天，我也會走向上網搞東西的發展，專搞遠方的對象。請特別注意！它們也必還是只能以檔案作為一個最基本的單元，才能成為被遠方來客處理的東西。

此範例雖在 Win32Forth 系統中不能跑，您若有心，使用 see 指令去看我在該系統中的設計，其內容，仍然是這個範例的設計，我的檔案操作設計，基本精神只有一個，永遠也不想再改。至於不同作業系統會產生不同處理現象的問題，不是重要問題。在微軟的作業系統環境中，可以產生虛擬檔案，用完就不見了。在 Linux 公益作業系統環境中，只要產生，就絕對不是虛擬檔案。這種問題，並不重要，重要的是，要知道檔案到底存放在那理？

我用這些設計做過許多事情，例如：

我的系統能使用語音輸出技術，自動讀出文字檔案的內容。
我的系統能動用好幾個檔案，叫系統合奏音樂。
我的系統能透過網路，自行發動送出百萬封信件的攻擊訊息。
我的系統能跟任何共存於系統中的所有其他軟體進行檔案溝通，並行運作。
我的系統能跑別人的可執行檔，也能再跑許多的另一個自己。
我的系統能開啟別人禁止開啟的檔案，細看究竟。

有說不完的功能，端賴於我自己想不想用而已。

百例中有許多範例的主題，不討論數學計算，換句話說，範例程式做了許多事情，不只是數學計算，這一個範例就是與數學計算無關的程式設計。
:

 
\ FILEOP.F

\ This file operation program works in Lina64 under CentOS7
\ Author: Ching Tang Tseng
\ Date  : 20160310, Hamilton, New Zealand

WANT ALLOCATE

\ -rw0rw0rw0 = 110110110b = 438d : all are R/W enable
438 CONSTANT R/W 
VARIABLE FileID
VARIABLE Fptr
VARIABLE Frem
VARIABLE Flen
VARIABLE Fsize   10240 Fsize !

\ (1)floating Fadr
\ : Fadr PAD 4096 +  ;

\ (2)allocate Fadr ???
\ Fsize ALLOCATE DROP CONSTANT Fadr

\ (3)fixed Fadr: 100 KB below EM
\ EM HERE - . --> get 33425420 --> 33 MB free spaces
: Fadr EM 102400 - ;

: SetUpFptrFrem   ( -- )
  Fadr Fptr ! 
  Flen @ Frem ! ;
: (FILE.)   ( -- addr len )
   Fadr Flen @ ;
: FileType  ( -- )
  CR (FILE.) TYPE ;
: FileDump ( -- )
  CR (FILE.) DUMP ;

\ Beware! only a R/W attributed file can be manipulated 
\ S" Filename.f" GetFile = "Filename.f" GetFile  
: GetFile ( addr len -- )
 Fadr Fsize @ 0 FILL
 R/W OPEN-FILE
 IF CR ABORT" OPEN-FILE error?" THEN FileID !
 CR ." File ID is : " FileID @ .
 Fadr Fsize @ FileID @ READ-FILE
 IF CR ABORT" READ-FILE error?" THEN
 DUP Flen ! 
 CR . ." Bytes has been read!"
 FileID @ CLOSE-FILE
 IF CR ABORT" CLOSE-FILE error!" THEN
 SetUpFptrFrem ;

\ use PAD area create all blanks
: NewFile ( addr len n -- )
  Flen !    PAD Flen @ 32 FILL
  R/W CREATE-FILE
  IF CR ABORT" CREATE-FILE error!" THEN FileID !
  PAD Flen @ FileID @ WRITE-FILE
  IF CR ABORT" WRITE-FILE error!" THEN
  CR Flen @ . ." Bytes has been written!"
  FileID @ CLOSE-FILE
  IF CR ABORT" CLOSE-FILE error!" THEN
  SetUpFptrFrem ;

\ Beware! Flen must be set, before you WriteFile
: SaveFile ( addr len -- )
  R/W CREATE-FILE
  IF CR ABORT" CREATE-FILE error!" THEN FileID !
  CR ." FileID is: " FileID @ .
  Fadr Flen @ FileID @ WRITE-FILE
  IF CR ABORT" WRITE-FILE error!" THEN
  CR Flen @ . ." Bytes has been written!"
  FileID @ CLOSE-FILE
  IF CR ABORT" CLOSE-FILE error!" THEN ;

\ for simple testing usage:
\ S" This is a simple test." SendText>F
\ : SendText>F ( adr n -- )
\   DUP Flen ! Fadr SWAP MOVE ;

\ Frem, Fptr are to be used for other testing.

第(73)個範例，仍是用來考驗 Forth 印出資料能力的另種典範。

在很前面的範例中，早就介紹過印出數字或資料的範例了，為什麼這裡又留了一個這樣的範例？當然是有理由才留的。

在電腦系統中，數字是數字，字元是字元，兩者的意義是完全不同的。只印出數字時，程式很容易寫，只印出文字時，程式也很容易寫。想要連起來印，中間不得有空格時，程式就不好寫了。所以才留這個簡化到最精簡程度的範例程式。

問題的來源是一個希臘生物研究中心的網頁，在該國經濟崩盤後，網頁也消失了，所以貼示網址已無意義，只剩這個例題還能參考。

生物學方面的研究，最需要藉著電腦作為工具進行研究的地方，可能就是處理大量的基因圖譜資訊。人體有十萬個以上的基因序列，而基因的最基本元素，則被歸納出只有四種花樣，簡稱 AGTC ，然後，所有的基因序列，就只根據這四個基本單元重複排列而成。大家對這種東西都不熟，只聽說過基因『比對』的術語，可能比較熟悉。這樣講，大家就比較熟了。

附註：生物術語
A ：Adenine 腺嘌呤
G ：Guanine 鳥糞嘌呤
T ：Thymine 胸腺嘧啶
C ：Cytosine 胞嘧啶

在大海中撈針，也像在網路上『比對』出指定資訊，也像在十幾萬個基因序列中『比對』出指定花樣。電腦可用的地方，除了計算很快，就是『比對』速度也很快，所以，近代生物科技的研究，電腦可以幫得上很大的忙。這個範例中沒有去搞『比對』之事，只是專搞『字連數目』的印出方法。電腦的比對運作，使用者是看不見的，您除了印出來看，以了解電腦程式做得對不對，別無它法。所以，想發展這方面的程式時，首先該搞的事情，就是把『字連數』的資料印出來看。我在自行發展浮點數計算系統時，也是從這樣的起點開始做起，先準備好印出浮點數的指令，才寫系統程式，測試程式設計的對不對時，只能把浮點數直接印出來看，才能知道自己設計的程式對不對。

範例分別用純粹的 Forth 與用我自創的 BASIC 式寫法設計程式，比對輸出結果後，功能完全一致。輸出時的關鍵，在印字時，使用 emit ，印數目時，使用 dot(『.』)。

在傳統 BASIC 式語言中，沒有 emit 可用，每印一筆資料後，都會跳列，就無法完成這類工作要求。

我自創的 BASIC 中設計了一個天下無敵的 RUN 指令，這在最原始的設計人 Charles H. Moore 提供的原始程式中，也沒有這樣的設計。他的 RUN 指令，只能跑上列單一個剛剛設計好的 BASIC 式程式。這一特點，是令我能夠青出於藍的最主要部份。我相信發明人與另一位改寫成 F83 規格的程式設計者，早就看過我網上發表過的訊息，這兩位前輩深明其中的道理，但我不曾跟他們聯絡過，就讓彼此的技術都埋在心裡。

關於基因『比對』的技術，此後可能會成為一個全世界想要研究的主題。不久前，我在 Euler’s 問題網頁見到一個題目，就是從一堆成千上萬個亂七八糟毫無規則的數陣中，用程式，比對出一個例如是十位數的數字花樣來。這樣的問題，就是『比對』應用的範例，我自己寫成過程式。也自己再重新使用 Fibonancci 級數產生夠亂的大數花樣，繼續考驗程式，設計得非常成功。驗證結果時，也是把比對完成後的位址留下，再根據此址，把數字當文字印出，人工觀察程式執行後看得到結果對不對？

問題若擴充到網頁上的應用，也是這樣。我從駭客技術中，讀到破解進入他人有對談網頁可以操作的方法。就是根據網址，下載該網頁的網頁用語 HTML 程式源碼，把它存進一個檔案，然後用程式去比對出要求輸入密碼時的固定指令在那裡？程式快速找出位址，予以定址後，就能單改該處的程式結構成為不問密碼，立刻回存修改後的源碼，存回同一網址，此時，就能讓想要入侵者，沒輸入密碼，就能進入網頁，遂行操作。快速比對修改源碼結構的動作，可以寫成程式，完全自動執行，侵入速度就能很快，而且程式可以到處快速測試，能侵入時就侵入，侵入完後就恢復沒入侵前的源碼狀態，不留痕跡。有防火牆不能侵入時，就快速離去。
:

 
\ PRINTstudy.f

: main1 ( -- )
  cr
  73 65
  do
     13 1
     do
          j emit i .
     loop
     cr
  loop ;

2 integers i j

: main2 ( -- )
  basic
10 run cr
20 for j = 65 to 72
30 for i = 1 to 12
40 run j emit i .
50 next i
60 run cr
70 next j
80 end ;

\s

main1 
A1 A2 A3 A4 A5 A6 A7 A8 A9 A10 A11 A12 
B1 B2 B3 B4 B5 B6 B7 B8 B9 B10 B11 B12 
C1 C2 C3 C4 C5 C6 C7 C8 C9 C10 C11 C12 
D1 D2 D3 D4 D5 D6 D7 D8 D9 D10 D11 D12 
E1 E2 E3 E4 E5 E6 E7 E8 E9 E10 E11 E12 
F1 F2 F3 F4 F5 F6 F7 F8 F9 F10 F11 F12 
G1 G2 G3 G4 G5 G6 G7 G8 G9 G10 G11 G12 
H1 H2 H3 H4 H5 H6 H7 H8 H9 H10 H11 H12 
 ok
main2 
A1 A2 A3 A4 A5 A6 A7 A8 A9 A10 A11 A12 
B1 B2 B3 B4 B5 B6 B7 B8 B9 B10 B11 B12 
C1 C2 C3 C4 C5 C6 C7 C8 C9 C10 C11 C12 
D1 D2 D3 D4 D5 D6 D7 D8 D9 D10 D11 D12 
E1 E2 E3 E4 E5 E6 E7 E8 E9 E10 E11 E12 
F1 F2 F3 F4 F5 F6 F7 F8 F9 F10 F11 F12 
G1 G2 G3 G4 G5 G6 G7 G8 G9 G10 G11 G12 
H1 H2 H3 H4 H5 H6 H7 H8 H9 H10 H11 H12 
 ok

第(74)個範例是專講數值積分的最基礎範例，雖基礎，卻已足供實際應用。

我大學時教我數值分析的王九龍老師曾經明確的告訴我們：『數值積分可以積出很精確的結果，數值微分的結果則比較不準。』

那時，我們沒有太多的機會跑大型電腦，對老師授業的重點沒有太強烈的感受，五十幾年後的今天，不一樣了，我自己能創作系統，更對這樣的教誨充滿了無限的感激，因為，我能仔細體會得出教師名言的實際意義。我先講實例說明，再討論程式。(貼出此文的日期是教師節 9 月 28 日）

量測中子時的積分技術。

新的原子爐，自己的偵測系統無法量得讀數非常微小的中子源微少的中子通量訊號，少到什麼程度？大約是每秒只能有個位數的中子量。但新爐子就只能用這個量來啟動，邊測量邊添加燃料入爐，一直要測量到爐子能夠自持臨界反應的零功率後，我再測出可以添加多出大約 6% 的燃料量，以便控制系統能涵蓋掉這多出來的 6% 反應能力，才終止實驗，停下爐子，拆除我的實驗設備。

由於爐子在臨界之後，便可以擁有很強的輻射，能持續的激發已放在爐中的中子源，此後的重新起爐，就不必再藉助於另外加裝非常靈敏的中子量測裝置了。我記得很清楚，當時，我啟動測量裝置時，每秒的中子讀數是 2.67 個，系統自己安裝之永久性三套中子偵測裝置，必須要等到我的測量讀數為每秒超過百萬個中子以後，指針才有指示。臨界過一次之後，三套儀器就都能自行工作了。

每秒鐘 2.67 個中子，是怎麼測量出來的？只憑內含武器級鈾(鈾濃縮度大於90%)的中子偵測器，每次累積測量半小時，等到顯示裝置長出一個 Gauss 分佈狀的小饅頭專有信號後，我把曲線下方的面積積分出來，總數除以 30x60=1800 秒鐘，才得到這樣的數目。實驗過程中，我要不斷地搞信號積分，把每秒鐘的中子通量讀數算出來。

這種所謂的個位數讀數，必須調整我的偵測器與中子源的距離來完成。太近了，讀數太大，雖比較靈敏，等到臨界時，我的偵測器讀數就會超飽和，滿標而失效了，就不能決定燃料的可多添加量。太遠了，讀數太小，因為不靈敏，低讀數的測量就得進行更多次，每次加添過一次燃料後，我的測量至少就都得等上半小時，這樣搞爐子，從早上開始搞起，等到天黑了，也搞不出來。所以，我得自己調整出起始讀數就是個位數，這樣的實驗，我做了兩星期。

這就告訴了大家，積分確實可以積得很準，準到能夠搞原子爐啟動，與搞原子彈爆炸。

微分呢？不太準，不準到什麼程度？新聞上有個例子可以解釋。

波音 737 MAX 老摔飛機，最後找出了問題所在，就是機身上的攻角測量儀器讀數有問題，才亂搶控制權，讓飛機栽下來的。攻角怎麼量？量的就是機體跟水平方向的夾角，靠空速管內流體的變量，比對變率後，提供數據，這儀器肯定靠微分(變率)產生信號，結果，微分天生就做不準，電腦再好，就是算不準，就出問題了。結論是，微分的不準，可能不準到會摔飛機。

在原子爐的控制實務中也有個微分很難取得準確值的例子。

控制原子爐時，強調要把中子的倍增週期控制在 10 秒鐘左右，而週期的換算公式是 1/T = dln(N)/dt ，其中： T 是週期，N 是當時量得的中子通量，中子通量之對數標值的變率就是週期的倒數。中子通量是數位式的計量值，如果採用把每秒幾個中子的讀數曲線進行這樣的微分，肯定是來不及反應出準確的週期值來，這樣操控原子爐就很危險。實際的控制系統採用二極體的電壓對電流的對數式特性曲線來對應到中子通量，其變率是類比電腦的反應速度，便可以即時指示出中子倍增週期為幾秒。這種控制信號不能採用對中子通量之數位訊號進行數值分析的技術來獲得微分變率，它會不準也不夠快。

回頭講我們的範例程式，積分很簡單，就是用微小的寬度 dx 乘上該點的函數值 f(x) ，然後從頭到尾，把所有的微小面積加總起來，就得到結果了。很不幸，函數曲線斜率變化較大時，亂算面積，所得結果，誤差就會很大，於是，我就把光算長方形，取低位左邊長、中央長、高位右邊長的結果全算一遍來比較，另外也搞個梯形面積算法也參加進來，全都不準。最後，用到了所謂的辛普森算法，才終於有了比較好的結果，於是，範例就留在這裡了。這些算法在 Rosetta code 網頁上有詳細說明，也有很多程式可以比較，有興趣者，請自行前往。範例(74)是純用 Forth 寫成的程式，附帶方法說明。範例(74-1)則是全用 BASIC 語法寫成的程式，附帶執行結果。

積分準確的程度，甚至於可以用來作為系統的函數值產生源，問題在要細分成萬點以上的面積計算才有效果時，算得太慢了，所以才沒有人這樣使用。例如： 1/x 的積分就是對數函數 ln(x) 之值，還有好幾個函數都能這樣設計。

積分原理，也不僅只有這幾個而已，我還知道例如：辛普森 1/3 法則、 3/8 法則、韋爾德(Weddle)法則...等等，都可以用，數值分析的課本中有更多的方法可以參考，我只能點到為止的提及它們的名稱。

本網頁 20180401 的文章，是全用英文貼出介紹數值積分(Numerical Integration)的文章，其內容就是取材自這個範例。
:

 
\ (74)SimpsonIntegration.f

fvariable step
 
defer method ( fn F: x -- fn[x] )
 
: left                    execute ;
: right  step f@       f+ execute ;
: mid    step f@ 2e f/ f+ execute ;
: trap
  dup fdup  left
      fswap right f+  2e f/ ;
: simpson
  dup fdup  left
  dup fover mid 4e f* f+
      fswap right f+  6e f/ ;
 
: set-step ( n F: a b -- n F: a )
  fover f- dup 0 d>f f/ step f! ;
 
: integrate ( xt n F: a b -- F: sigma )
  set-step
  0e
  0 do
    dup fover method f+
    fswap step f@ f+ fswap
  loop
  drop fnip
  step f@ f* ;
 \ testing similar to the D example
: test
  ' is method  ' 4 -1e 2e integrate f. ;
 
: fn1  fsincos f+ ;
: fn2  fdup f* 4e f* 1e f+ 2e fswap f/ ;
 
7 set-precision
test left    fn2  \ 2.456897
test right   fn2  \ 2.245132
test mid     fn2  \ 2.496091
test trap    fn2  \ 2.351014
test simpson fn2  \ 2.447732

\S
Numerical integration
Write functions to calculate the definite integral of a function f(x)     
using all five of the following methods:

  rectangular
  left
  right
  midpoint
  trapezium
  Simpson's

Your functions should take in the upper and lower bounds (a and b),   
and the number of approximations to make in that range (n).

Assume that your example already has a function that gives values for f(x).

Simpson's method is defined by the following pseudo-code:

h := (b - a) / n
sum1 := f(a + h/2)
sum2 := 0

loop on i from 1 to (n - 1)
sum1 := sum1 + f(a + h * i + h/2)
sum2 := sum2 + f(a + h * i)

answer := (h / 6) * (f(a) + f(b) + 4*sum1 + 2*sum2)
Demonstrate your function by showing the results for:

f(x) = x3,  where x is [0,1], with 100 approximations. The exact result is 1/4, or 0.25.
f(x) = 1/x, where x is [1,100], with 1,000 approximations. The exact result is the natural log of 100, or about 4.605170
f(x) = x,   where x is [0,5000], with 5,000,000 approximations. The exact result is 12,500,000.
f(x) = x,   where x is [0,6000], with 6,000,000 approximations. The exact result is 18,000,000.

\ (74-1)RectangleIntegration.f
\ Author: Ching-Tang Tseng
\ 20161107, Hamilton, Hamilton, NZ

9 reals a b h sum1 sum2 x f(x) f(a) f(b)
2 integers n i

: function ( -- )
basic
10 let { f(x) = 1 / x }
20 end ;

: Setting ( -- )
basic
10 let n = 1000
20 let { a = 1 }
30 let { b = 100 }
40 end  ;
\ All parts above are adjustable setting.

\ All parts below are fixed.
: SimpsonInit ( -- )
basic
10 run Setting

20 let { h = ( b - a ) / I>R ( n ) }

30 let { x = a + h / 2 }
40 run function
50 let { sum1 = f(x) }

60 let { sum2 = 0 }

70 let { x = a }
80 run function
90 let { f(a) = f(x) }

100 let { x = b }
110 run function
120 let { f(b) = f(x) }

130 end ;

: Simpson ( -- )      \ 辛普森法
basic
10 run SimpsonInit

20 for i = 1 to n - 1

30 let { x = a + h * I>R ( i ) + h / 2 }
40 run function
50 let { sum1 = sum1 + f(x) }

60 let { x = a + h * I>R ( i ) }
70 run function
80 let { sum2 = sum2 + f(x) }

90 next i

100 print "       Simpson integration = " ;
    { ( h / 6 ) * ( f(a) + f(b) + 4 * sum1 + 2 * sum2 ) }

110 end ;

: L-Rectangle ( -- )    \ 長方形左側值法
basic
10 run Setting

20 let { h = ( b - a ) / I>R ( n ) }
30 let { x = a }
40 let { sum1 = 0 }

50 for i = 1 to n
60 run function
70 let { sum1 = sum1 + h * f(x) }
80 let { x = x + h }
90 next i

100 print "   L-Rectangle integration = " ; { sum1 }

110 end ;

: R-Rectangle ( -- )   \ 長方形右側值法
basic
10 run Setting

20 let { h = ( b - a ) / I>R ( n ) }
30 let { x = a }
40 let { sum1 = 0 }

50 for i = 1 to n
60 let { x = x + h }
70 run function
80 let { sum1 = sum1 + h * f(x) }
90 next i

100 print "   R-Rectangle integration = " ; { sum1 }

110 end ;

: M-Rectangle ( -- )   \ 長方形中間值法
basic
10 run Setting

20 let { h = ( b - a ) / I>R ( n ) }
30 let { x = a }
40 let { sum1 = 0 }

50 for i = 1 to n
60 let { x = x + h / 2 }
70 run function
80 let { sum1 = sum1 + h * f(x) }
90 let { x = x + h / 2 }
100 next i

110 print "   M-Rectangle integration = " ; { sum1 }

120 end ;

: Trapezoidal ( -- )   \ 梯形法
basic
10 run Setting

20 let { h = ( b - a ) / I>R ( n ) }
30 let { sum1 = 0 }

40 let { x = b }
50 run function
60 let { f(b) = f(x) }
70 let { x = a }
80 run function
90 let { f(a) = f(x) }
100 let { sum1 = h * ( f(b) + f(a) ) / 2  }

110 for i = 1 to n - 1
120 let { x = x + h }
130 run function
140 let { sum1 = sum1 + h * f(x) }
150 next i

160 print "   Trapezoidal integration = " ; { sum1 }

170 end ;

18 sigdigits !
: main ( -- )
cr cr ." Integration of 1/x from 1 to 100 = Ln(100) " cr
L-Rectangle
R-Rectangle
M-Rectangle
Trapezoidal
Simpson
cr cr ." Ln(100) = " 100.0e0 fln f. cr ;

main

\S

Integration of 1/x from 1 to 100 = Ln(100)

  L-Rectangle integration =        4.65499105751467615  
  R-Rectangle integration =        4.55698105751467615  
  M-Rectangle integration =        4.60476254867837518  
  Trapezoidal integration =        4.60598605751467615  
      Simpson integration =        4.60517038495714217

Ln(100) = 4.60517018598809137
 ok

第(75)個範例是以高斯(Gauss)法求解多元一次聯立方程式的實用程式。

我在 1980 年代就已經很熟這種專門用來求解下列矩陣方程式的技術：

[A] * [X] = [Y], 已知 [A] 及 [Y] 求 [X] 。

方法也很簡單，我都能背下來了，根據線性方程式的定義，基本運算原理就是：『任何一列乘以或除以任意常數，加或減到另一列去，整個矩陣方程式的本質不變』。

於是，像前面的範例展示過的矩陣分解法ㄧ樣，把 [A] 搞成上三角矩陣後，就能從右下角最後一個元素的直接對等關係，得到該 X(n) 的解答，然後，反推上去，次一個 X(n-1) 也有了答案，依次類推，最後，所有的 X(i) 的答案就都得到了。其中， i = 1, ... n 。

話雖如此，想寫成程式，並不簡單。一方面是運算時，後面跟著的 Y(i) 值也得跟著變。另方面是為了維持浮點數運算的最佳精確度，矩陣在操作前，最好先挑出數值為最大的元素來，將其轉置到 [A] 矩陣的左上角當樞紐元素，然後才開始演算，這樣，所得的結果才會最精確。還有，免不了有時會碰上無意義的(trivial)線性相依情況，設計程式時，也得一併解決，成果才能萬用。

講完了原理，改為強調一下這個程式的特色，才是今天介紹範例的重點。

請注意程式的寫法，不是 BASIC 語法，每一列前，不用標號，沒有 BASIC 指令。也不是純 FORTH 語法，你看不到堆疊操作指令，例如： dup, over, swap, drop ... 等等。

程式中所有的算式，都可以寫成純粹中算符(infix)的格式，只須記住，在『{{ ..... }} 』內，只對浮點數運算，在『 [[ ..... ]] 』內，只對單整數運算，記住這樣的規矩就夠了。如果另外還想進行複數的中算符運算，那就多記一個，在『{[ ..... ]} 』內，只對複數運算。

這樣的系統性能，在我提供的 ABC Forth 系統中，永遠都有，大家可以直接載入此程式來看結果，就能明白。

事實上，我早在 1980 年代，就已經發展出這種單用中算符寫程式的技術了，那時還沒有搞出 BASIC 程式語法中所需要的指令。最大的差異，就在邏輯語法， Forth 與 BASIC 兩者互相相反。固定迴路寫法，也完全不同， Forth 用 DO ..... LOOP ， BASIC 用 FOR ..... NEXT ，上下限用法完全顛倒。

這個程式改寫出來的日期，標示為 2009 年，更原始的純 Forth 格式程式，也是我親手設計的，在『FORTH期刊』第 9 期，民國七十六年(1987)年一月份出版， p.58 中就有我寫的程式，該文作者項舉勁在文末參考文獻中，提到使用我的這個設計，是純用傳統 Forth 程式語言寫成的，程式早就在許多地方試用過，都沒有問題。

我從使用純粹的 Forth 開始，於 1980 年代，發展出中算符技術，就像這個範例中所留下來的記錄一樣。 1990 年代以後，我改走 BASIC 語法技術，演變出了今天的成果，就是現在大家所見與所用的 ABC Forth 。

我自創的系統中，還有許多其他的隱藏特性，沒有用到之前，我就沒講，看到了再說。我自知，就算我介紹完ㄧ百個例題，也不可能討論得完所有我自創系統的特點。更何況，我一直還在發展，有些新東西，還都是配合現行作業系統才發展出來的新特性。這無限多種的特殊技術，是講也講不完的東西，我們只能見一個討論一個，這也是 Forth 程式語言的本質，一直都能繼續發展。
:

 
\ (75)高斯消去法解多元聯立一次方程式2009-08-22

INTEGER #ROW INTEGER ITEMP             \ #ROW為方陣量，本例為3X3矩陣
INTEGER M INTEGER M1 INTEGER M2
REAL FTEMP
10 ARRAY Y    10 ARRAY X               \ 矩陣上限宣告，暫設為10X10
10 10 MATRIX AA

: 輸入數據
  [[ #ROW = 3 ]]
  {{ AA ( 1 1 ) = 3 }} {{ AA ( 1 2 ) = 2 }} {{ AA ( 1 3 ) = 1 }}
  {{ AA ( 2 1 ) = 2 }} {{ AA ( 2 2 ) = 4 }} {{ AA ( 2 3 ) = 5 }}
  {{ AA ( 3 1 ) = 4 }} {{ AA ( 3 2 ) = 1 }} {{ AA ( 3 3 ) = 1 }}
  {{  Y ( 1 ) = 17  }} {{  Y ( 2 )  = 41 }} {{  Y ( 3 ) = 16  }}  ;

: 高斯處理
  #ROW 1
  DO
    [[ M = I + 1 ]]
    #ROW 1+ M
    DO
      #ROW 1+ M
      DO
  {{ AA ( J I ) = AA ( J I ) - ( AA ( K I ) * AA ( J K ) ) / AA ( K K ) }}
      LOOP
  {{ Y ( I ) = Y ( I ) - AA ( I J ) * Y ( J ) / AA ( J J ) }}
    LOOP
        M #ROW =
        IF  {{ FTEMP = AA ( M M ) }}
            FTEMP F0=
            IF ." Equation disparity! " QUIT THEN
        THEN
   LOOP ;

: 得到解答
  [[ M = #ROW + 1 ]]
  #ROW 1+ 1
  DO
    [[ M1 = M - I ]]
    {{ X ( M1 ) = Y ( M1 ) / AA ( M1 M1 ) }}
    [[ ITEMP = ( M - I ) - #ROW ]]
    ITEMP DUP
    0>
    IF ." Trival! " QUIT THEN
      0<
      IF #ROW 1+ 1
        DO
        [[ M2 = M - I ]]
        {{ X ( M1 ) = X ( M1 ) - AA ( M1 M2 ) * X ( M2 ) / AA ( M1 M1 ) }}
        [[ ITEMP = ( M - I - 1 ) - ( M - J ) ]]
        ITEMP
        0=
        IF LEAVE THEN
        LOOP
      THEN
  LOOP ;
: 輸出答案
  #ROW 1+ 1
  DO
   {{ FTEMP = X ( I ) }}
   CR ." X( " I . ." )= " FTEMP F.
  LOOP ;
: GAUSS
  輸入數據
  高斯處理
  得到解答
  輸出答案
;

2024年9月23日星期一

一百個例題 ( 66 ~ 70 )

Ching-Tang Tseng

Hamilton, New Zealand

24 September 2024

第(66)個範例，是兩個展示叫用自己的自用副程式(recursive)如何使用的簡例。

第一個所舉之例，可用來求解出兩個數字的最大公約數。利用輾轉相除法，互相除到 0 前的最後數字，就是最大公約數。

Forth 的程式寫法，就是這麼簡單，只需一列程式，簡單到非常容易看懂，簡單到幾乎不需要寫說明來解釋的必要。

但是，如果我不指出一些叫用自己之自用副程式寫法上該注意的要點，大家可能光看這個程式，還是很難理解。

第二個所舉之例，是以運用自用副程式的方式，設計『一個從 1 加到 n 之總和為多少？』的程式，例如：從 1 加到 10 得 55 ，如何設計這種程式？。

使用自用副程式時，該注意的關鍵要點為：必須有結束叫用自己的機制，否則程式的執行會循環不已而令系統崩潰。

這兩個例題中終結叫用自己的機制都被安排在邏輯指令組中的 if ,,,,, then 之前的 ?dup ，一旦達到終止條件的要求，就可以不再以叫用自己的方式結束程式。 if 之前的 ?dup 出現 0 時，會產生終結條件，程式便不再執行 if ..... then 之間的程式，也就結束了 recurse 。
:

 
\ (66)RecursiveGCD.f

: GCD ( n1 n2 -- n3 )
  ?DUP IF SWAP OVER MOD RECURSE THEN ;

variable TotalSum
: n-  ( n -- )
  ?dup if dup TotalSum +! 1- recurse then ;
: ss ( n -- )
  0 TotalSum !
  n-
  TotalSum @ . ;

第(67)範例是一個新冠疫情初起時用來預卜次日確診人數的指數回歸(exponential regression)分析程式。

2020 年初，新冠病毒疫情之事，吵得沸沸揚揚，網上言論無數，嚇壞了大家，搞得全台灣連口罩都被戒嚴處理，無處可買，實在好笑。我對這些事情沒有興趣，但也不甘寂寞，對涉及濫用數學理論曲解問題的胡扯言論，不以為然。這次事件給了我一套可貴的實際數據，能用於程式測試，這才是我真正感興趣的地方。平時，我若研究數學問題，很難獲得實際有意義的數據來配合測試，數據都是一長串的數字，自己亂造，肯定不太合理，官方發佈的數據，就絕對可取。收集數據是長期的工作，早收集也沒有，晚收集會消失，只能每天固定在早上起床後立刻做，做了好幾天，都是有用的東西，非常可貴。

事件剛開始時，台灣反中，政府藉機恐嚇群眾，用了一個指數函數的傳染病增長模式來恐嚇大家，照那套講法，一個月內就會有幾千萬人的災難發生。過了一個多月，中國沒有發生胡扯的結果，可見年初台灣發表恐嚇言論者的可惡。不久，台灣也染上了病毒，按照同樣的講法，一個月內全台灣都會感染，可能嗎？不可能。

本來，描訴疫情應該使用的數學模式叫 Logistic function 不是 exponential function 。但在疫情初起時，沒有真正的數據可供進行曲線調適來算出 Logistic function ，當時所能預估的事情，也只能以指數回歸算法預測次日的可能確診人數。

上述兩種數學模式的初起階段，變化近似，所以在現象尚未進行到超過數學上所謂的曲線之拐點(inflection point)以前，兩者近似，過了拐點，結論就會完全不同。我在收集了十天數據之後，開始設計這個指數函數的回歸分析程式，並用來直接預測次日結果，雖未中的，但趨勢不變，誤差不大。

可以預知之更後來的情形，就是一旦指數回歸分析預卜出來的數字逐漸高出實際確診人數越來越多時，就能確定疫情的拐點已經出現。後半部的演變，會是第三象限曲線以原點為樞紐表現於第一象限的對稱曲線，兩者合併起來，就是 Logistic function 。

指數函數的形式是這樣的：

y(t) = a * exp ( b * t ) 其中

t 表時間，以天為單位。n 為所取數據組的數量，此例為 10 組。

a 及 b 均為係數，在取得此函數式的 n 組相關數據後，藉統計學上的回歸(regression)分析技術，就能算出這兩個係數的值，相關公式如下：

b = [ n * SUM { Xi * ln ( Yi ) } - ( SUM Xi ) * { SUM ln ( Yi ) } ]
/ [ n * SUM ( Xi ^ 2 ) - ( SUM Xi ) ^ 2 ]

a = exp [ (1/n) * { SUM ( ln Yi ) - b * ( SUM Ei ) } ]

相應的程式如下：
:

 
\ (67)exponential regression

2 integers n i 
9 reals s1 s2 s3 s4 s5 x y a b
20 array xx
20 array yy


: initialize
basic
10 let { s1 = 0.0e0 }
20 let { s2 = s1 }
30 let { s3 = s1 }
40 let { s4 = s1 }
50 let { s5 = s1 }
60 end 
;

: InputData
basic
10 let n = 9
20 let  { xx ( 0 ) = 0.0e0 } :: { yy ( 0 ) =   291.0e0 }
30 let  { xx ( 1 ) = 1.0e0 } :: { yy ( 1 ) =   324.0e0 }
40 let  { xx ( 2 ) = 2.0e0 } :: { yy ( 2 ) =   571.0e0 }
50 let  { xx ( 3 ) = 3.0e0 } :: { yy ( 3 ) =   830.0e0 }
60 let  { xx ( 4 ) = 4.0e0 } :: { yy ( 4 ) =  1287.0e0 }
70 let  { xx ( 5 ) = 5.0e0 } :: { yy ( 5 ) =  1975.0e0 }
80 let  { xx ( 6 ) = 6.0e0 } :: { yy ( 6 ) =  2744.0e0 }
90 let  { xx ( 7 ) = 7.0e0 } :: { yy ( 7 ) =  4599.0e0 }
100 let { xx ( 8 ) = 8.0e0 } :: { yy ( 8 ) =  6078.0e0 }

200 end
;

: regression 
basic
10 run initialize InputData
20 for i = 0 to n - 1
30 let { x = xx ( i ) }
40 let { y = ln ( yy ( i ) ) }
50 let { s1 = s1 + x }
60 let { s2 = s2 + y }
70 let { s3 = s3 + x * x }
80 let { s4 = s4 + y * y }
90 let { s5 = s5 + x * y }
100 next i
 
110 let { b = ( i>r ( n ) * s5 - s2 * s1 ) / ( i>r ( n ) * s3 - s1 * s1 ) }
120 let { a = exp ( ( s2 - B * s1 ) / i>r ( n ) ) }
130 print " exponential regression: y = a * exp ( b * t ) "
140 run  cr ." a = "  a f. 
150 run  cr ." b = "  b f. cr

160 let { y = a * exp ( b *  8.0e0 ) }
170 run cr ." y(8) = " y  f. cr
180 let { y = a * exp ( b *  9.0e0 ) }
190 run cr ." y(9) = " y  f. cr 

200 end
;

regression

\s
程式的執行結果如下：

exponential regression: y = a * exp ( b * t ) 
a = 254.72339048
b = 0.4020241052

y(8) = 6351.0221997

y(9) = 9493.8088062

第(68)個範例，是一個能自動將數字的因數全部分解出來，改寫成質數連乘等式的程式。

程式中列示源程式的來源，是一個法國人列寧(Lehning，跟馬克斯沒有關係)，於 1982 年捐贈給巴黎數學機構的公益程式。我自創的 ABC Forth 很輕易的就能改寫完成，而且到處應用。我的個人網頁網文中，有一篇『因數分解美國總統』的諷刺性技術貼文，就用這個程式。

深究我自創系統的性能，有一項功能是一般傳統 BASIC 所不能比的。也就是在固定循環指標量非常大時，如果循環內有邏輯判斷判定出可以不必再循環下去的條件時，我的系統有辦法以非結構化的方式跳出迴路，可以不必繼續執行剩餘的天文數字循環指標量，可以大量節省執行時間。

一般人設計系統，為了迎合『結構化』的條件，就絕對無法辦到我設計的這項功能。此前，全世界除了組合語言的天生可以亂跳之功能也可以如此之外，我沒見過任何程式語言能有這種性能。

這套判斷數字是否為質數的原理很普及，但我們在學校裡的教科書全部不寫、不教，我是在大學畢業後，買了全套三本的『數學之內容方法及意義』書籍，從其中的第十章之『數論』章節內學來的，大家可能沒這些書，細節我就不寫了，只告訴大家，這是一個古希臘的數學家 Eratosthenes (西元前 276 年到西元前 194 )，所提出的『篩選』法。也就是大約就在孔子(西元前 551 年到西元前 479 年)死後沒多久， Eratosthenes 就知道怎麼搞了，他還知道如何利用太陽光線直射的角度，來測量地球的圓周有多長？頭腦真不簡單。

篩選法就是：先用 2 試，再用 3 試，如果都不能除得盡，就再用 6i-1 與 6i+1 繼續試，其中 i 是從 1 開始的整數，一直試到被測數的開平方為止，如果全都除不盡，那麼，這個數字本身就是個質數。

沒人能證明質數的總量是有限個，我們的電腦可以使用的單整數越來越大，現在 64 位元可有 18 位數，已大到非常驚人，我為了想證明這個方法找出的質數依然有效，我就寫過程式，讓我的電腦跑 18 位數整數的測量，還真有效呢！我也能寫 18 位數以上的測量程式，會讓電腦跑到冷卻風扇一直全速運轉，電腦就累死了，也熱死了，為了不把電腦搞壞，只好壓 ctrl+C 結束系統的執行了。

這個程式的道理很簡單，在許多方面也都可以繼續延伸性的使用，它不僅是留下來的範例，常被我挑出來活用，我程式中的說明也很清楚，每次使用時，都可以不花腦筋的套進我的新程式中，這就是 BASIC 式程式語言的最大好處。
:

 
\ factors.f
((
'*********************************************
'*    Factorization of an integer number     *
'* ----------------------------------------- *
'* Sample run:                               *
'*                                           *
'* ? 394616                                  *
'*                                           *
'* 2 2 2 107 461                             *
'*                                           *
'* ----------------------------------------- *
'* Ref.: "Mathmatiques par l'informatique   *
'*        individuelle (1) By H. Lehning     *
'*        et D. Jakubowicz, Masson, Paris,   *
'*        1982" [BIBLI 06].                  *
'*********************************************
))

4 integers N D I L

: test
page
basic
10 print " Enter integer number to be factorized "
20 inputi N

\ 'Test if multiple of 2
40 let D = N mod 2
50 if D = 0 then 70
60 goto 100
70 print  2
80 let N = N / 2
90 goto -40

\ 'Test if multiple of 3
100 let D = N mod 3
110 if D = 0 then 130
120 goto 200
130 print 3
140 let N = N / 3
150 goto -100

\ 'Test of divisors 6i-1 and 6i+1 up to sqr(N)
\ 'Prime numbers are of the form 6i-1 or 6i+1
200 let L = sqrt ( N ) + 1
210 for I = 6 to L step 6

220 let D = N mod ( I - 1 )
230 if D = 0 then 250
240 goto 300
250 print ( I - 1 )
260 let N = N / ( I - 1 )
270 goto -220

300 let D = N mod ( I + 1 )
310 if D = 0 then 330
320 goto 400
330 print ( I + 1 )
340 let N = N / ( I + 1 )
350 goto -300

400 next I

500 if N > 1 then 520
510 goto 600
520 print N

600 end ;

第(69)個範例是源自於他種程式語言有此用法的 Forth 等效設計，叫『局部副程式』(quotations)設計技術。

局部副程式只是我個人的譯名，取其意義而譯。

這套技術，是遲至前幾年才出現的術語，只因有些慣用其他程式語言的闖入者，拿這種術語考驗 Forth 老手，於是，在國際論壇上展開過一陣子論戰。末了，有人捐獻了這麼一個精簡設計，我錄下了這份技術。當然，取用之初， Win32Forth 是跑不通的。我看得懂它在談什麼，就依理改出了這個檔案，當作範例留參。

所謂『局部』，也就是這種用法，只能在單一個指令的設計範圍內有效，出界就不能再用，不能給別的指令設計時使用。『副程式』就是可以像巨集指令那樣，在本指令範圍內，不斷地被重複叫用，叫用幾次都可以執行。

為什麼要強調這樣的功能？

別種程式語言的程式設計方法與格式，通常都非常龐大、複雜。如果發現整個程式中，有許多片段，都是同樣的語法及內容，就會希望同樣片段只需寫一次程式就夠了，後續若再需要時，只需取得該片段程式的起始執行位址，填在各個一再使用之處，就能完成設計。

這麼簡單的執行邏輯，是難不倒 Forth 熟手的，能解決的先決條件，就是得非常熟悉 Forth 系統內部的結構，就能輕易實現等效功能。

它在 Forth 系統中的用法，是將特定片段的起始執行位址，放在數據堆疊上，待後續想再用到之處，將此執行位址，推放到回返堆疊上去，就完成了工作。待您執行設計好的指令時，重複使用的片段，都能一一被執行出指定的結果。如果您還不能體會這種意義，那麽，請在最後一個 >exec 前面，多添加一列『 dup >exec 』後，重新編譯再執行 greet 指令試試看，您會見到印出了兩次『 world 』這個字，不是只有一次，這就是效果。

實際上，我在 FigTaiwan 論壇上，捐獻過其他類似的設計，只因我熟知系統結構，設計都很簡短。有人提問，能否在 Win32Forth 系統中實現 in-line assembly 的功能，問題本該由 Forth 老手提供解答服務，我看到等了好幾天沒人理會，就動手寫了一個，公開貼出，這種問題對我而言，不是大問題。但我不喜歡這樣用這種技術，那會把程式搞得非常混亂，難以看懂。所以我也不鼓勵大家這樣寫程式，除非一定需要在某處高速、高效執行程式，否則不必硬要這樣寫程式。當然，我也能設計出能在組合語言中突然轉用高階指令插入程式設計之處的古怪程式，也是沒有必要這樣做。

事實上，這個局部副程式的設計方法，根本沒必要存在。因為想這樣用時，還不如就在指令外面另外設計必須反覆使用的部分，給它一個名稱就可以了。後續指令想怎麼用？想用多少次？就都可以悉聽尊便。這本就是 Forth 的基本精神，何須多此一舉又搞一個 quotation 來煩自己？
:

 
\ quotation.f

: ends> r> ; 
: end   >r ; immediate

\ quotations: 局部副程式

: >EXEC >R ;  \ for Win32Forth only

: [: 
       postpone ahead here -rot 
       ends> postpone exit postpone then postpone literal 
   ; immediate 

: ;]  >r ; immediate  \ or an alias of end 

: greet   [: ."  world" ;] 
          ." hello"
          >exec ;      \ for Win32forth
\          execute ;   \ for other forth

第(70)個範例，總共有 8 個依序發展出來的等效程式。

前四個是單純用來判斷一個數字是否為質數的程式，判斷方法，第(68)個範例中已經說明過了。

第一個是最原始的寫法，不求精簡，也不求快速，容易看懂，放諸四海而皆準，不管系統為多少位元，全都有效。

第二個寫法，把所有的邏輯判斷方式都精簡化了，但全仍照根本原則設計程式，適合用來理解語法的精要。

第三個程式，開始添加了 outrun 指令，強調『逸出定量迴路法』。這是目前所有其他種程式語言所沒有的性能。主要的關鍵，是我把迴路指標設計成放在普通的數據堆疊上，所以，我自創的系統才能這樣做。請注意，傳統的 Forth 系統，把『DO ..... LOOP』迴路的指摽，放在回返堆疊上，也能這樣做，但要熟悉系統的人，才有本領照樣做。隨著系統的不同，放法不一，所以， Forth 也不能寫成通用教材來教育大家。我自創的系統，可以。

第四個程式，發展出了對指定範圍內所有數字進行測試的程式。 32 位元情況時之測試記錄，留在程式下方，供作參考。

第五到第八個程式，名稱都改為 PRange 了，但是，他們都是後來適用於 64 位元環境的大型測試用程式。主要運作的系統，是 Lina64 。我不喜歡別種不透明系統的環境，所以只採用完全透明的 Lina64 。

不要小看這一系列的發展進程，我是下了功夫仔細研究過許久的人，知道其中的意義與功效。舉例來說，我在研究純粹以四則運算法設計亂數時，很需要使用一個許多位數的質數來配合應用，所產出的亂數品質才能夠好。若沒有這一系列的工具，你什麼事也別想進行。以前，中山科學研究院請日本技師來為 CDC 大型電腦主機安裝亂數產生程式時，據說，當時是論多少百萬台幣付費的，安裝時，使用硬體裝備灌進記憶體裡燒死，那個關鍵質數放在技師的口袋裡，沒有人在旁的時候，他才 1010 的敲進去。這個質數，若沒有經過大量、長時間的評鑑測試，不敢提出來賣這麼多錢。現在，我的程式，就可以做這件事情。

另外，我曾經在田納西大學的質數網站發現，這所大學花了不少錢，為 54 (非64) 位元以下的質數，建立了不知有多少頁的質數表格，供全世界查用。需要這樣嗎？也許以前紙本書的時代確實需要，不要怪這所名校浪費金錢，畢竟，並不是人人都能像我這樣用 Forth 自創系統，自寫程式，隨時可以得到這方面想要的東西。就算我現在把程式都在網上公開貼出了許久，會用 Forth 的人越來越少，我又不想再搞推廣，能用的，還有幾人？

為了驗證第(68)個範例中所說的所有質數都能寫成 6i+1 或 6i-1 (或 6i+5 )的格式，我也下了功夫，以程式核實，讓系統從頭跑到底，把 64 位元以下所有的質數查證一遍，自動告訴我，有沒有那一個不能寫成這兩種格式的質數？這就形同使用電腦進行全程自動執行的『數學歸納法』驗證數學理論，結果，真是如那個希臘數學家所言：『確實是都行』。留下的這些範例中，最後一個範例 PRange4.f ，就是在做這件事情。

20240927 我在 Twitter 網站上看到這個貼文：100,000,001 can be divided by 17. 我立刻載入因數分解出質數的程式，得知 1,0000,0001 ＝ 17 * 5882353 。別說一億零一了，一千兆零一，我也能立刻得到答案： 1000,0000,0000,0001 = 7 * 11 * 13 * 211 * 241 * 2161 * 9091 。學而時習之，不亦樂乎？

探討質數的理論可以延伸出無數的問體，以後，我還可能遭遇到與質數理論相關的更多問題，這些東西也都還有用處，實際發展出來的程式就代表了我研究問題時的精神。
:

 
\ prime1.f

4 integers D N L I

: test
basic
10 print " Input an integer number N = "
20 inputi N
30 let D = N mod 2
40 if D = 0 then 220
50 let D = N mod 3
60 if D = 0 then 220
70 let L = sqrt ( N ) + 1
80 for I = 6 to L step 6
90 let D = N mod ( I - 1 )
100 if D = 0 then 210
110 let D = N mod ( I + 1 )
120 if D = 0 then 210
130 next I
140 print N ; " is a prime number."
150 goto 300
210 run 2drop
220 print N ; " is not a prime number."
300 end ;

\ prime2.f

3 integers N I L

: outrun 2drop ;

: test
PAGE
BASIC
10 print " Input an integer number N = "
20 inputi N
30 if ( N mod 2 = 0 ) OR ( N mod 3 = 0 ) then 220
40 let L = sqrt ( N ) + 1
50 for I = 6 to L step 6
60 if ( N mod ( I - 1 ) = 0 ) OR ( N mod ( I + 1 ) = 0 ) then 210
70 next I
110 print " Yes, " ; N ; " is a prime number."
120 goto 300
210 run outrun
220 print " No, " ; N ; " is not a prime number."
300 end ;

\S

test

Input an integer number N : ? 600
No, 600 is not a prime number. ok

test

Input an integer number N : ? 1999
Yes, 1999 is a prime number. ok

\ prime3.f

6 integers N I J L C D

: outrun 2drop ;

: test
PAGE
BASIC

10 print " Input start and stop numbers C and D = "
20 inputi C , D

30 for J = C to D

40 let N = J
50 if ( N mod 2 = 0 ) OR ( N mod 3 = 0 ) then 220
60 let L = sqrt ( N ) + 1
70 for I = 6 to L step 6
80 if ( N mod ( I - 1 ) = 0 ) OR ( N mod ( I + 1 ) = 0 ) then 210
90 next I

110 print N
120 goto 220

210 run outrun
220 next J

300 end ;

TEST

\ prime4.f

6 integers N I J L C D

: outrun 2drop ;

: test
PAGE
BASIC

10 print " Input start and stop numbers C and D = "
20 inputi C , D

30 for J = C to D

40 let N = J
50 if N mod 2 = 0 then 220
60 if N mod 3 = 0 then 220

70 let L = sqrt ( N ) + 1
80 for I = 6 to L step 6
90  if N mod ( I - 1 ) = 0 then 210
100 if N mod ( I + 1 ) = 0 then 210
110 next I

120 print N
130 goto 220

210 run outrun
220 next J

300 end ;

TEST

\S
testing range must be <= 2147483647 ( get it from -1 1 rshift . )  
Input start and stop numbers C and D = ? 1999999000 2000000000
      1999999003 
      1999999013 
      1999999049 
      1999999061 
      1999999081 
      1999999087 
      1999999093 
      1999999097 
      1999999117 
      1999999121 
      1999999151 
      1999999171 
      1999999207 
      1999999219 
      1999999271 
      1999999321 
      1999999373 
      1999999423 
      1999999439 
      1999999499 
      1999999553 
      1999999559 
      1999999571 
      1999999609 
      1999999613 
      1999999621 
      1999999643 
      1999999649 
      1999999657 
      1999999747 
      1999999763 
      1999999777 
      1999999811 
      1999999817 
      1999999829 
      1999999853 
      1999999861 
      1999999871 
      1999999873 
      1999999913 
      1999999927 
      1999999943 
      1999999973 

\ prange1.f use \lina64\nfp2\f5102

6 integers N I J L A B

: outrun 2drop ;

: test
BASIC

10 print " Input start and stop numbers A and B = "
20 inputi A , B

30 for J = A to B

40 let N = J
50 if ( N mod 2 = 0 ) OR ( N mod 3 = 0 ) then 220
60 let L = sqrt ( N ) + 1
70 for I = 6 to L step 6
80 if ( N mod ( I - 1 ) = 0 ) OR ( N mod ( I + 1 ) = 0 ) then 210
90 next I

110 print N
120 goto 220

210 run outrun
220 next J

300 end ;

TEST

\S
Input start and stop numbers A and B = 
? 1000000000001000 1000000000002000

    1000000000001003 
    1000000000001027 
    1000000000001063 
    1000000000001089 
    1000000000001117 
    1000000000001209 
    1000000000001269 
    1000000000001293 
    1000000000001347 
    1000000000001371 
    1000000000001413 
    1000000000001491 
    1000000000001503 
    1000000000001551 
    1000000000001617 
    1000000000001623 
    1000000000001669 
    1000000000001741 
    1000000000001749 
    1000000000001819 
    1000000000001839 
    1000000000001867 
    1000000000001897  OK

\ prange2.f use \lina64\nfp2\f5102

6 integers N I J L A B

: outrun 2drop ;

: test
BASIC

10 print " Input start and stop numbers A and B = "
20 inputi A , B

30 for J = A to B

40 let N = J
50 if N mod 2 = 0 then 220
60 if N mod 3 = 0 then 220

70 let L = sqrt ( N ) + 1
80 for I = 6 to L step 6
90  if N mod ( I - 1 ) = 0 then 210
100 if N mod ( I + 1 ) = 0 then 210
110 next I

120 print N
130 goto 220

210 run outrun
220 next J

300 end ;

TEST

\S
Input start and stop numbers A and B = 
? 1000000000001000 1000000000002000

    1000000000001003 
    1000000000001027 
    1000000000001063 
    1000000000001089 
    1000000000001117 
    1000000000001209 
    1000000000001269 
    1000000000001293 
    1000000000001347 
    1000000000001371 
    1000000000001413 
    1000000000001491 
    1000000000001503 
    1000000000001551 
    1000000000001617 
    1000000000001623 
    1000000000001669 
    1000000000001741 
    1000000000001749 
    1000000000001819 
    1000000000001839 
    1000000000001867 
    1000000000001897  OK

\ prange3.f use \lina64\nfp2\f5102

7 integers N I J L A B C

: outrun 2drop ;

: test
BASIC

10 print " Input start and stop numbers A and B = "
20 inputi A , B

30 LET C = 0

40 for J = A to B

50 let N = J
60 if N mod 2 = 0 then 220
70 if N mod 3 = 0 then 220

80 let L = sqrt ( N ) + 1
90 for I = 6 to L step 6
100 if N mod ( I - 1 ) = 0 then 210
110 if N mod ( I + 1 ) = 0 then 210
120 next I

130 LET C = C + 1
140 print N
150 goto 220

210 run outrun
220 next J
230 PRINT " total = " ; C

300 end ;

TEST

\S
Input start and stop numbers A and B = 
? 1000000000001000 1000000000002000

    1000000000001003 
    1000000000001027 
    1000000000001063 
    1000000000001089 
    1000000000001117 
    1000000000001209 
    1000000000001269 
    1000000000001293 
    1000000000001347 
    1000000000001371 
    1000000000001413 
    1000000000001491 
    1000000000001503 
    1000000000001551 
    1000000000001617 
    1000000000001623 
    1000000000001669 
    1000000000001741 
    1000000000001749 
    1000000000001819

2024年9月21日星期六

一百個例題 (61 ~ 65)

Ching-Tang Tseng

Hamilton, New Zealand

22 September 2024

第(61)個範例，展示一套為檔案資料加密的方法。

一百個例題包羅萬象，只因為 Forth 有能力做出所有電腦程式語言能夠辦到的事情，我才能夠同樣地寫出各種千奇百怪的範例程式。

程式名稱很簡單，但是程式內容涉及之事比較廣泛，而且，一般初期發展而成的 Forth 系統，都還不具有所有現成可以直接使用的指令來完成此事，所以要寫出這個範例。

加密的方法有很多種，我只找一個最簡單的方法舉例。其他的任何方法，也都必須遵循同樣的處理程序，才能達到將檔案資料加密的目的。所以，仔細看完程式處理步驟是必要的，也才能體會得出將檔案資料加密後的效果、意義與目的。

這一套簡單的加密方法，是單只採用把每一位元組的資料依次取來，與依序取得密碼之單個相應位元組的資料，進行一是一非則得真的XOR邏輯運算，所得結果，重新放回檔案內的同一位置，短的密碼要被循環著使用，完成整個操作後，將檔案回存，就達到加密的目的了。

為了驗證整個效果，我設計了展示指令： test 。執行後，您可以看到檔案原始的內容，再印出加密後顯示的效果，然後再根據同樣的密碼，反解密回去，重新顯示檔案內容，整個效果就全部顯現了。

這麼簡單的原理，想寫成程式，就沒那麼簡單，整個程序則絕對的完整。如果有一天，您學會了另一套加密的方法，想試一試本領，有此範例程式，就會很容易改寫。加密關鍵只在做 XOR 運算，您學會新的加密運算時，也只需要在此處代換掉 XOR 所做之事，就夠了。

一般 Forth 系統內所缺乏的現成指令，從資料結構的宣告開始，就缺乏了，加密、解密所需的密碼，是被儲存在一個必須留存住長度的字串資料結構內的，這個結構必須自己設計，我給的，也只是個典範，它也可以被設計成是一個非常長的檔案，具有長篇大論的文數字，讓人打字都打不出來。我的設計，只採用最長為 80 個字元的方式儲存密碼。

密碼輸入的對談程式: EnterPassword 也是一個必須自己設計的指令，我的展示技巧，在使用一個憑 CR 跳行(13)，作為最終限定用字的方法，如此一來，密碼就可以是一句話，中間帶有空格也無所謂的設計格式。

存取檔案的對談程式： EnterFileName 也是一個必須自己設計的指令，通常，檔案的整個名稱，不得帶有空格，因此，使用一個 BL 空格(32)，作為輸入時的最終限定用字。事實上，幾乎現行的 Forth 系統，無論是商售或免費的公益系統，都沒有像我系統中所提供的檔案存取現成指令，此事，前曾討論過，有了他們，這裡才能自由使用。我也已經討論過這一套檔案存取的機制，其中就有 FileType 指令可以直接使用，非常方便。

操作執行 test 前，您得先想好自己的密碼，例如：我就用『i like forth』一句話當密碼。也得先想好將被測試的純文字檔案的全名，例如：假設有一個在資料夾中不太重要的 lnx.txt 純文字檔案可供測試。這樣執行 test 後，您就可以看到所有的結果，操作不慎，毀壞了這個 lnx.txt 檔案，也比較無所謂。也許有一天，您能夠用上它，只是要記得，傳送的雙方都必須要先行知道密碼才行，如果密碼與資料都在一起傳送，那加密就是多餘的了。想搞清楚這套用法的精神，自己先用腦筋想一想整套運作的合理協定應該為何？

這套方法雖簡單，想破解密碼也並不容易，誰能那麼輕易的得知隨便給出的密碼？當然並不容易。設定過密碼的電腦，沒密碼就不好開啟。只有你我才知道的密碼，別人當然也就沒那麼容易解碼，其理在此。
:

 
\ (61)Encrypte.f

\ [---head---][lm|ln|---string---]
: STRING ( lm -- )
  CREATE 1 MAX MAXSTRING MIN DUP C, 0 C, ALLOT   ( lm -- )
  DOES> 1+ COUNT ;                               ( -- addr ln )

: S! ( addr1 cnt1 addr2 cnt2 -- )
  DROP DUP 2 - C@
  ROT MIN >R DUP 1- R@ SWAP C! R> CMOVE ;

80 STRING FILEPASSWORD

\ password build here
\ Usage: S" CTT20160310" FILEPASSWORD S!
: EnterPassword ( -- )
  CR ." Enter encrypte password : "
  CR QUERY 13 WORD COUNT FILEPASSWORD S! ;

: EnterFileName ( -- )
  CR ." Enter file name : "
  CR QUERY BL WORD COUNT GET-FILE
;

: FILEENCRYPTE ( -- )
  FLEN 0
  DO
  I FADR + C@
  I    FILEPASSWORD NIP MOD    FILEPASSWORD DROP + C@
  XOR
  I FADR + C!
  LOOP ;

: test ( -- )
  EnterPassword
  EnterFileName
  CR CR FILETYPE
  FILEENCRYPTE
  CR CR FILETYPE
  FILEENCRYPTE
  CR CR FILETYPE
;

第(62)個範例，是一個獨特的整數函數設計。

我一直在研究數學運算系統的學問，搞了四十幾年，心中難免就積聚了許多類比的數學觀念，每當在常用的實數數學體系中，實現各種常用函數的設計時，我就會問自己，別的數學體系中有沒有？或需不需要也有等效的函數？自問的效果，絕對是正面的，身體力行，再去實踐之後，我就能發現更為深入的道理。我也隨時提醒自己，這些東西都很有用處，不要浪費掉這些辛苦研究出來的成果。這個整數體系內以 2 為底的對數函數，就是這麼來的。

函數都有反函數，設計出來的整數函數，當然也就可以有整數的反函數， Log2 通常可以寫成 logb ， ALog2 則為 alogb ，但我在後來的設計中，不再使用這樣的名稱，改為 nlb 的用名了，以強化表示，它是只適用於整數體系環境內的對數函數。相應的其他函數，也因為同樣的類推而產生出來，例如： nlg 為以 10 為底的整數對數函數。我不擬談得太多，範例中也不打算把所有的相關函數設計全都包括進來，那樣做，就會把範例程式搞得非常凌亂，失去了中心要義。

我熟悉數學函數的設計方法，往上、往下延伸研究時，會發現這些東西，都有其連貫性。也就是說，往上，我在自行設計出實數的整套對數函數群時，所有的成果，也都是根基於這個以 2 為底的整數對數函數來設計的。往下，我拜讀過不少硬體數學運算處理器的論文，發現那些能提供對數函數的硬體元件，根本也是基於先有以 2 為底的對數單元，才組合出後來常用的以自然對數 e 為底與以 10 為底所有的其他對數函數。更甚的是，請注意我在 Log2 中的設計，它也是只靠一個向右移位的移位器(rshift)就能設計出來的程式。

研究這樣的軟體程式，是想自行設計製造硬體元件的基礎，但是，中國人還無能設計好用的數學運算處理器(coprocessors)。我們的系統環境從 32 位元進入 64 位元後，系統硬體內的這個元件，並沒有因為位元基礎量翻倍而跟著翻倍，也就是， 32 位元系統中使用 80 位元的 coprocessor ， 64 位元環境中，沒有全面改裝 160 位元的 coprocessor ，而仍是 80 位元的結構。我看過 iForth 系統設計者貼文說過，有這種 160 位元的 coprocessors ，加裝一片要價 250 歐元。中國人將來想不想搞出自己的 coprocessors？沒有類似這個範例程式的研究、設計、與應用，那就難成。

整數對數函數的實際應用，也不僅只是上列的討論而已，我在設計系統時，免不了要為最終成果修邊剪枝，讓函數的最終值能適用於真實環境，就會產生該如何處理最後一位不準位數的問題，怎麼辦？好的系統就有所謂的捨位與截項(round off and truncate)技術應用。請注意！這些技術不單只是講演算式子的截項方法而已，還包括了該如何四捨五入與棄除位數的數字處理。做這些美化成果的工作時，不是見到數字就直接處理這麼簡單，要顧及位數。那麼，取得位數的最快方法，就是利用這個範例中的整數對數函數，我設計的系統中，幾乎所有的常用函數提供最終結果前，都用到這樣的設計。其他的應用，就不多說了。

關於我設計之系統中的各種常用數學函數，在我的個人網頁一篇文章中有具體的圖形展示，文章是 20190216 貼出的：『測試函數』，那裡也談到上述細節中的問題，顯示了我如何好好地把這些研究設計發揮出痛快效果的歷史記錄。
:

 
\ (62)Log2.f

\ Lb(n) = Log2(n) , Lb(8)=3  means 8=2^3
\ ALb(n)= ALog2(n) = 2^n , ALb(3)=8 means 2^3=8

: Log2     ( n -- Log2[n] )
  -1 SWAP
  BEGIN
     ?DUP
  WHILE
     SWAP 1+ SWAP
     1 RSHIFT
  REPEAT
;

: ALog2    ( n -- 2^n )
  1 SWAP 0
  ?DO 2 * LOOP ;

第(63)個範例是一個把字串內容反向印出來的優質程式。

字串處理在 Forth 領域，算是比較複雜一點、屬於二級難度的問題。我在長期使用 Forth 的歷史中，卻經常會碰到有此需求，必須把字串反向印出來。在第(25)個範例中，就曾用到這種技術，但是在那個範例中，我仍然設計了將字串反轉過來的操作指令叫作 turn 。現在這個範例，就能解決不反轉字串，直接反向印出字串的結果，所以我稱它是優質程式，特予收集。

程式的原作者是個住在加州附近的美國朋友，國際論壇上的署名是 Hugh ，網上發言非常彪悍，長期力戰群雄，毫不退縮，大家都對他不好，也很感冒，他卻對我很好，每次寫信給我時，都恭恭敬敬地問問題，包括問中國人的文化。若不論性格， Hugh 有許多可取之處，他在公開場合，大方貢獻出不少漂亮的設計，尤其是字串處理方面的成就，非常優秀。我在個人網頁的貼文中，凡是用到他的設計，必指名道謝，大概也是因為這樣，他才非常尊重我。

Hugh 的專長，本是為數位機械設計計算尺(slide ruler)式的顯示尺，後來，全美工作母機的生產事業崩潰後，他就沒有再做此工作，改為替別人寫 Forth metacompiler 的職業，那時， Forth 最著名的期刊：四度空間(FORTH DIMENSIONS)出版到第 19 年最後一冊，在最後一期內，還可以看得到 Hugh 貢獻出來的文章。所以他也算是 Forth 界的老手了，我一直以禮相待，尊重他。他後來的成就，是專注於具有字串堆疊的 Forth 系統設計，他擁有自己獨創的 novice Forth 系統，把所有的基本指令都設計成專搞字串，所以指令都是以 $ 符號開頭，如： $DUP , $DROP , $SWAP ..... 之類的用法。許多專家對 List 式的資料結構操控能力，沒他熟悉。

他在 XP 時代送過我他所創作之系統的源程式，我則回敬他用組合語言寫成之叫用 coprocessor 完成工作的浮點系統源程式，並指導他利用免費的 MVP Forth 系統，反組譯出指令內容的使用技巧。後來，他在國際論壇上，經常憑此方法與群雄論劍，直到 64 位元普及之後，那套方法解不出 64 位元的碼了，他才停手。

後來的聯繫，都是他為了請教數學問題而來信，由於他只上過高中，我們能談的事情就無法深入。超出高中程度的問題，則是他代替他大學畢業的哥哥向我請教的問題，我都好好的答覆，等到大家都退休了，來信也就少了。

我收集他所貢獻的程式，當然應該清楚交代我跟他交往的事情，表示尊重。以後想再對他有所要求，他才會誠懇對待，我確信，我若開口向他再要 novice Forth 最新版的源程式，他一定會給我，但我若無法深入參與共同研發，最好別去打擾，留著他當一輩子的朋友，就像我現在與大家交往的情誼一樣。

我在研究首個 64 位元 Forth 系統時，使用純用 C 寫成之模擬式的 Sod64 Forth 系統，因為那時還沒有公益性的 64 位元 Forth 系統可用，只有這一套。但因它係以 C 在 32 位元環境中模擬 64 位元設計而成，全世界都沒有人對它有興趣，我只好自己精研。我在其中加設計浮點系統時，才發現系統將待印出的數字都擺在 PAD 以下，向 HERE 的方向成長，反正沒有 256 個數字的東西要印出來，碰不到 HERE 點，不會破壞系統，這樣擺就無所謂。可是，對我的發展影響很大，每個數字都得反過來印出，不先反轉就不能直接使用 type 指令，還真是不好處裡，為了這事，我苦惱了很久，如果那時就有 Hugh 公開捐贈的這個 backtype 程式，問題就很容易解決了，但這個程式是後來才取得的。別篇網文中也曾提及這個 backtype 指令的應用，這個範例則詳述了他的來源。

我在設計無限位數的除法設計時，也碰到除得的商數，是一個一個依序擺進字串的，想印出結果時，就得反過來印，就得使用這個 Hugh 的貢獻，才會方便，後面的範例中也許會再論及這套技術。請注意，他沒有把字串先翻轉過來後才印，是就地處理，設計方式算得上是比較高明。
:

 
\ (63)BackType.f

\ 20160314 Hugh Aquilar posted on c.l.f.

: backtype ( adr cnt -- )
    over +  1- 
    2dup u<
    if
        do  I c@ emit  -1 +loop 
    else 
        2drop
    then ;

: $>pad ( adr cnt -- )
  dup pad c!
  pad 1+ swap cmove ;

s" this is a test." $>pad
cr cr pad count type
cr cr pad count backtype

第(64)範例是一個 Forth 的特殊應用程式，讓電腦操控喇叭，演出單音音樂。

以前的個人電腦，可以自行將控制值放進指定記憶體內，直接驅動喇叭，演出音樂。後來的個人電腦，只能透過作業系統提供的功能程式來驅動喇叭。前者有它的方便性，後者有它的多樣性。我們在這個範例中，只利用一個很單純的功能，只提供頻率與震盪時間長度來發出聲音，效果當然不會很好，但是能告訴大家，軟體系統是如何辦到這些控制硬體之事的，也就是說能夠告訴您，這種程式該如何設計？

完美可叫用的音效功能程式，還可以包括波譜的選擇，便能模擬出各種樂器的聲音，甚至於可以模擬出名聲樂家的話語，還能調聲調高低，令其唱歌，清唱也行，唱歌劇也行。現在這個範例，只是最基礎、最初步的發聲方法，但，再複雜的功效，程式也得照同樣的方法來設計，才能達到目的。我也沒有增加能夠配合作業系統同時演出兩個檔案程式的設計，來辦到能像交響樂那樣演出的效果。這樣的設計，我也能夠辦得到，只是需要花費很多時間，卻與數學計算無關。

發聲的整個關鍵在給音長與頻率兩個參數，在 Win32Forth 系統中有內建的 Beep 指令可以叫用，在 Wina32 Forth 系統中，以及在 FasmForth 系統中，也都可以先連結起作業系統中的動態連結程式(.dll)後，再叫用其中指定名稱的指令，來達到目的。另外兩個附屬的 (64-1) , (64-2) 範例，只展示它們的叫用方法，編寫音調程式的部份則不列出。所有 Forth 系統叫用作業系統功能程式的設計原理，都完全一樣，搞通這個範例，就能放諸四海而皆準。

實際上，好的音樂程式，沒這麼單純，但要設計得好，就得在音樂的樂理上有點素養，我有比較完整的設計，並將成果貼文於我的個人網頁，而且全用英文寫成的文章，文內附貼我自己使用那一套卡拉 OK 式的程式，唱出『茉莉花』中文歌的短片錄影。

該文有許多音樂方面的特色，單憑兩個參數，能調整出任意的音調與速度，曲譜是我自創的文字譜，曲調也可以隨我自己的聲調來調整，速度板拍的宣告，就用音樂專用的術語來表示，整套音樂功能盡量類似於五線譜的樂譜表達方式，文字譜則選用『首調』唱法，而不用『固定調』唱法來寫文字譜。樂理，我只點到為止，我略懂。

歐洲及俄國人擅長把音樂搞成全國普及的程度，所以他們參訪該網文的數目很多，看完後，也在國際論壇網頁上推介。因為他們知道，那才是比較接近正規音樂的程式設計方式。程式的內容比較大，不宜當範例教材，大家若有興趣聽我真人實唱，請自行前往參考。該篇文章可以永遠留在網上，也就是本人唱歌的聲音，可以遍傳全世界許多年，這樣做比較有意義。
:

 
\ (64)tone.f

\ 20160402
\ 192 = 2*3*4*8

192 3 * VALUE Duration
\ n1:Duration   n2:Tone Frequency
: sing ( n1 n2 -- ) CALL Beep drop ;
: silence ( n1 -- ) MS ;
: | ; immediate

\ Get duration ( -- n1 )
: 4T   Duration 4 * ;
: 3T   Duration 3 * ;
: 2T   Duration 2 * ;
: 1T   Duration     ;
: T/2  Duration 2 / ;
: T/3  Duration 3 / ;
: T/4  Duration 4 / ;
: T/8  Duration 8 / ;
: T/16 Duration 16 / ;

\ Get tone frequency ( -- n2 )
: MDO 523 ;
: MRE 587 ;
: MMI 659 ;
: MFA 698 ;
: MSO 784 ;
: MLA 880 ;
: MSI 988 ;

\ Sing it ( n1 -- )
: Dl MDo sing ;   : LDo MDo 2 / sing ;   : HDo MDo 2 * sing ;
: Re MRe sing ;   : LRe MRe 2 / sing ;   : HRe MRe 2 * sing ;
: Mi MMi sing ;   : LMi MMi 2 / sing ;   : HMi MMi 2 * sing ;
: Fa MFa sing ;   : LFa MFa 2 / sing ;   : HFa MFa 2 * sing ;
: So MSo sing ;   : LSo MSo 2 / sing ;   : HSo MSo 2 * sing ;
: La MLa sing ;   : LLa MLa 2 / sing ;   : HLa MLa 2 * sing ;
: Si MSi sing ;   : LSi MSi 2 / sing ;   : HSi MSi 2 * sing ;

: Waltzing-Matilda
   1T SO   T/2 SO  T/2 SO  1T SO   1T MI
   1T HDO  T/2 HDO T/2 SI  1T LA   1T SO
   1T SO   T/2 SO  T/2 SO  1T LA   1T SO
   1T SO   T/2 FA  T/2 MI  1T RE
   T/2 Dl T/2 RE  1T MI   1T MI   1T RE  1T RE
   T/2 Dl T/2 RE  T/2 MI T/2 Dl  T/2 LLA T/2 LSI  1T Dl
   1T LSO  T/2 Dl T/2 MI  1T SO   T/2 FA T/2 MI
   1T RE   T/2 RE  T/2 RE  1T Dl  ;
: Jasmine
  1T MI T/2 MI T/2 SO T/2 LA T/2 HDO T/2 HDO T/2 LA 1T SO T/2 SO T/2 LA 2T SO
  1T MI T/2 MI T/2 SO T/2 LA T/2 HDO T/2 HDO T/2 LA 1T SO T/2 SO T/2 LA 2T SO
  1T SO 1T SO 1T SO T/2 MI T/2 SO 1T LA 1T LA 2T SO
  1T MI T/2 RE T/2 MI 1T SO T/2 MI T/2 RE 1T Dl T/2 Dl T/2 RE 2T Dl
  T/2 MI T/2 RE T/2 Dl T/2 MI 1T T/2 + RE T/2 MI 1T SO T/2 LA T/2 HDO 2T SO
  1T RE T/2 MI T/2 SO T/2 RE T/2 MI T/2 Dl T/2 LLA 2T LSO
  1T LLA 1T Dl 1T T/2 + RE T/2 MI T/2 Dl T/2 RE T/2 Dl T/2 LLA 2T LSO ;
: I-am-a-little-bird
  1T DL 1T DL 1T DL 1T T/2 + MI T/2 RE 1T DL
  1T MI 1T MI 1T MI 1T T/2 + SO T/2 FA 1T MI
  1T SO 1T FA 1T MI 3T RE
  2T RE T/2 DL T/2 SI 1T DL 1T RE 1T MI
  2T FA T/2 MI T/2 RE 1T MI 1T FA 1T SO
  T/2 SO T/2 FA 1T MI 1T RE 2T DL ;
: Pocalicalina
  1T SO T/2 HMI T/2 HMI T/2 HRE T/2 HRE 2T HDO 3T HMI
;

: TEST 1 0
  DO
  Waltzing-Matilda 2T silence
  Pocalicalina 2T silence
  Jasmine 2T silence
  I-am-a-little-bird 2T silence
  LOOP ;

\s
The following table shows the relationship of notes and their frequencies
in one octave.
_______________________________________________________
C       D       E       F       G       A       B
Dl      Re      Mi      Fa      So      La      Si
261.63  293.66  329.63  349.23  392.00  440.00  493.88
523.25  587.33  659.26  698.46  783.99  880.00  987.77
-------------------------------------------------------

By doubling or halving the frequence, the coinciding note values can be
estimated for the preceding and following octaves.

\ (64-1)play music in wina forth

: Z ( sc -- adr)
  0 , DROP ;                                    
                                                  
: make-constant ( n adr -- )                                                
  BODY> >R    
  R@ >DFA  !   
  'BL >CFA @   
  R> >CFA ! ;           
                                                   
: LOAD-DLL: ( sc -- u ) 
  CREATE $, DROP 
  DOES>  DUP >R $@  LOAD-DLL           
         DUP R> make-constant ;                                     
                                               
: DLL-ADDRESS:  ( sc xt -- adr ) 
  CREATE , $, DROP     
  DOES>  DUP >R   CELL+ $@  
         R@ @ EXECUTE DLL-ADDRESS   
         DUP R> make-constant ;          

"kernel32.dll" LOAD-DLL: K32 
"Beep" 'k32 DLL-ADDRESS: BBB
: sing ( time frequency -- )
  BBB call drop ;

\ (64-2)FasmForth64Tone.F
\ REQUIRE API_2: lib\WAPI.4
\ REQUIRE API_2: WAPI.4

: DLL_CREATE CREATE  PARSE-NAME FILE-BUFF ASCII-Z DLL_S
  DUP 0= ABORT" api unavailable" , ;

: API_0: DLL_CREATE DOES> @ API_0 ;
: API_1: DLL_CREATE DOES> @ API_1 NIP ;
: API_2: DLL_CREATE DOES> @ API_2 NIP NIP ;
: API_3: DLL_CREATE DOES> @ API_3 NIP NIP NIP ;
: API_4: DLL_CREATE DOES> @ API_4 NIP NIP NIP NIP ;
: API_5: DLL_CREATE DOES> @ API_5 NIP NIP NIP NIP NIP ;

\ CR .( LIB LOAD=)  S" KERNEL32.DLL" DROP DLL_L DUP H.
\ DUP 0= [IF] .( KERNEL32.DLL LOAD ERROR) ABORT [THEN]

S" KERNEL32.DLL" DROP DLL_L
CONSTANT KERNEL32DLL
 
\ Get API Beep
KERNEL32DLL API_2: Beep Beep
: sing ( duration freqency -- )
  Beep DROP ;

第(65)個範例程式，能像發送摩爾斯(Morse)電碼一樣，產生出發送電報時的聲音。

這個程式的原始碼，可以追溯到 1980 年代。早期 IBM XT 個人電腦時代，在 Poly Forth 系統中就附贈這樣的 Demo 程式。

它的源程式就是附隨在這個範例程式最前面的部份，但已不能執行，因為源程式中使用了 pc! 這個指令，那是對 isolated I/O 記憶體進行直接操控時所使用的指令，在視窗軟體環境中，對執行軟體進行了所謂保護模式狀態下才能工作的限制，就不能再使用這種指令了。

為了驗證 W10 還能工作於此範例，我特地仔細測試了一下，才發現，當初，我是為了測試 Wina32 系統而使用此例發展程式，不是為了 Win32Forth 系統而改寫的範例。但沒有關係，只需在編寫 Morse 電碼程式的前面，插入一段擷取來自 Wina32 的設計，重新載入程式就能執行了。令電腦發出聲音的方法與前一(64)範例相同，也是叫用系統已有的 Beep 功能程式而實現的，所需的 dit , dah 兩個聲音，由頻率決定。

從最古老的設計，推進到最現代的環境，我一直沒有中斷過使用 Forth ，所以熟知問題所在。現在，回頭去介紹古老 8086 CPU 的工作機制，是沒有必要的事情。若要仔細談起來，這種搞電腦對外控制的技術，還有像 6502 CPU 那樣的 memory map I/O 的學問要談，會越扯越多。這些東西我都熟，在以 Forth 為主之軟硬體對應關係上的發展技術，因為我追用時間很長、又未曾中斷過使用 Forth 的關係，還算熟悉，能夠快速的改出程式。

能適用於 W10 作業系統的這個範例，整個程式修正如下：
:

 
: tone ( time frequency -- )
  CALL Beep drop ;

\ morse demonstation begins here 
880 constant freq    \ 440 --> 880,         value --> constant
 45 constant adit    \ 1 dit will be 45 ms, value --> constant 

: ns 4 * ;           \ convert adit to duration value, without nano second meaning

: dit_dur     adit ns ; 
: dah_dur     adit 3 * ns ; 
: wordgap     adit 5 * ns ;
: off_dur     adit 2/ ns ; 
: lettergap   dah_dur ; 

: sound   ( duration -- )   freq tone ;
: silence ( duration -- )      0 tone ;

: MORSE-EMIT  ( char -- )               \ send morse and echo to console 
        dup  bl =    \ bl = blank = 32
        if 
             wordgap   silence drop 
        else 
             PAD C!                     \ write char to buffer 
             PAD 1 EVALUATE             \ evaluate 1 character 
             lettergap silence          \ pause for correct sounding morse code
        then ; 

: BOUNDS ( ADDR LEN -- Last Init )
  OVER + SWAP ;

: TRANSMIT ( ADDR LEN -- ) 
           CR                          \ newline, 
           BOUNDS                      \ convert loop indices to address ranges 
           DO 
              I C@ dup emit            \ dup and send char to console 
              MORSE-EMIT               \ send the morse code
           LOOP ;

: .   ( -- ) dit_dur sound  off_dur silence ;
: -   ( -- ) dah_dur sound  off_dur silence ;


\ define morse letters as Forth words. They transmit when executed 

: A  . -  ;     : B  - . . . ;   : C  - . - . ;    : D  - . . ; 
: E  . ;        : F  . . - . ;   : G  - - . ;      : H  . . . . ; 
: I  . . ;      : J  . - - - ;   : K  . - . ;      : L  . - . . ; 
: M - - ;       : N  - . ;       : O  - - - ;      : P  . - - . ; 
: Q  - - . - ;  : R  . - . ;     : S  . . . ;      : T  - ; 
: U  . . - ;    : V  . . . - ;   : W  . - - ;      : X  - . . - ; 
: Y  - . - - ;  : Z  - - . . ; 

: 0  - - - - - ;     : 1  . - - - - ; 
: 2  . . - - - ;     : 3  . . . - - ; 
: 4  . . . . - ;     : 5  . . . . . ; 
: 6  - . . . . ;     : 7  - - . . . ; 
: 8  - - - . . ;     : 9  - - - - . ; 

: '  - . . - . ;     : \  . - - - . ;    : !  . - . - . ; 

: ?  . . - - . . ; 
: ,  - - . . - - ; 
: /  . . . - . - ;  ( SK means end of transmission in int'l Morse code) 
: .  . - . - . - ; 

S" CQ CQ CQ DE VE3CFW VE3CFW / " TRANSMIT

2024年9月18日星期三

一百個例題 (56 ~ 60)

Ching-Tang Tseng

Hamilton, New Zealand

19 September 2024

第(56)個範例是一個強行使用 Forth 來模擬出一個益智問題的程式。

問題為：一瓶啤酒 2 元，兩個空瓶可以換一瓶啤酒，四個瓶蓋也可以換一瓶啤酒，請問 n 元可以喝到幾瓶啤酒？

這個問題據說是蘋果公司招考員工時的考試題目，我不想考究。網上也有不少人談論這個問題，看了都只覺得那些網貼是話太多、胡扯的多、強辯無益者更多，所以，大可不必去網搜，我也保證此前沒有以 Forth 程式解出此問題之答案的網貼。

問題取材的來源，是臉書上一位武陵中學高我八屆的莊姓學長所提出的問題，這位學長知道我用電腦設計系統，於是轉貼了這則網上買啤酒的熱門話題，問我有何看法？也問我單用程式來解題的話，程式要怎麼寫？我想了一想，就先以解析法給了他一個答案：就說它是一種物質不滅定律的延伸，姑且就叫作價值不滅定律吧。用10元除以啤酒淨值0.5元/瓶，得20瓶，但奇數塊錢不能這樣解析。至於程式解法，我隔天再試寫給他看，於是就寫出了這則範例程式。程式執行的結果，獲得的答案是 15 瓶，不是 20 瓶。我當這件事是個可貴的經驗，於是將程式設計的結果，留作範例而收集。

不要小看這則問題，也不要小看這個程式，程式是典型的模擬設計，直接全程模擬如何按照條件買啤酒，買到最後，就知道能買得幾瓶了，這裡面沒有解析性的投機演算，直接採用暴力算法(brute force algorithm)模擬實際問題。這是有了電腦以後，人類才可以妥善使用電腦計算能力的解題用法：直接進行模擬。

詳細解法的中文說明，範例程式中已經記錄，此 Forth 程式的特點則在使用 value 定義變常數，而不使用 variable 來宣告出純變數的資料結構。好處是可以少用 fetch(@) 指令，讓程式易於閱讀。 Forth 系統中泛用的可以加存 +to 或 +! 指令，都對解這類問題而設計程式有很大的幫助。注意迴路結束條件，就在模擬到買不下去時，就得終止。我們分析問題、建立模式、設計程式、獲得解答，整個過程，在口語式的問題中並沒有提供，你得有點本事自己想得出來才行。

模擬程式中的運算，主要依賴一個 Forth 中特有的除餘(/MOD)運算指令，它執行兩個輸入數字相除後會在堆疊上留下商數與餘數兩個數字的結果。80 年代，只有 Forth 程式語言採用除法運算後能直接獲得餘數的指令，後來，其他程式語言才跟進使用。

範例程式中使用了 begin ..... until 的迴路結構指令組， until 前面的迴路終止條件，就是再也無法取用空瓶或瓶蓋換出新啤酒來了的狀況。

程式除了模擬這個問題外，還必須能適用於各種狀況：無論啤酒、空瓶、瓶蓋各自的價值如何變化，也無論問題換用多少錢作為起始，程式都要能夠涵蓋所有的狀況才行。
:

\ (56)Beers.f

\ 一瓶啤酒 2 元，兩個空瓶可以換一瓶啤酒，四個瓶蓋可以換一瓶啤酒。
\ 請問 n 元可以喝到幾瓶啤酒？

\ t : 可得啤酒總數
\ r : 循環回收次數
\ a : 暫態計算結果
\ b : 空瓶數
\ c : 瓶蓋數

5 values t r a b c

\ 只對偶數的起始費用有效，奇數元必須先扣除1元
: Init ( n -- )
  0 to r
  2 /mod dup dup to t  to b  to c
  if cr ." 1 dollar redundant." then ;

: Recycle  ( -- )
  0 to a
  b 2 /mod +to a  to b
  c 4 /mod +to a  to c
  a +to t  a +to b  a +to c
;

: main ( n -- )
  Init
  begin 1 +to r Recycle  a 0= until
  cr ." Total = " t .
  cr ." 回收次數 = " r .
  cr ." 剩餘空瓶數 = " b .
  cr ." 剩餘瓶蓋數 = " c .
;

\ 經過分析可得下列結論
\ 啤酒的純脆價值為每瓶0.5元，每個空瓶為1元，每個瓶蓋為0.5元。
\ 因此，若以 n 元購此啤酒，兌換到最後，必定只剩 1 個空瓶，3 個瓶蓋。
\ 以此結論，各樣東西的價值不會無中生有，也不會無中消失，故
\ ( n - 1 - ( 3 * 0.5 ) ) / 0.5 便為可獲得之所有啤酒瓶數的等效金錢。
\ 但此式只對偶數的起始金錢有效，奇數金錢，一開始就必須結餘1元，不列入計算。
\ 例如：花 23元，可獲得 ((23-1)-2.5) / 0.5 = 39 (瓶)

\S
10 main 
Total = 15 
回收次數 = 6 
剩餘空瓶數 = 1 
剩餘瓶蓋數 = 3  ok
23 main 
1 dollar redundant.
Total = 39 
回收次數 = 9 
剩餘空瓶數 = 1 
剩餘瓶蓋數 = 3  ok
57 main 
1 dollar redundant.
Total = 107 
計回收次數 = 13 
剩餘空瓶數 = 1 
剩餘瓶蓋數 = 3  ok
100 main 
Total = 195 
回收次數 = 15 
剩餘空瓶數 = 1 
剩餘瓶蓋數 = 3  ok

第(57)個範例是個解不等式問題的典型範例。

此範例的收集，與我在 20170301 個人網頁貼出的 ”Daily using in ABC Forth” 文章有關，那篇文章用來向全世界推介我設計的系統，所以全用英文。文章內題目的來源，都是從一個叫做 Brilliant.org 的網頁內取得的網貼內容。該文貼出時，一口氣展示了 13 個題目的程式解法，內容已經非常多了，不宜再添加，所以這一題沒有補進該文，我將其收集進百例此處，以供永久留參。

還有一個特別的因素，我才沒有將此範例補貼於上述該文。那篇文章內所有的程式解法，我都採用 BASIC 式的語法設計程式，如果是使用純 Forth 語法設計出來的程式，我就不補進該文，改為自己收集。

從上述的差異說明，大家可以仔細比對，看看用 BASIC 語法或用 Forth 語法設計解數學問題的程式時有什麼不同？這個範例只需用到一個非常簡單的邏輯判斷來取得結果，連一個變數的資料結構都不需要先行宣告出來，就能完成任務了。

可是算式的寫法：

( n^2 - 2 ) ( n^2 - 20 ) < 0

就相當傷腦筋了，我自己若要回查，不看這一列說明，也得花不少時間才能重新從程式倒寫得出這個不等式來。

這一個例題是使用純脆的 Forth 語法寫成的。您能自己改用我給的 BASIC 語法方式另外寫出解這個問題的程式嗎？寫成之後，可以比較一下，就能感覺出它們的不同處到底在那裡了。

我完成系統設計後，長期絞盡腦汁設想過許多問題來考驗系統。沒有網路的時代，只能靠讀書來找問題，有了網路之後，方便多了，但早期可用網頁的資料，還是非常欠缺，需要時間來累積。我在步向這條路途的過程中，也見過許多網頁起起伏伏的興衰，許多網頁早就從網上徹底消失了。目前，雖仍還有不少網頁仍能參考，但總覺得他們都缺乏系統性的規劃，能用材料都屬於點到為止的東西，不是可以猛下功夫去全面追隨，然後重新實現於 ABC Forth 的理想對象。

我接觸過的 Brilliant.org 例題，大該是最平易近人題目的來源。類似 Rosetta code 的網頁例題，則是比較通俗性的題目來源，類似的網頁也有很多個，有偏重于生物學方面的範例網頁，也有偏重于天文學方面的範例網頁，我都大量造訪過，然後取材來磨練自己。題目比較高深難解的網頁，則是專論理論數學方面的網頁，例如 Euler’s 問題網頁，我也曾一口氣下載了幾十個問題一一解出答案，但因這個網頁要加入該組織、登錄個人資料後才能自由出入，我不想惹事，就放棄了，沒去參與。百例的取材大部份是這樣形成的。

事實上，我並不鍾情於以上述資源來設計程式，取得成果，它們都屬於只是為了向大家舉例，才不得不為之的事情。

我個人真正的喜好，是從優質論文中取得的資源。已經寫過不少了，也還有不少題目尚未實現，題目都很大，技術也不容易全面仔細學得通，但值得做。

我個人真正的願望，是找幾項專業工程，精細的實現整套的應用實務。例如：設計出全套的機械工程套裝程式，或是完成整套的數值分析實用庫存程式。

回顧這個範例，令我想到 Forth 的未來，確實是越來越難再走下去。介紹百例能有效果的大前題，是大家應該都已熟知了 Forth 基本指令與程式語法了才行。可是，我自知，這已是苛求了，將來會更難。所以，推介也好、喜好也好、願望也好，都不要期待太高，永久這樣使用 Forth ，是為自己而用它一輩子。
:

\ LessThan.f
\ How many integers n satisfy this inequality?
\ ( n^2 - 2 ) ( n^2 - 20 ) < 0
\ Can you find them all?
\ use FORTH code solve it as following:

: main ( n -- )
  dup dup * dup >r
  2 - r> 20 - *
  0 <
  if   cr .
  else drop
  then ;

: TellMe ( to from -- )
  do I main loop ;

\S
\ Usage:
100 -100 TellMe
-4 
-3
-2 
2 
3 
4  ok

第(58)個範例是用來檢討暫態數據可以被放置在系統內何處的問題。

只有 Forth 程式語言允許使用者能夠看透整個系統，也允許使用者能夠更動系統的內容。當然，後果得自行負責。

不能達到上述目的的 Forth 系統，我是一概不用，因為我必須在 Forth 系統上建立自己的延伸性附屬系統，能夠擁有上述功能者，我才願意使用這樣的 Forth 系統。 AIM65 Forth, Micromotion Forth-79, F83, eForth, FPC-Forth, Win32Forth, ciForth ..... 等等，都是符合上述規格的系統，所以我用。其他的系統，我若願意討論，都只是輕描淡寫地玩一玩就算了，絕不浪費時間去深入探討或使用。

能把 Forth 系統玩通，有兩件事情必須掌握，一個就是指令結構，另一個就是規劃出來的系統記憶體圖譜。目前市面上販售的 Forth 系統與歐洲推行的系統，都不鼓勵使用者這樣子用系統，如： Swift Forth, MVP Forth, iForth, gForth ..... 等都是。這些系統我都用過，但不公開討論，也不鼓勵任何人去用，因為系統不透通，用了就等於被害、被綁死，沒有發展的價值。其他許多上不了檯面的殘障 Forth ，我就不提了，連碰也懶得碰，碰了就是浪費時間。

為什麼先批評一段 Forth 的各種系統再談今天的範例？因為這個範例牽扯到系統的概念，必須先搞懂系統結構才行。

我在長期使用 Forth 之後，知道使用者經常會有暫態數據必須處理的問題。於是，就會問自己，如果不宣告出變數來用時，數據可以擺在那裡？又因 Forth 程式語言強調使用堆疊來運轉系統，堆疊上數據的擺放秩序，嚴重影響執行需求，設計程式時，經常就得把堆疊上的數據進行翻來覆去的處理，被處理的數據就是暫態數據，不找個地方暫時放著，就難以翻來覆去的搞出正常的秩序。還有，每當我見到 Forth 指令或變數的頭部長得這麼龐大、耗掉這麼多記憶體時，我就會想用經濟可靠的記憶體使用方式處理暫態數據，而不想規規矩矩的使用宣告出變數名稱的方法來儲存暫態數據了。

關於儲存暫態數據的技術，有很多可行辦法，方法不是只有這一種，這個範例只是其中比較通用的一種，也是比較傳統的使用方法。我在程式說明中介紹了一個可以存放於三處的使用範例，並且只拿處理堆疊上數據放置秩序當例題，展示現象，都能得到正確的結果。然後，我強調了直接設計出沒有儲存上限的浮動陣列式使用方法解決問題，給了一個通用性的 [PAD] 陣列結構，能被萬用。我後來再設計系統、加添許多固定功能時， [PAD] 被永久固定在我的系統中了，只是名稱可能有點不同，每個陣列單元記憶體的控留量，也會因環境的不同而有點差異，如此而已。

使用本範例中的記憶體位置來儲存暫態數據時，必須注意一項要點，由於儲存位置都是浮動的，因此，只宜在系統沒有編譯新的程式進入系統時使用，系統編譯納入新指令時，這些儲存位置所依賴的 HERE , PAD 指標，都會跟著浮動變化掉。而回返堆疊指標，也必須在編譯完一個指令後歸原，所以程式中的暫態數據，只能臨時放一放，編譯完成前，回返堆疊的指標必須恢復正常。

我設計過許多儲存暫態數據的不同程式，這是個龐大的題目，沒有碰到實際程式時，我們就不討論，碰上一種就談一種。如果您也想自行設計系統，各種技巧就都得深入探討。而且，這類技術也算是保護技術的一種方式，當別人不知道暫態數據被放在何處時，程式就難以追蹤，想要搞反編譯來破解系統時，就會遭遇到重大的阻礙。
:

\ (58)PadHere.f
\ 20151201
\ 利用 PAD, HERE, Return Stack存放變數的技巧
\ 只適用於 4 個數字時使用
: test ( 1 2 3 4 -- 4 3 2 1 )
  here ! pad ! swap >r >r
  here @ pad @ r> r> ;

\ 均適用於整數、實數、複數，故一個單元就取 4 個 cells
\ : [PAD] ( n -- addr )
\   4 * cells pad + ;
\ 比較理想的設計為包括也考慮 bytes/float = 10 的情況
\ 如果記憶體用量根本不是問題時，每一 單元乾脆就取 80 Bytes.字串也能使用

: [PAD] ( n -- addr )
  80 * PAD aligned + ;

: [PAD]N.S       ( from to -- )
  CR 1+ SWAP DO I [PAD] @ . LOOP ;
: [PAD]F.S       ( from to -- )
  CR 1+ SWAP DO I [PAD] F@ FS. LOOP ;
: [PAD]Z.S       ( from to -- )
  CR 1+ SWAP DO I [PAD] Z@ Z. LOOP ;

: ntest ( 1 2 3 4 -- 4 3 2 1 )
  CR .S
  1 [pad] ! 2 [pad] ! 3 [pad] ! 4 [pad] !
  1 4 [PAD]N.S
  1 [pad] @ 2 [pad] @ 3 [pad] @ 4 [pad] @
  CR .S
  ABORT
;

: ftest ( f: 1e 2e 3e 4e -- 4e 3e 2e 1e )
  CR F.S
  1 [pad] f! 2 [pad] f! 3 [pad] f! 4 [pad] f!
  1 4 [PAD]F.S
  1 [pad] f@ 2 [pad] f@ 3 [pad] f@ 4 [pad] f@
  CR F.S
;

\ Usage:
\ 1 2 3 4 ntest
\ 1e 2e 3e 4e ftest

第(59)個範例也是展示 Brilliant.org 網頁問題的程式，採用暴力算法時，介紹另種用法，若另掺入分析技巧，便能將程式擴展應用至非常大數字之問題去了。

問題為： n^n 級數所有方次加起來後的總和之最後一位數到底是多少？

運算後溢量(overflow)的問題，永遠都是個困擾電腦計算技術的問題，尤其是在利用電腦執行暴力算法時更甚，使用於 32 位元環境時，大約是十億左右， 64 位元環境，雖能處理得更大，大到十八位數。實際上，對研究大數問題而言，容量仍是相當有限。這個範例中的程式，就在展示能表現數字位數受限時的影響。我故意不強調我的系統能算大數，才能讓溢量問題凸顯出來。範例程式後面，列有執行結果，小心比對最後計算出來的大數字，就能明白系統在方次多大時，顯示出來的結果是溢量數字了。

問題的主要精神，也是在強調應該分析問題後才解題目， 10 的 10 次方，就是十億了，使用暴力算法硬解此題，顯然是會弄爆可用整數上限的，改換 64 位元來處理，也不過是只能夠提高幾位數而已，我展示了算到 15 次方還不會有問題。這些展示，只是提醒自己，要注意數字處理上限，真正的解題技巧，不能完全只用暴力，採用分析，就能減輕負擔。因為問題只問您 n^n 級數所有方次加起來後的總和之最後一位數到底是多少？

只問最後一位數，我們就可以回頭考慮方次與加法的運算，不管您怎麼算，最後一位數是可以完全納入掌握的數字，於是，我們就能設計自己的程式，單純的只專注於最後一位數，可被留供後續演算，便能解決方次很大時溢量的問題。

demo1 與 demo2 兩套暴力算法，都只能處理到 11 次方。我不以可解決小方次問題時的能力為滿足，所以擴充設計， test 便能讓您使用到非常大的方次，因為每次算完都只留最後一位數，問題便能迎刃而解。

分析後才設計程式的處理要點，都以註解說明的方式寫在程式後面了，本文不再說明。

我在分析問題時，另外想知道最後一位數是不是也有限制規格？於是，又設計了一個 main 指令，測試的結果，顯示這個所謂的最後一位數，從 0 到 9 都有可能。所以，分析問題時，也不必再從這方面找規則來簡化程式了。電腦、系統、思考、樂趣都是不花錢的組合，多寫幾個程式滿足自己的認知，不亦樂乎？
:

\ (59)LastDigit.f

((
https://brilliant.org/practice/level-2-3-operations/

What is the last digit of this sum?

Find the last digit when

1^1 +  2^2 + 3^3 + .... + 9^9 + 10^10

is written out as an integer.
))

6 integers i j k l p s

: demo1   	( n -- )   \ maximum less than 11
[[ j ]] !
basic
10 let k = 0
20 for i = 1 to j
30 print " " ; i ; " ^ " ; i ; " = " , i ^ i
40 let k = k + i ^ i
50 next i
60 print " Sum( " ; j ; " ) = " , k
70 end
;

: demo2   	( n -- )   \ maximum less than 11
[[ j ]] !
basic
10 for k = 1 to j
20 let s = 1

30 for i = 1 to k
40 let s = s * k
50 next i

60 print s
70 next k

80 end
;

: test   ( n -- )               \ maximum can up to very huge
[[ j ]] !                       \ 最高方次為 j (highest power is j)
basic
10 let s = 0

20 for k = 1 to j               \ k 表從 1 開始，做到最高方次 j (k terms)
30 let l = 1

40 for i = 1 to k
50 let l = ( l * k ) mod 10     \ 每自乘一次均只留單個尾數 l (last digit)
60 next i

80 let s = s + l                \ s 為所有尾數之和 (summation of last digits)
90 next k

100 let s = s mod 10            \ 取得尾數總和之最後一位尾數
110 run s .
120 end
;

: main ( -- )
basic
10 for p = 0 to 20             \ 算出最高方次 p(power) 為從 0 到 100 的所有結果
20 print " Last power = " ; p ; " last digit is "
30 run p test
40 next p
50 end ;

\ page main

\ 執行 main 後，結果顯示，0 到 9 都有可能。

\S
 ok
10 demo1 
1 ^ 1 =                1
2 ^ 2 =                4 
3 ^ 3 =               27 
4 ^ 4 =              256 
5 ^ 5 =             3125 
6 ^ 6 =            46656 
7 ^ 7 =           823543 
8 ^ 8 =         16777216 
9 ^ 9 =        387420489 
10 ^ 10 =     1410065408  ----> overflow
Sum( 10 ) =   1815136725  ok

ching@ching-H81M-S2H:~$ cd lina64
ching@ching-H81M-S2H:~/lina64$ ./f5082

AMDX86 ciforth 5.2.1 
include lastdigit.f
I : ISN'T UNIQUE                                                
J : ISN'T UNIQUE                                                
 OK
15 demo1

                 1
                 4
                27
               256
              3125
             46656
            823543
          16777216
         387420489
       10000000000
      285311670611
     8916100448256
   302875106592253
 11112006825558016
437893890380859375  OK

第(60)個範例是一個我個人發展系統程式時的工具，用來觀察整數及其平方後的二進制數值位元花樣。

我花過幾千塊錢台幣，買過 NC4000 Forth CPU 發展套件，後來從未用過，只是為了支持當時的 Forth 發展而買，因為台灣的易符公司下了很大的資本，每個月透過丁陳老師給 Charles Moore 五千美金，拖了好幾年，才得以有這個第一套 Forth CPU 面世，並在台灣出售，當然，隨套件附贈的文件有一大堆，我沒發展，所以也沒細看。

聽過丁陳老師講過，Forth CPU 內的開平方，是一個 CPU 的工作時脈就能完成工作的現成指令，我看不到，也沒摸成過，只把此事記在心中。

我學過數位電路，對數位運算的邏輯電路有點概念，聽丁陳老師說過，那個開平方指令，是用平移電路完成的，在 CPU 內部的震盪時脈，可以是 CPU 外在工作時脈的幾十或幾百倍，因此，善用移位電路，就能讓 CPU 以一個工作時脈執行出開平方指令。

2016 年 3 月，我從網上訊息，看到 Intel 公司公告過，新造的 CPU 改善了 division 與 sqrt 指令的執行速度，但 sqrt 指令作用於不同的暫存器時，也還仍然需要 4.5c 至 7.5c 的耗用時脈，實際運作時的潛伏期(latency)時脈，更長，從 13.1c 到 18.1c 左右。顯然，一般再好的 CPU ，也沒辦法僅耗單一個時脈，執行出開平方指令。

談這些硬體細節，可能太深奧了，也不是我這種人能接觸得了的東西，所以我僅拿一點關於硬體搞開平方的事情起頭，告訴大家為什麼我要自己設計這個範例，用來作為了解以右移(rshift)指令設計開平方程式的研究工具。

前曾述及，我長期收集過各種開平方的程式，也深入探討過這個問題，並且實際測量過各種方法的真實速度，最後，才將結論固定於自己設計的系統中。

我的系統中，隨著環境的不同，也都各有不同之整數開平方的執行辦法，不是固定的只用一套，有的時候，我只考慮自己的需要，不求速度最高，把值得留下來的技術保留在我自己設計的系統中，永久留參。

因此，以 rshift 達到開平方目的所設計程式，也仍留在我後來開發的新系統中。以軟體 rshift 模擬而成的結果，當然沒有以真正做 rshift 的硬體速度快，但速度快不是我搞研發的目的，就算我一輩子也沒機會接觸能罩刻出開平方所需元件的硬體單晶生產圖檔，只要能搞懂它的道理，不能接觸又有什麼關係？

這個範例，是為了滿足我個人求知慾望時所設計的工具，能用在模擬軟體被執行時，看到真正位元花樣的變遷狀態。我只能告訴大家， Forth 在研究軟硬體介面的過程中，確實是個好東西，這麼簡單的程式，就能看清數位訊息的存儲花樣，我用過它。

若要簡談開平方的道理，我可以講一套最簡單的暴力算法得到解答，事實上也真能這樣得到答案，就是不管怎樣，就從 1 開始自乘，當結果比所給值大時，前一個整數就是開平方的所得，只用乘法與比較邏輯。

所以，設計出開平方的程式，不是難題。程式精簡，算得快，功能必須覆蓋整個定義域，算得正確，才是要求重點。

想要純用硬體實現一個工作時脈就能執行出開平方的指令，得先搞出可以模擬硬體的軟體，然後去燒 FPGA 試一試。最後，才花錢向台積電買罩刻晶體的圖檔，生產您心目中的 CPU 。這事，我一輩子都幹不了，只能滿足自己的好奇心，以軟體模擬技術測試性能，這個範例就是研究過程中用到的工具。
:

\ (60)SquarBitServey.f

: test ( n -- )
  dup dup * >r >r
  cr decimal r@ 10 .r binary r> 20 .r
  cr decimal r@ 10 .r binary r> 20 .r
  decimal
  ;

: main (  -- )
  page
  11 0 do
  cr i test cr
  loop ;

2024年9月15日星期日

一百個例題 (51 ~ 55)

Ching-Tang Tseng

Hamilton, New Zealand

15 September 2024

第(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 時很基本的概念，我不擬浪費時間講得這麼細。

我在介紹百例超過一半後，大家應該可以體會到，這些範例，實際上有許多個，都是我在設計系統時，所必須牽扯到的技術。所以，若有人也想設計系統，就不能不搞通這些技術。想要善用新增添出來的系統功能，那就更應該仔細了解這些分項獨立出來的範例討論。我們還有一半的範例可以探討。
:

 
\ 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 個單元的陣列，範例中也有對應的宣告方式可以參考。
:

\ 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 次方也印出來外，完全一樣。

我在抄改寫舊程式時，經常會碰到這類問題，有時覺得問題好像就是原作者故意安排的陷阱，學習者若不仔細實作，根本不可能知道問題發生在那裡？

我在發展這些技術時，為了能讓自己永久留參，牢記有問題的地方，忠實留下發展記錄，絕不盲從。錯誤最能考驗系統性能的正常與否？這是這個範例被留下來的主要特色。
:

\ 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.

:

\ (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 ，對吧？試大的輸入數字就不行了。

至於計算出三角塔式的係數，是死算，是把固定算法寫成程式的範例練習。
:

\ (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)

2024年9月30日星期一

一百個例題 (76 ~ 80)

2024年9月27日星期五

一百個例題 (71 ~ 75)

2024年9月23日星期一

一百個例題 ( 66 ~ 70 )

2024年9月21日星期六

一百個例題 (61 ~ 65)

2024年9月18日星期三

一百個例題 (56 ~ 60)

2024年9月15日星期日

一百個例題 (51 ~ 55)

網誌存檔

勵志語選

2024年9月30日 星期一

一百個例題 (76 ~ 80)

2024年9月27日 星期五

一百個例題 (71 ~ 75)

2024年9月23日 星期一

一百個例題 ( 66 ~ 70 )

2024年9月21日 星期六

一百個例題 (61 ~ 65)

2024年9月18日 星期三

一百個例題 (56 ~ 60)

2024年9月15日 星期日

一百個例題 (51 ~ 55)

網誌存檔

勵志語選

2024年9月30日星期一

2024年9月27日星期五

2024年9月23日星期一

2024年9月21日星期六

2024年9月18日星期三

2024年9月15日星期日