一百個例題 (30 ~ 34)

Ching-Tang Tseng

Hamilton, New Zealand

31 August 2024

一百個例題 (25 ~ 29)

Ching-Tang Tseng

Hamilton, New Zealand

30 August 2024

一百個例題 (21~24)

Ching-Tang Tseng

Hamilton, New Zealand

26 August 2024

ilikeforth@gmail.com

第(21)個範例是能夠叫用自己的自用副程式功能展示。

前面求行列式值的範例曾經說過， Forth 天生具有叫用自己的能力，使用時的指令名稱就叫做 RECURSE 。它的工作原理也很簡單，就是在現行編譯位址內填入現正編譯指令的執行位址，就能達到目的了。如果用 Fig-Forth 的標準來解釋，最為簡單：

: RECURSE LATEST , ; IMMEDIATE

但是，這樣叫用自己的方式，與其他傳統程式語言的自用副程式執行方式，有著很大的不同。 Forth 不處理參數的傳遞，別種程式語言則在系統內耗用大量的記憶體，來進行大量參數如矩陣類的傳遞。因 Forth 只留單一個自己的起始執行位址來達到叫用自己的目的，若還要包括許多參數，就得另行設法。使用自用副程式設計程式時，外觀上的層次性影響在認知方面有點複雜，並不很受歡迎，我也不喜歡使用。所以，這個程式只是為了展示功能與用法，才將其列入範例，僅供留參。

這個範例中有三個例題，都是以叫用自己的方式設計的。第一個是求階乘的函數， 32 位元的整數系統，能表現的最大正整數，大約只能算到 12 的階乘。

使用 RECURSE 時最需要牢記的要點，就是程式中必定被分成兩個部份來處理，一個部份就是不能每次都被執行的事情，另個部份就是每次都得執行的部份。求階乘時， 0! 恆等於 1 ，是特例，就是不能每次都被執行的部份。而每次都得執行的部份又有另一個特點，就是 RECURSE 不能無限次的執行下去，否則當機。因此，程式中必定有個必須停止再度叫用自己的邏輯判斷機制，讓執行有機會能夠停止下來。這三個範例中，我都在程式後面加註了註解，仔細閱讀就能明白。

第二個例題是搬運 Hanoi 三柱塔上，只能小片壓大片之碟片搬運規矩的展示程式，這原本是一個早期討論人工智慧時，以專家系統的立場來解題的範例。程式內有個特殊設計，用到了局部變數(local variables)的技巧，所以，我刻意將其譯寫成 BASIC 的格式，考驗一下我的 ABC Forth 系統內，能不能也具有使用局部變數的功能？程式證明了，可以！所以，十二年前，我曾將這種具有局部變數與叫用自己的自用副程式能力之雙重特殊性能的 BASIC 系統，公開在國際論壇上展示貼文，引起許多人的注意。大家可以仔細看看這兩項同時存在的功能，不是其他一般 BASIC 所能具有的性能。

第三個例題是將 1 至指定數字之間，所有間隔為 5 的數字，以一五、一十的方式，全部印出來的範例。問題簡單，很容易懂，但在邊界處，要添加將 i0 指定為 0 的工作，才能正確執行，這是我們在設計程式時，經常會碰到的起始值設定問題，所以，我也將其留參。

我再強調一次，叫用自己的自用副程式功能，不是個程式語言中絕對需要的性能，我很少使用。大家喜歡強調，我就能以具體事實回應，我的系統也行，如此而已。 :

 
\ Recursive in BASIC
\ 能執行自用副程式之BASIC性能展示
\ 注意！這只是一個在BASIC環境中強行使用recurse的示範程式。
\ 有效計算範圍只在12以下有效，13以上都是錯的。

2 integers u u!

: (u!) ( -- )
BASIC
10 if u = 1 then 50              \ 結束叫用自己時，所需要的終止條件
20 let u! = u! * u
30 let  u = u - 1
40 run recurse                   \ 以上為叫用自己時，每次都必須執行的內容
50 end
;

: BFact ( u -- u! )
[[ u ]] !
BASIC
10 let u! = 1
20 if u = 0 then 40              \ 以上為不能在每次叫用自己時都執行的內容
30 run (u!)
40 end
u!
;

: Ffact ( u -- u! )
  DUP 0=
  if   DROP 1
  else DUP 1- recurse *
  then ;

\s
\ **********************************************************

\ 原 hanoi 程式將左邊(left)的疊片搬到中間(middle)，BASIC則改為left搬到right
\ 程式的特色在展示局部變數與叫用自用副程式的功能均可在BASIC式程式中使用

CREATE peg1 ." left "   
CREATE peg2 ." right "                \  peg2換成right
CREATE peg3 ." middle "               \  peg3換成middle

: .string ( addr -- )
  COUNT TYPE ;

: (hanoi)  ( n peg1 peg2 peg3 -- )
   LOCALS| vvia tto ffrom n |        \ 局部變數
   BASIC
10 if n = 1 then 30                  \ 結束叫用自用副程式的條件在此
20 goto 50
30 run cr ." Move disk from " ffrom .string ." to " tto .string
40 goto 100
50 run  n 1- ffrom vvia tto recurse  \ 可以多次使用叫用自用副程式
    =>  1    ffrom tto vvia recurse
    =>  n 1- vvia tto ffrom recurse
100 end ;

: hanoi ( n -- )
  peg1 peg2 peg3 (hanoi) ;

((
\ 原始參考程式

CREATE peg1 ," left "   
CREATE peg2 ," middle " 
CREATE peg3 ," right " 

: .string ( addr -- )
  COUNT TYPE ;

: MOVE-DISK ( n peg1 peg2 peg3 -- )
  LOCALS| vvia tto ffrom n |
  n 1 =
  IF   CR ." Move disk from " ffrom .string ." to " tto .string
  ELSE n 1- ffrom vvia tto RECURSE
       1    ffrom tto vvia RECURSE
       n 1- vvia tto ffrom RECURSE
  THEN ;

: test ( n -- )
  peg1 peg2 peg3 MOVE-DISK ;
))

\ *********************************************************

2 integers i0 i9

: (ex) ( -- )
   BASIC
10 if i0 > i9 then 100             \ 結束recurse的條件
20 let i0 = i0 + 5
30 run i0 .
40 run recurse                     \ 以上為recurse的部份
100 let i0 = 0                     \ 歸零才可以供下一次繼續使用
110 end
;

: ex ( n -- )
  [[ i9 ]] !                       \ 不可以recurse的部份
  BASIC
10 run (ex)
20 end
;

cr cr
.( Usage: ) cr
.( u Bfact u. ) cr
.( u Ffact u. ) cr
.( n hanoi ) cr
.( n ex ) cr
\ .( n peg1 peg2 peg3 MOVE-DISK ) cr

第(22)個範例是修正 Win32Forth 系統中之除法運算的程式。

最近所用的範例，都有其精彩的創作性，不要小看這些材料。

如果您知道開源式作業系統，能夠一口氣就公開提供將近 400 多個可自由叫用的功能程式時，您就不會認為第(20)個程式是個沒什麼內容的東西了。

如果您知道 Forth 與生俱來就存在有自用副程式功能，如第(21)個範例所示，您也就不會認為別人能執行得出來的程式功能 Forth 會辦不到了。

同樣的，今天要討論的範例，是一個算法上的原則問題，當您發現這麼大的一個系統，連最基本的整數除法運算都與眾不同時，怎麼辦？這個問題非常根本，我拿 Forth 研究數學計算，跟全世界的 Forth 流行趨勢不一樣，要搞得出名堂，就得深究問題，才能創作得出好的結果。

數學是一門非常嚴謹的科學，不能接受我在範例程式中列舉的運算規則，整數除法以 -15 除以 4 得到近似的 -4 而不是正確的 -3 。這個問題，是我在自行創作分數計算系統進行測試時發現的，但它並不是一個錯誤的設計，原創者有他硬要這樣設計的理由，我不想進一步說明。現在，我只強調這與數學運算基本要求不符，必須修正，因為我專注於發展數學計算應用，別種應用，我沒興趣。

發現問題之初，我僅採用高階定義方式解決問題，所以程式中用到了一個 ANSI 標準的 SM/REM 指令，它強調除完後的餘數必須與被除數同號，就能用來協助我解決問題。(相對的，另外有一個叫做 FM/REM 的指令，它強調除完後餘數必須與除數同號)。

待我完成整個分數系統的創作與測試後，我自覺以高階方式設計的指令，執行速度必定較慢，於是，再自行詳讀組合語言的說明書，將 /MOD 及 / 兩個指令，都修正成符合數學運算基本要求的低階組合語言程式了。為了避免與系統原始設計指令用名衝突，我改用 asm 前引於這兩個指令名稱之前，當作新的指令名稱，並固定於新長成的系統中。所以，一旦您再度載入這個程式時，您會發現系統告訴您指令用名重複了，但不影響後續運作。我在自創的分數運算系統中，確實使用了這個範例中的低階組合語言，設計這兩個指令。

這個範例是一個很好的 Forth 組合語言應用範例，我在網上見過名家介紹全世界組合語言工具時，竟然發現，也把 Win32Forth 中的組合語言功能，當作一個項目來表彰。但請注意！雖然這個系統中的組合語言功能，是作者 Jim Schnieder 於 1995 年捐贈出來的程式，必須注意！它有版權。公開發表用到它時的創作，就必須道謝，不能亂用它來營利，它只是一個公益捐獻。我若公開貼文用到了它，也必定會照規矩做出宣告。

搞 Forth 的人，為了解決任何所有的問題，免不了絕對必須接觸低階的組合語言，這個範例也在強調這種精神。

如果我對某個 Forth 系統有興趣，我就必定會深入使用到其最根本的運行工具：組合語言。任何 Forth 系統都能跑組合語言，有的時候，還得從源程式起，就強行使用組合語言，才能把系統搞得十全十美。我在接觸 Lina Forth 後，更有這種感覺。如果您沒有能力在 Forth 系統中使用組合語言來設計程式，那就必定會有許多問題無法解決。

這個範例教大家一丁點 Forth 的組合語言，實際上是告訴大家，必須熟悉組合語言，否則搞不出 ABC FORTH 系統。因為系統裡面，有許多地方都必須使用組合語言來解決問題。 Forth 本身就具有組合語言的使用能力，就看您怎麼用，這個範例就是一種標準用法，是一個很好的範例。

也許大家都很想學數學計算系統整套的設計方法，但我只談到了20幾個範例程式，就已經涉及了許多額外話題，數學上的算法學、系統運行的哲理、系統的性能、設計的技巧、理論問題 …. 等等題外話。這百個範例，後面幾乎都是我的個人創作，要不然，就是我改寫後的成果，改寫程式的技術，並不是件簡單事情。我並沒有把大型程式列入百例，因為太大了，就無法以短文討論，只能一天一頁。想更深入搞清問題者，必須自己執行程式，得出結果，才能全面了解。 :

 
\ (22)Asmdiv.f
\ 20140313

code asm/MOD   ( n1 n2 -- r q )
     mov   ecx, edx
     pop   eax
     cdq
     idiv  ebx
     push  edx
     mov   ebx, eax
     mov   edx, ecx
     next c;
\ : asmMOD ( n1 n2 -- r ) asm/MOD DROP ;
\ : asm/ ( n1 n2 -- q ) asm/MOD NIP ;
code asm/ ( n1 n2 -- q )
     mov   ecx, edx                    \ save UP
     pop   eax                         \ 被除數
     cdq
     idiv  ebx                         \ 除數
     mov   ebx, eax                    \ 商
     mov   edx, ecx                    \ restore UP
     next  c;

\S
Win32Forth原系統中的除法指令( / )，執行結果與數學計算之基本規則不符。
例如：執行 -15 4 / . 時會得到不合理的 -4，合理結果應該是 -3。
大部份Forth系統及C語言系統，均輸出-3而非-4，因此，必須重新設計/及/MOD。
臨時以高階定義方式完成之/指令，變通性設計如下：
: / ( n1 n2 -- q ) >r s>d r> sm/rem nip ;
參考VFX Forth系統及Win32Forth系統自身內的SM/REM內容便可以仿照設計，資料如下：
in VFX forth
see /
code / ( n1 n2 -- n3 )
     mov  eax, [ebp]
     cdq
     idiv ebx
     mov  ebx, eax
     lea  ebp, [ebp+04]
     next,
其中：
CDQ
指令為有號數除法專用指令，執行時會將EAX暫存器的最高位元值填滿整個EDX
   EAX ==> EDX:EAX
IDIV  來源運算元
指令用於有號數除法運算。上例為 IDIV EBX。
EBX放的是除數，被除數則放在EDX:EAX。
執行後，商放在EAX，餘數放在EDX。
in Win32Forth
SEE SM/REM 
SM/REM IS CODE 
  ( $401C28 8BCA )            mov     ecx, edx 
  ( $401C2A 5A )              pop     edx 
  ( $401C2B 58 )              pop     eax 
  ( $401C2C F7FB )            idiv    ebx 
  ( $401C2E 52 )              push    edx 
  ( $401C2F 8BD8 )            mov     ebx, eax 
  ( $401C31 8BD1 )            mov     edx, ecx
測試結果：
-15 -4 asm/ . 3  ok
-15  4 asm/ . -3  ok
 15  4 asm/ . 3  ok
 15 -4 asm/ . -3  ok

第(23)個程式是解常數係數一元三次方程式根的全套標準程式。

這個程式已經是一個 ABC Forth 發展到最後階段時的標準應用典範了。程式所實現的許多基本要求，足以顯示這個系統可以解決所有的數學計算問題，此範例就在驗證整套系統的性能。

數學發展到十六世紀，那是 400 多年前，義大利數學家卡丹，才解決了任意一元三次方程式以公式方式求根的問題，分析與所得公式原本很重要，後來，因為數學群論、體論的健全發展，進一步分析解決了所有更高次方程式的問題，這套公式才失去了它的重要性。但不管怎麼說，它仍然是一套很好的工具，將其寫成庫存程式，需要時，仍能派得上用場。我拜讀群書，取出解題精要，以其為例設計程式，更用它來試用系統，再把這些特性逐個提出來討論。

數學計算系統不可能容納所有的數學函數，只安裝必要、常用的函數就夠了，那麼，若有少數人自己認為必須再加裝新的常用函數時，要怎麼辦？我的系統也必須提供這種功能，所以，範例中，您可以見到一列宣告範例：

‘ fcbrt rdef cbrt

這一列宣告，就等同於傳統 BASIC 或者是 FORTRAN 系統中，用來定義出新函數的指令 DEFINE 。整數有整數的宣告方法，複數有複數的宣告方法，大整數有大整數的宣告方法，分別以 IDEF, RDEF, ZDEF, BDEF 為之。新函數可以純用 Forth 產生，也可以強用 BASIC 的設計格式產生，這種性能非常奇異，但也有用法上的哲理，探討起來會非常廣泛。我是運用函數執行時有優先度(precedence)要求的道理來形成這種功能，例如：先加減、後乘除、括弧內最優先.....等等規矩。而優先度執行秩序至少可以區分出 27 個層次以上，太多了，此處無法討論下去，您學會使用就好。

範例中，我也給出了直接隨便使用浮點運算的現成功能，以三分之一次方的方式，隨便解得開立方。註解也說明了，誤差很大。所以，我特意去 Forth 科學庫存程式庫(FORTH Scientific Library)中找出 Julian V. Noble 捐獻出來的公益程式 fcbrt ，裝入系統。沒裝入前的試用方式，就如上述。

根的判別式分析法，請自行參考數學文獻，網上很多，此處不進行討論，如果不想去找資料，範例中的文字註解，也足夠參考了。程式最後也包括了所獲得答案的驗算，得到了不等於零但數值很小的驗算結果，就代表根值是正確的。網上雲端的資料，時間久了都可能會消失，要是有興趣，現在(20240828)，去瀏覽這個中國人的網頁還看得到：

一元三次方程求根公式
https://cloud.kepuchina.cn/newSearch/imgText?id=7163950974689681408

有了理論依據，設計程式只需配合。但是，如果程式語言的性能不佳，設計起來就會很痛苦，我指的是如果強用 Forth 的後算符規矩寫程式，人人都會放棄。能有 BASIC 的程式語法功能，就能全面簡化問題，快速、正確的寫出這種程式。您自己閱讀這個程式便能了解效果，高中程度的學生，就能看得懂這個程式，我也不用多加解釋。

為了美觀，也為了強調連續使用 let 指令時可以減少打字上的麻煩，我特地將類似 30 let -----, 40 let ----- 的使用狀況，簡化成可以連續使用兩個重複的冒號 :: 來取代，這在想要於程式中直接輸入數據時非常有用。同理，我也設計了可以在連續出現類似 70 run -----, 80 run ----- 的使用狀況時，可以使用 => 來取代。

最後面的程式用法說明內有輸入係數的三種用法，一是只有一組固定輸入值時的直接測試，就用 test 指令。二是對談式用法，就用 hi 指令。三是必須進一步討論的重點，使用 main 指令，它能執行出從檔案輸入許多組數據的功能。

能執行從檔案輸入數據，是一種很重要的設計，一般人很難設計得出這種複雜性能。由於牽扯到作業系統的叫用方法，隨作業系統不同就會不同，不能通用，所以很難發展。另外，若非必要，例如：還沒有心理準備去強用檔案輸入數據前，大可不必浪費時間搞清這件事情。可是，在未來的大數據分析應用上，可能必須用到此一性能，若不事先探討清楚，來日面對問題時，就會毫無頭緒。我在 Win32Forth 系統上先行實現了這套設計，覺得有其存在的價值。這也是為什麼在這個系統稍嫌過時的當頭，我仍願不厭其煩的與大家探討 Win32Forth 及百例應用的主要原因。

試想，網路上要如何傳遞大量的數據？必須要有大家相約成俗的協定(protocol)才行。什麼是協定？如果能用單一個指令或程式解決的問題，就不需要什麼特別的協定，如果解決問題的步驟有點複雜，需要好幾個步驟才能完成，這樣解決問題的規矩，就叫做協定。大數據的傳遞，基本協定，應該是使用檔案，因為你不可能隨便亂闖別人在網路上、某一節點上的電腦，亂拿單一個或許多個數據。電腦之間的溝通使用檔案傳輸數據，是現成的傳輸功能，系統都能辦得倒。那麼，我設計 ABC Forth 系統時，該怎麼配合檔案傳輸？首先就扯到開擋、關擋的問題，只能叫用作業系統的功能程式來實現。很不幸，一般 Forth 系統雖然都提供了開、關檔的基本指令，但很少人會用，也不能直接地使用，我下了很大的功夫解決這個問題，但此處只能概述。

我設計了開檔後立刻載入文字檔案程式到一個系統內指定緩衝區(buffer)的功能，起始位址在 Fadr ，然後就立刻關檔，免得影響後續操作，也不至於影響作業系統後續的運行。一切數據，只要有個起始位址，又有幾個使用步驟與規矩作為協定，您要拿檔中數據就容易了。我設計的做法是：開檔(Get-File)後，每次只讀取一列輸入文字(GetOneLineData)。只有一個數據，就只寫一列程式處理掉此單一個數據。單一列中若有很多個數據，就寫一個定量迴路，取得整列的多個數據，看是該放進陣列，或者是該直接運算掉，都可以。這就是從檔案輸入數據的協定運作方式，我特別將其展示在那個叫做 ReadCoef 的指令內，請仔細閱讀與體會，這是系統設計時的重要精神所在。

實際上，要完成這樣的設計，並不簡單，設計內容內有高度的處理技巧與完整的哲理。我能這樣實現設計，自己想過很多問題，參考過別種傳統系統的外觀用法，也拜讀過許多系統設計名家的見解，才搞出上列我個人的協定用法，可以保證任何人都未曾在網上見過。延伸性的見解，則是我的系統能應用到現正流行的大數據分析上去。我在 V4.2 的 Win32Forth 版本中，實現過定時自動傳送檔案到網上另一端電腦內指定資料夾內去的程式設計，您看看上述協定的用法，就能明白我講的大數據分析能力，單用我設計的系統，就能辦得到。

第(23-1)CubicEqCoef.f 數據檔案內，有配合這個範例程式使用時，必須遵守的數據安放格式，也是個固定規矩，仔細看看，也可以自己添加、修改，執行後就能見到新的結果。以後，您再見到所謂的『大數據分析』術語時，您就可以了解大數據分析是怎麼運行的？必須具有這種系統設計出來的功能，才能辦得到。
:

 
\ (23)CubicEquation.f
\ Cube roots of a real coefficients cubic equation: a0x^3+a1x^2+a2x+a3=0
\ Author: Ching-Tang Tseng
\   Date: 20140321

((
\ 以下列簡易設計之導得函數fcbrt直接計算，會得到極大的誤差。
\ 實測驗證後更可凸顯出正確設計fcbrt函數的重要性。
\ 自20140430起，fcbrt已建進V654以後的版本，函數名稱固定為CBRT。

: fcbrt    ( F: N -- N^1/3)
  fabs 1e 3e f/ f** ;

: fcbrt    ( F: N -- N^1/3)
        fdup  f0= if fdrop 0e0 exit then
        FDUP  F0<  FABS                 ( F: -- |N|)  ( -- f)
        FDUP  FSQRT                     ( F: -- N x0 )
        BEGIN
        FOVER FOVER
        FTUCK   FDUP F* F/  FSWAP  F2*  F+  3E0  F/
        FTUCK
        F-   FOVER   F/   FABS   1.0E-8  F<
        UNTIL
        FTUCK   FDUP F* F/  FSWAP  F2*  F+  3E0  F/
        IF   FNEGATE   THEN
;

' fcbrt RDEF cbrt
))

\ Cubic Equation: a0x^3 + a1x^2 + a2x + a3 = 0
\ 1.三實根    1, -3, 5
\ 2.三重實根  3
\ 3.一實根，共軛複數根

11 reals a0 a1 a2 a3 q r s t R1 Zr Zi

: GetCoefficients
basic
\ 10 let { a0 = 1 } :: { a1 = -3 / a0 } :: { a2 = -13 / a0 } :: { a3 = 15 / a0 }
\ 10 let { a0 = 1 } :: { a1 = -9 / a0 } :: { a2 = 27 / a0 } :: { a3 = -27 / a0 }
10 let { a0 = 1 } :: { a1 = 20 / a0 } :: { a2 = 600 / a0 } :: { a3 = 1200 / a0 }
20 end
;

2 reals x f(x)

: CheckRealRoot
basic
10 let { f(x) = a0 * x * x * x + a1 * x * x + a2 * x + a3 }
20 end
;

7 complexs z za1 za2 za3 f(z) f(z1) f(z2)

: CheckComplexRoot
basic
10 let [ za1 = r>zr ( a1 ) ] :: [ za2 = r>zr ( a2 ) ] :: [ za3 = r>zr ( a3 ) ]
20 let [ f(z) = z * z * z + za1 * z * z + za2 * z + za3 ]
30 end ;

: CheckConjugateComplexRoots
basic
10 let [ z = r>zr ( Zr ) + r>zi ( Zi ) ]
20 run CheckComplexRoot
30 let [ f(z1) = f(z) ]
40 let [ z = r>zr ( Zr ) - r>zi ( Zi ) ]
50 run CheckComplexRoot
60 let [ f(z2) = f(z) ]
70 end ;

: Compute
basic
10 let { q = ( 3 * a2 - a1 * a1 ) / 9 }
    :: { r = ( 9 * a1 * a2 - 27 * a3 - 2 * a1 * a1 * a1 ) / 54 }
    :: { s = cbrt ( r + sqrt ( abs ( q * q * q + r * r ) ) ) }
    :: { t = cbrt ( r - sqrt ( abs ( q * q * q + r * r ) ) ) }
    :: { R1 = s + t - a1 / 3 }
    :: { Zr = negate ( ( s + t ) / 2 + a1 / 3 ) }
    :: { Zi = abs ( ( sqrt ( 3 ) / 2 ) * ( s - t ) ) }
20 end
;

4 reals dlt dlt1 dlt2 sigma

\ Discriminatory analysis
\ delta > 0 one real root, two conjugate complex roots 一實根，二共軛複數根
\ delta = or < 0 three real toots     三個實根

: delta
basic
10 let { dlt1 = ( ( a1 * a1 * a1 ) / ( 27 * a0 * a0 * a0 ) )
              + ( a3 / ( 2 * a0 ) )
              - ( a1 * a2 ) / ( 6 * a0 * a0 ) }
    :: { dlt2 = ( a2 / ( 3 * a0 ) )
              - ( a1 * a1 ) / ( 9 * a0 * a0 ) }
    :: { dlt  = ( dlt1 * dlt1 ) + ( dlt2 * dlt2 * dlt2 ) }
20 print { " delta = " ; dlt }
30 end
;

: delta>0
basic
10 run cr
20 print
   " delta > 0 means there are one real root and two conjugate complex roots."
30 run cr
40 let { x = R1 }
50 print { " Root x1 = " ; x }
60 run CheckRealRoot
70 print { " f( x1 ) = " ; f(x) }
80 run cr
90 print { " Root x2 = " ; Zr ; " + " ; Zi ; " i" }
100 run CheckConjugateComplexRoots
110 print [ " f( x2 ) = " ; f(z1) ]
115 print { " ABS(f(x2)) = " ; ZABS ( f(z1) ) }
120 run cr
130 print { " Root x3 = " ; Zr ; " - " ; Zi ; " i" }
140 print [ " f( x3 ) = " ; f(z2) ]
145 print { " ABS(f(x3)) = " ; ZABS ( f(z2) ) }
150 run cr
160 end
;

: delta=<0
basic
10 run cr
20 print { " delta = or < 0 means there are three real roots. " }
30 run cr
40 print { " Root x1 = " ; R1 }
50 print { " f( x1 ) = " ; f(x) }
60 run cr
70 let { x = Zr + Zi }
80 print { " Root x2 = " ; x }
90 run CheckRealRoot
100 print { " f( x2 ) = " ; f(x) }
110 run cr
120 let { x = Zr - Zi }
130 print { " Root x3 = " ; x }
140 run CheckRealRoot
150 print { " f( x3 ) = " ; f(x) }
160 run cr
170 end
;

: CubeRoots
basic
10 run cr
20 run Delta
30 run Compute
40 if { dlt > 0 } then 70
50 run delta=<0
60 goto 80
70 run delta>0
80 end
;

: test
basic
10 run GetCoefficients
20 print " Typic cubic equation is: a0x^3 + a1x^2 + a2x + a3 = 0"
30 print { " a0 = " ; a0 ; " a1 = " ; a1 ; " a2 = " ; a2 ; " a3 = " ; a3 }
40 run CubeRoots
50 end
;

: hi
basic
10 print " Typical cubic equation is: a0x^3 + a1x^2 + a2x + a3 = 0 "
20 print " Please enter 4 coefficients a0 a1 a2 a3 "
30 inputr a0 , a1 , a2 , a3
40 let { a0 = a0 / a0 } :: { a1 = a1 / a0 }
    :: { a2 = a2 / a0 } :: { a3 = a3 / a0 }
50 run CubeRoots
60 end
;

2 integers i n
20 4 matrix Coef

: ReadCoef
basic
10 run S" (23-1)CubicEqCoef.f" Get-File
20 run GetOneLineData
30 let n = INT ( BRA ( 1 ) )
40 for i = 1 to n
50 run GetOneLineData
60 let { coef ( i 1 ) = BRA ( 2 ) }
70 let { coef ( i 2 ) = BRA ( 3 ) }
80 let { coef ( i 3 ) = BRA ( 4 ) }
90 let { coef ( i 4 ) = BRA ( 5 ) }
90 next i
100 end
;

: main
basic
10 run ReadCoef
20 run cr
30 for i = 1 to n
40 print " Typic cubic equation is: a0x^3 + a1x^2 + a2x + a3 = 0"
50 print " For given data set( " ; i ; " )"
60 let { a0 = Coef ( i 1 ) } :: { a1 = Coef ( i 2 ) }
    :: { a2 = Coef ( i 3 ) } :: { a3 = Coef ( i 4 ) }
70 print { " a0 = " ; a0 ; " a1 = " ; a1 ; " a2 = " ; a2 }
80 print { " a3 = " ; a3 }
90 let { a0 = a0 / a0 } :: { a1 = a1 / a0 }
    :: { a2 = a2 / a0 } :: { a3 = a3 / a0 }
100 run CubeRoots
110 print " ===================================================================="
120 next i
130 end
;

cr cr
.( 程式用法: ) cr
.( 1. test : 固定的輸入數據時使用。 ) cr
.( 2. hi   : 交談式輸入數據時使用。 ) cr
.( 3. main : 由檔案輸入數據時使用。 ) cr

\S
for 8 B/FLOAT ABC V654 system
18 sigdigits !
main

Fileid is: 1080 
96 Bytes are read into the Fadr RAM buffer!

Typic cubic equation is: a0x^3 + a1x^2 + a2x + a3 = 0
For given data set( 1 )
a0 = 2.00000000000000000 a1 = 0.00000000000000000E0 a2 = 1.00000000000000000
a3 = 4.00000000000000000

delta = 4.03703703703703744 

delta > 0 means there are one real root and two conjugate complex roots.

Root x1 = -1.37879670012954912 
f( x1 ) = 1.15463194561016272E-14 

Root x2 = .689398350064774528 + 1.55750128578313152 i
f( x2 ) = -5.32907051820075136E-15 - 1.26565424807267840E-14 i  

Root x3 = .689398350064774528 - 1.55750128578313152 i
f( x3 ) = -4.44089209850062656E-15 + 1.28785870856518160E-14 i  

====================================================================
Typic cubic equation is: a0x^3 + a1x^2 + a2x + a3 = 0
For given data set( 2 )
a0 = 1.00000000000000000 a1 = 9.00000000000000000 a2 = 5.00000000000000000
a3 = 1.00000000000000000

delta = 5.62962962962967616 

delta > 0 means there are one real root and two conjugate complex roots.

Root x1 = -8.42030108523868416 
f( x1 ) = -2.13162820728030048E-14 

Root x2 = -.289849457380658 + .186407861904320 i
f( x2 ) = -1.33226762955018784E-15 + 2.22044604925031328E-16 i  

Root x3 = -.289849457380658 - .186407861904320 i
f( x3 ) = -1.55431223447521920E-15 - 1.11022302462515664E-16 i  

====================================================================
Typic cubic equation is: a0x^3 + a1x^2 + a2x + a3 = 0
For given data set( 3 )
a0 = 1.00000000000000000 a1 = 2.00000000000000000 a2 = 3.00000000000000000
a3 = 4.00000000000000000

delta = 1.85185185185185152 

delta > 0 means there are one real root and two conjugate complex roots.

Root x1 = -1.65062919143938784 
f( x1 ) = 1.77635683940025056E-15 

Root x2 = -.174685404280306 + 1.54686888723139616 i
f( x2 ) = 1.33226762955018784E-15 + 0.00000000000000000E0 i  

Root x3 = -.174685404280306 - 1.54686888723139616 i
f( x3 ) = 2.22044604925031328E-15 + 0.00000000000000000E0 i  

====================================================================
Typic cubic equation is: a0x^3 + a1x^2 + a2x + a3 = 0
For given data set( 4 )
a0 = 1.00000000000000000 a1 = -3.00000000000000000 a2 = -13.0000000000000000
a3 = 15.0000000000000000

delta = -151.703703703703680 

delta = or < 0 means there are three real roots. 

Root x1 = 1.00000000000000000 
f( x1 ) = 1.77635683940025056E-15 

Root x2 = 5.00000000000000000 
f( x2 ) = 0.00000000000000000E0 

Root x3 = -2.99999999999999936 
f( x3 ) = 1.42108547152020032E-14 

====================================================================
Typic cubic equation is: a0x^3 + a1x^2 + a2x + a3 = 0
For given data set( 5 )
a0 = 1.00000000000000000 a1 = -9.00000000000000000 a2 = 27.0000000000000000
a3 = -27.0000000000000000

delta = 0.00000000000000000E0 

delta = or < 0 means there are three real roots. 

Root x1 = 3.00000000000000000 
f( x1 ) = 1.42108547152020032E-14 

Root x2 = 3.00000000000000000 
f( x2 ) = 0.00000000000000000E0 

Root x3 = 3.00000000000000000 
f( x3 ) = 0.00000000000000000E0 

====================================================================
Typic cubic equation is: a0x^3 + a1x^2 + a2x + a3 = 0
For given data set( 6 )
a0 = 1.00000000000000000 a1 = 20.0000000000000000 a2 = 600.000000000000000
a3 = 1200.00000000000000

delta = 4982222.22222222208 

delta > 0 means there are one real root and two conjugate complex roots.

Root x1 = -2.13581907909972992 
f( x1 ) = -1.36424205265939232E-12 

Root x2 = -8.93209046045013504 + 21.9559354517858656 i
f( x2 ) = -2.72848410531878464E-12 + 1.81898940354585632E-12 i  

Root x3 = -8.93209046045013504 - 21.9559354517858656 i
f( x3 ) = -1.81898940354585632E-12 - 1.81898940354585632E-12 i  

====================================================================
Typic cubic equation is: a0x^3 + a1x^2 + a2x + a3 = 0
For given data set( 7 )
a0 = 1.00000000000000000 a1 = -5.00000000000000000 a2 = 8.00000000000000000
a3 = -4.00000000000000000 

delta = 1.04083408558608432E-17 

delta > 0 means there are one real root and two conjugate complex roots.

Root x1 = 1.00000000000000000 
f( x1 ) = 0.00000000000000000E0 

Root x2 = 2.00000000000000000 + 0.00000000000000000E0 i
f( x2 ) = 0.00000000000000000E0 + 0.00000000000000000E0 i  

Root x3 = 2.00000000000000000 - 0.00000000000000000E0 i
f( x3 ) = 0.00000000000000000E0 + 0.00000000000000000E0 i  

==================================================================== ok

(23-1)CubicEqCoef.f 的內容：

7
1 2 0 1 4
2 1 9 5 1
3 1 2 3 4
4 1 -3 -13 15
5 1 -9 27 -27
6 1 20 600 1200
7 1 -5 8 -4

第(24)個程式是求解整數的立方根之範例程式。

這個範例，最早在 16 位元的時代，我就開發了。進展到 32 位元時，也還能使用。直到 64 位元的時代，被全面仔細的測試過後，才發現了許多問題，結論是不宜採用。

凡是浮點函數，大部分都能對應的寫出整數形式的函數，有些很有用，有些就不太有用。例如：整數的開平方，整數的對數函數，都很有用。但整數的開立方，就不太有用。我在發展十進制浮點系統時，也發現使用整數開方的方式來獲得對應之浮點數的開方很不可行。

這個範例中也寬列了使用浮點數方式獲得開立方、然後轉換成整數來獲得答案之程式寫法。我用它們來測試純整數程式的執行結果，那個 test 指令，就在做這件事情。結果，在很多地方，發現都不對勁，最後乾脆不用了，內容則適合作為練習寫程式時的教材，並展示開發系統時都做過什麼事情。

測試系統是件很重要的工作，我在前一個範例(23)程式中，也列有驗算方式的測試程式，從測試中，您可以發現，許多地方的驗算結果，都沒有得到絕對值應該是 0 的答案，反而印出了許多 10 的負十幾次方的微小數值，您就可以了解，實際工程與學術理論上的許多差異。在學校裡，老師教您、考您的題目，數字答案都很整齊，那種題目，是業經挑選過的理想安排，實際工程中，無數的問題，都是像我在範例中解得的結果一樣，答案數字的花樣都亂七八糟。所以，測試程式非常重要，能輕鬆設計驗算程式時，不要省略，多做無妨。

設計這套 cbrtf1 的方法，採用牛頓法(Newton’s method)，是憑牛頓法開平方的設計原理推衍出來的，我若不解釋引用的理論式子，您可能不容易看懂，今天這篇貼文就做此事。

程式的開頭處，我列出了推衍出來的兩個根據式子，乍看之下，那有開立方的意義？我解釋一下，您就能明白。

恆等式所描述的意義是：如果 A 是 n 的開立方所得，那麼，這兩個式子，都能用來寫出程式，一直算到獲得了 A’ ，就是解答。

認知式子時，要這樣看： 當 n 剛好為 A 的三次方時，您把 A 的三次方代入上面兩個式子中的 n ，最後就能化簡出 A’ = A 的結果，這就是牛頓法的根據。

第一個式子雖然也能用來設計出程式，但因為要用到除以 3 ， 這在只能做整數除法的狀況下，誤差必大。也沒有除以 3 的移位技巧，可以用來優化程式。所以，還是選用第二個恆等式設計程式比較好。列出第一個式子的目的，是表示推衍時用過這個參考式子。

在程式結束前，我使用了一個 \s 取消了兩個比較好的另套設計，來求開平方與求開立方，您可以刪除 \s 指令，也讓這兩套設計可被載入系統進行測試。這個後來新添加的 cbrt ，有個特點，也是用了牛頓法，收斂很快，但可以不經邏輯判斷來結束環路，而是採用強行執行定量迴路獲得結果，在 32 位元中，只需計算 19 次，就絕對能收斂出答案，這個 19 次與後面註解說明的 64 位元時為 36 次，是我實測後獲得的答案。

這些資料留存在此處，均僅供參考，不宜正式使用。我在全面落定系統設計前，陸陸續續又發現這些用法都有問題。有些是在 0 或最大值情況，算得的答案是錯的，但您不能設計一個系統去強迫使用者避開這些數目不用。有些錯誤是在工作範圍內就出現循環而不收斂，例如：就在從 0 算到 100 的全定義域內全都測試時，發現在 86.5 附近，只收斂到小數點後面三位數就開始循環了，得到的答案當然不對。

最後強調：我在設計我的浮點系統時，採取了唯一且絕對收斂的『二分法』設計出開方函數，不採用這些設計。 :

 
\ (24)整數之cbrt.f

\ A' = ( ( n / A * A ) + 2 * A ) / 3
\ A' = ( ( n / A * A ) + A ) / 2

0 value Turnning
0 value Result

: cbrtf1 ( n1 -- n2 )
  0 to Turnning  0 to Result
  1                                 \ n1 A
  begin
        >r dup r@ r@ * / r@ + 2/    \ n1 A'    rs: A
        r> 2dup >=
        if   1 +to Turnning   to Result
        else drop
        then                        \ n1 A'
        Turnning 19 >  \ at least 5 times for condition = in 32 bits system
  until
  2drop Result
;

((
: fcbrt1   ( f: |x| -- |x|^[1/3] )
  fln 3e0 f/ fexp ;

: cbrtf1
  s>f fcbrt1 f>s ;
))

: F**2  
  FDUP F* ;

3 S>D D>F FCONSTANT F=3

: X' ( F: N x -- x')
  FTUCK   F**2  F/  FSWAP  F2*  F+  F=3  F/  ;

\ The magic number 1E-8 needs no change, even when extended (80-bit) precision
\ is needed.
: CONVERGED?    ( F: x' x x' --) ( -- f)
  F-   FOVER   F/   FABS   1.0E-8  F<   ;

: FCBRT2    ( F: N -- N^1/3)
  FDUP  F0<  FABS                 ( F: -- |N|)  ( -- f)
  FDUP  FSQRT                     ( F: -- N x0 )
  BEGIN    FOVER FOVER  X'   FTUCK   CONVERGED?   UNTIL
  X'   IF   FNEGATE   THEN  ;

: cbrtf2
  s>f fcbrt2 f>s ;

0 value cnt

: test
  0 to cnt
\  2147483647 2140000000
  9000000 8000000
  do
\     cr i .
     i cbrtf1 i cbrtf2 <>
     if 1 +to cnt
        cr ." i = " i .
        ." , cbrtf1= " i cbrtf1 dup . ." , i^3 = "  dup dup * * .
        ." , cbrtf2= " i cbrtf2 dup . ." , i^3 = "  dup dup * * .
     then
  loop
  cr cnt . ." errors!"
  cr ." finished!" ;

\s
\ 以此牛頓法開平方程式為根據： A' = ( n / A + A ) / 2
: sqrt9 ( square -- root )
  DUP 1 <
  IF DROP 0 EXIT THEN
     dup                            \ 起始猜測A為n
     begin
           2dup / over + 2/         \ 計算(N/A+A)/2
           >r
           dup r@ >                 \ A>(N/A+A)/2 ?
     while
           drop r>                  \ A>(N/A+A)/2時，則A換成(N/A+A)/2。
     repeat
     nip r> drop                    \ 當A<=(N/A+A)/2時，A即為所求。
;
\ 開發出下列牛頓法開立方程式： A' = ( ( n / A * A ) + A ) / 2
: cbrt ( n1 -- n2 )
  DUP 1 <                   \ These two lines guard against a negative input.
  IF DROP 0 EXIT THEN
  1                          \ n --- un, for n>0 only.
  19 0 DO                    \ 64 位元的系統須為至少36，32位元的系統19便可。
            >r dup r@ r@ *   \ n n A*A
            / r> +  2/       \ n (n/(A*A)+A)/2, 注意！A若為0，n/A會不能執行。
       LOOP                  \ n A'
       DUP DUP DUP * *       \ n A' A'^3
       >R SWAP R> U<         \ A' f           check
       IF ABS 1- THEN        \ aquare root,   correct it, if needed.   
  ; 

2024年8月15日 星期四

一百個例題 (11~20)

Ching-Tang Tseng

Hamilton, New Zealand

16 August 2024

ilikeforth@gmail.com

第(11)個程式是亂數產生程式。

產生的方法採用同餘法則(congruence method)，數學表示式的新亂數 U 為：

U(n+1) = [(magic number) * U(n) + (relatively prime)] mod M

從程式中可以看出，十六位元的最大正整數 M = 32767 而
magic number = 31421
relatively prime = 6927

關於同餘法則的原理與亂數評鑑技術，此處不深究，有興趣者請自行查找資料進行了解。

這個程式是 16 位元時代的產品，使用舊 CPU 單板電腦的人，還會有此需要。這個程式精簡扼要，也有利於用來講解亂數。進入 32 位元時代後的亂數產生程式，仍可以採用此法則來設計。有興趣的人，也可以前往 Forth 科學程式庫(FSL)找到編號為 #60 由 Julian V. Noble 捐獻的浮點系統亂數產生程式進行研究。

進入 64 位元時代後，我不再採用同餘法則來設計亂數產生程式，改採組合語言設計的位元炒作法(bits manipulation method)，左右平移少量位元後作邏輯比較，以如此方式炒動幾次位元，來得到新亂數，速度較快，也容易驗證與評鑑。

另外，種子(seed)亂數是個起始值，放在變數 RND 內。新時代的 CPU 旁邊有個時鐘式的自動計時器(timer)，隨時可以直接取得其值，就能提供一個很好的硬體隨機起始種子設定值，我設計系統時，已經這樣使用。這個程式的起始值，是編譯這個亂數時，系統當時耗用記憶體的位址 HERE 之值。
:

	
\ Example 11. Random numbers

VARIABLE RND                            ( seed )
HERE RND !                              ( initialize seed )

: RANDOM ( -- n, a random number within 0 to 65536 )
        RND @ 31421 *                   ( RND*31421 )
        6927 +                          ( RND*31421+6926, mod 65536)
        DUP RND !                       ( refresh he seed )
        ;

: CHOOSE ( n1 -- n2, a random number within 0 to n1 )
        RANDOM UM*                      ( n1*random to a double product)
        SWAP DROP                       ( discard lower part )
        ;                               ( in fact divide by 65536 )

\s
Usage:
To test the routine, type

        100 CHOOSE .
        100 CHOOSE .
        100 CHOOSE .	

第(12)個程式，是整數的開平方程式。

目前介紹的幾個範例程式，雖然都必須被列為整數系統的基本函數，但在一般能夠免費獲得的 Forth 原始系統中，都不被建成固有函數。不搞數學計算的人，可以不用這些設計，搞數學計算的發展者，才需要這些設計。這些設計也絕非一成不變的東西，常常需要因時因地而制宜，今天介紹的這個開平方程式範例也是如此。

學會這些函數的設計，是培養設計數學計算系統能力的先決條件。想設計系統，必須先知道函數的設計技術，我的 ABC FORTH 系統內，就都是以設計這類函數程式的本領起家，但各套的內容未必相同。可能算法不同，格式不同，語法不同，還可能非用低階組合語言設計不可。搞系統，就得直接面對這些方法。這也就是我並不太喜歡談『如何建造我的系統？』之主因，誰願意花那麼大的功夫去搞通十幾種方法設計出單一個函數的道理與優劣？這個開平方程式也是，我手頭上收集、試用、評鑑、修改過十幾套開平方的程式。這年頭，已經不會有人想看了。介紹百例的貼文，以引介 Forth 為始，以介紹曾用過的工具及函數為輔，終結於使用 ABC FORTH 的格式來發展應用程式。

我在完成十進制浮點系統(decimal floating point system)前，浮點的開平方函數，確實曾採用過雙整數的開平方方式設計系統，16 位元時沒有關係，32 位元時也還勉強能用，64 位元以後，就不能用了，為什麼？我把 32 與 64 位元能表示的整數範圍列在此處，您就能明白，就算現在的電腦跑得這麼快，也無法令人滿意的接受某些古時候的開平方函數算法及其程式。

在 Win32Forth 系統中，測出單整數及雙整數能表示之最大範圍的方法，可執行下列指令來獲得：
-1 1 rshift . 2147483647 ok
-1 -1 ud. 18,446,744,073,709,551,615 ok

在 64 位元的 Forth 系統中，執行同樣指令，您可以得到這樣的顯示：
-1 1 rshift .
9223372036854775807 OK
-1 -1 ud.
340282366920938463463374607431768211455 OK

你看看，算開平方時，不管你執行下一次回算能跳過多少增量，每次都要從 0 開始測試，然後，測試它是否大於迴路設定的上限了？ 是，就得到答案了，不是，就還得增量再算。今天的範例程式，就在做這件事情，核算的關鍵點，在 +loop 指令被執行時增量是多少？前面函數的輸入待求值，會被轉放到回返堆疊上去，環路執行到 +loop 時，每次的增量會是 2n+1 。就算算到後面，增量隨著 n 的增大而會越來越大，什麼時候才能到達上列的天文數字最大值？ 我測試過了，會耗時很久。系統不能要求使用者只能用它來開小數目字的平方，大數目字的平方不管。所以，這個方法在 64 位元的環境中，就不能再用了。用來當作生產主機板時的預燒測試(burn-in test)程式，也許還能用用。

Forth 名家 Wil Baden 曾經給過全世界的勸告是：always use FSQRT. 意思是就算您只想算整數的開平方，最佳的方法，仍是使用硬體數學運算處理器。先求得浮點數的開平方後，才將其轉換成單整數來使用，方法最好，速度最快。我實測過，所言不假，但是，不是所有使用 Forth 的場合，都有現成的 coprocessor 可用，Forth 愛用者，還是必須自行解決這種問題，這個範例程式就必須留下來。

單論開平方的算法，也有很多種，為了能讓計算的收斂速度加快，另有許多其他算法技術，例如：牛頓法、二分法...等等方法，可以列出十幾種。各種方法除了速度快慢上的優劣點外，也都還有無法解決的問題，例如：奇點、發散、或陷入循環等等數學上的理論問題。我有完成過這些研究後寫成的長文，我想，現代人不看長文了，是否貼出？以後看狀況再說吧。

貼文不適合再談下去，延伸問題會有很多，因此，就此打住。不熟 Forth 的人，應該趁著在此處的討論，拿範例當教材，學會自己流暢使用 Forth 設計基本指令的程式方法。如果只是看看貼文就算了，不去實際使用 Forth ，我想，最後的結果，必定還是不會用 Forth 。我傳播 Forth ，我也每天只用 Forth ，而且是只用我發展妥當後的 Forth 來解決所有的電腦程式問題。學而時習之，能喚醒許多自己原本已培養成的技術，對我自己有益。
:

	
\ Example 12. Square root

: SQRT ( n1 -- n2, n2**2<=n1 )
        0                               ( initial root )
        SWAP 0                          ( set n1 as the limit )
        DO      1 + DUP                 ( refresh root )
                2* 1 +                  ( 2n+1 )
        +LOOP                           ( add 2n+1 to sum, loop if )
        ;                               ( less than n1, else done )

\s
Usage:

n sqrt .

第(13)個程式是最大公因數(GCD)與最小公倍數(LCM)的函數程式。

此例題是一個很好的結構化指令群 BEGIN ----- WHILE ----- REPEAT 之使用範例，其中也展示了邊界條件為 0 時就不執行迴路的 ?DUP 指令之恰當用法，原始程式係參考了 gForth 系統內的設計。

程式的算法，採用輾轉相除法，大家都熟，所以沒有什麼好解釋的。但學一學如何將『輾轉相除，直到餘數為0時，就能得到答案』的口語，將其寫成程式的方法，則是一項很好的程式設計練習題目。目前，所有的程式語言，都沒有把這兩個函數固定成啟動系統時就是現成可用的函數，因為它們很容易設計，一列程式就能完成。

我積極納入這兩個整數數學函數將其固定於系統的目的，是為了要實現一項全世界都還沒有人實現過的新型數學體系計算能力，也就是『繁分數的可程式化』計算能力。這樣的系統，能對所有的分數進行四則運算的計算，而且能寫成程式。這樣的系統完全根基於這兩個函數 ： GCD 與 LCM 。

在我公告我的系統性能之前，全世界只能找到真分數的運算系統，也就是例如：1/2 + 1/3 = 5/6 的運算能力。但我們的實際生活中，與小學的教材內，搞分數，就得搞到繁分數的四則運算，才會有實際用途。我拿這個問題來研究最基礎的數學系統，體會獨立數學體系的設計方法，最後，培養出了自己能夠設計出任何獨立數學體系進行獨立運算系統的能力。後來也設計出了無限位數的浮點運算系統，所有能力都源自於此。這些系統的源程式，均未曾對全世界公佈過，由於內容龐大，數學體系種類繁多，想要以反向工程倒推出源程式來，已不容易。

現在給大家的系統，裡面有五個獨立的數學體系：整數、分數、浮點數、複數、大整數。我沒有把大浮點數加載進去，但可以在要用的時候才加載。以後，若有其他數學體系如基因運算數學體系、或量子運算數學體系被開發出來時，也可以同樣的方法載入系統。我所謂的體系運算功能，不是僅止於純粹只用傳統 Forth 設計出來的體系就算了，必須加增能力成為能用 ABC Forth 格式寫成程式，才是我的奮鬥目標。我舉實例說明：

例如 : 只做兩個繁分數的除法，並印出結果。

繁分數就包括了所有種類的分數，它需要三個整數單元來表示一個元素，例如： 1 又 1/2 要寫成 ( 1 1 2 )，可以不用括號，其餘類推。印出時的格式則是 1 & ( 1 / 2 )。你可以直接試著執行：

1 1 2 /3./

看看結果。

凡是分數體系的程式或指令，都必須使用一對斜線 /.../ 來包裹。我寫個簡單的小例子，展示一下。
:

	
3 fractions x y z

: test 
Basic
10 let /{ x = ( 3 1 2 ) }/
20 let /{ y = 1 3 7 }/
30 let /{ z = x / y }/
40 run z /3./
50 end 
;

您請自行試試看，驗證執行 test 後顯示出來的結果，會印出 2 & ( 9 / 20 )。

您若不怕打錯字，可以在螢幕上直接輸入程式，或乾脆 copy 此處的程式後，直接貼進您的系統，就能看到結果。

我沒介紹這些功能前，您沒見過能這樣執行這種數學體系運算的系統吧？請特別注意系統的可程式化性能，也包括能夠正確處理負的繁分數。系統不接受不合理的輸入，例如：( 3 1 0 )，因為，分母不得為 0。

這就是可程式化計算的能力。請注意！這些能力必須基於先行設計出正確的 GCD 與 LCM 後，才能完成。繁分數的運算過程與印出能力，都要大量用到這兩個整數函數。 :

	
\ Example 13. 	The Greatest Common Divisor (GCD) and Least Common Multiple (LCM)

: GCD ( n1 n2 -- n3 )
  BEGIN
      OVER MOD
      ?DUP
  WHILE
      SWAP
  REPEAT
  ABS ;

: LCM ( n1 n2 -- n3 )
  2DUP GCD */ ABS ;

\s
Usage:

n1 n2 GCD .   
n1 n2 LCM .

第(14)個程式用來算出 Fibonacci 級數的數字。

Fibonacci 級數的口語定義為：下一個數字，等於前兩個的和。它的意義、起源與用途，請自己查找網上資料。

範例程式中， Fib1 的程式還可以讓您看到一個加號計算，然後印出所有小於 50000 的一序列級數。不熟 Forth ，就很難感受出下一個是前兩個之和的計算是怎麼安排的？數學問題，採用純粹傳統 Forth 的語法寫程式時，因為只能用堆疊傳遞數字，必然影響這種直覺。

Fib2 根據前引輸入的數字，把所有低於此指定數字的 Fibonacci 數字印出來。程式中更看不出那裡才在做加法計算了，計算是隱性的，在執行 +loop 之後，系統會為迴路指標 I 執行加法，加出來的新指標量，會改放到回返堆疊去，不放在數據堆疊上。印出的數字，是 I 指標從回返堆疊取得的下標數字，這裡玩弄了一點程式技巧。
:

\ Example 14.	The Fibonacci Sequence

: Fib1 ( -- , print all Fibonacci numbers less than 50000 )
        1 1                             ( two initial Fib numbers )
        BEGIN   OVER U.                 ( print the smaller number )
                SWAP OVER +             ( compute next Fib number )
                DUP 50000 U>            ( exit if number too large )
        UNTIL                           ( else repeat )
        2DROP                           ( discard the numbers )
        ;

: Fib2 ( n -- , display all Fibonacci numbers smaller than n )
        1                               ( initial number )
        SWAP 1                          ( set up range )
        DO      DUP U.                  ( print current number )
                I                       ( the next Fibonacci number )
                SWAP                    ( prepare the next to come )
        +LOOP                           ( add current to index, if the )
                                        ( repeat until sum>n )
        U.                              ( print the last Fib )
        ;

\s
Usage:
To test the routines, try:

      Fib1
10000 Fib2

我只需發幾分鐘時間，就能根據口語的定義，寫出能夠執行出同樣結果的 ABC FORTH 程式， Fibonacci 的定義式子就表現在列號為 60 的程式中。比較程式的寫法後，就能感受得出我為什麼要創作 ABC FORTH 來寫數學計算程式了。 :

 
5 integers i n f(n) f(n+1) f(n+2)

: Setn ( -- )
basic
10 print " n = "
20 inputi n
30 end ;

: main
basic
10 run Setn cr
20 let f(n) = 1
30 run f(n) .
40 let f(n+1) = 1
50 run f(n+1) .
60 let f(n+2) = f(n) + f(n+1)
70 if f(n+2) > n then 1000
80 run f(n+2) .
90 let f(n) = f(n+1)
100 let f(n+1) = f(n+2)
110 goto -60
1000 end
;

取材自 16 位元環境中的範例，只把上限數字算到 65535 以下，最多只有 24 個。同樣的程式，在 32 位元環境中可以執行到大約第 46 個 Fibonacci 數字，可為十八億以上。在 64 位元系統中可以得到 91 個，使用者可以自己測試一下以體會可用數字的範圍。

更甚的是，我提供給大家的系統，具有大數計算能力，就拿這個 Fibonacci 級數舉例，它能夠算出第一萬個 Fibonacci 數字是多少？ 可用規模大約有四萬七千多個，遠大於 91 個。

程式受輸出位數最多印出 10008 位數的限制，所以下列程式的使用上限之 N 才會大約是四萬七千多個。

程式及執行結果如下：
:

 
\ Fibonacci numbers


2 INTEGERS I N
0 BIGVARIABLE F(N)   10000 ALLOT
0 BIGVARIABLE F(N+1) 10000 ALLOT
0 BIGVARIABLE F(N+2) 10000 ALLOT

: SetN ( -- )
   BASIC
10 PRINT " Enter n =: "
20 INPUTI N
30 END ;

: MAIN ( -- )
   BASIC
10 RUN SetN
20 LET B{ F(N)   = I>BIG ( 1 ) }B
30 LET B{ F(N+1) = I>BIG ( 1 ) }B
40 IF N < 3 THEN 100
50 FOR I = 1 TO N
60 LET B{ F(N+2) = F(N) + F(N+1) }B
70 LET B{ F(N)   = F(N+1) }B
80 LET B{ F(N+1) = F(N+2) }B
90 NEXT I
100 PRINT " Fibonacci number( " ; N ; " )="
110 RUN CR F(N) BIG.
120 END ;

cr cr
.( Usage: main ) cr

main

Enter n =:
? 10000

Fibonacci number( 10000 )=

2090 digits
54438373113565281338734260993750380135389184554695   :50
96702624771584120858286562234901708305154793896054   :100
11738226759780263173843595847511162414391747026429   :150
59169925586334117906063048089793531476108466259072   :200
75936789915067796008830659796664196582493772180038   :250
14411588410424809979846964873753371800281637633177   :300
81927941101369262750979509800713596718023814710669   :350
91264421477525447858767456896380800296226513311135   :400
99297627266794414001015758000435107774659358053625   :450
02461707918059226414679005690752321895868142367849   :500
59388075642348375438634263963597073375626009896246   :550
26687461120417398194048750624437098686543156268471   :600
86195620146126642232711815040367018825205314845875   :650
81719353352982783780035190252923951783668946766191   :700
79538847124410284639354494846144507787625295209618   :750
87597272889220768537396475869543159172434537193611   :800
26374392633731300589616724805173798630636811500308   :850
83967495871026195246313524474995052041983051871683   :900
21623283859794627245919771454628218399695789223798   :950
91219943177546970521613108109655995063829726125384   :1000
82420078971090547540284381496119304650618661701229   :1050
83288964352733750792786069444761853525144421077928   :1100
04597990456129812942380915605503303233891960916223   :1150
66987599227829231918966880177185755555209946533201   :1200
28446502371153715141749290913104897203455577507196   :1250
64542523286202201950609148358522388271101670843305   :1300
11699421157751512555102516559318881640483441295570   :1350
38825477521111577395780115868397072602565614824956   :1400
46053870028033131186148539980539703155572752969339   :1450
95860798503815814462764338588285295358034248508454   :1500
26446471681531001533180479567436396815653326152509   :1550
57112748041192819602214884914828438912417852017450   :1600
73055389287178579235094177433833315068982393544219   :1650
88805429332440371194867215543576548565499134519271   :1700
09891980266518456492782782721295764924023550759555   :1750
82056475693653948733176590002063731265706435097094   :1800
82649710038733517477713403319028105575667931789470   :1850
02411880309460403436295347199746139227479154973035   :1900
64126330742308240519999961015497846673404583268529   :1950
60388301120765629245998136251652347093963049734046   :2000
44510636530416363082366924225776146828846179184322   :2050
4793434406079917883360676846711185597501
=========1=========2=========3=========4=========5
12345678901234567890123456789012345678901234567890
=========(c) 2018 Copyright, Bottom Ruler=========
 OK

這麼大的 Fibonacci 數字有什麼用？

有人用它來產生亂數。理由是：請看顯示出來的數字結果，幾乎完全沒有規則可言，夠亂，隨便取一段或一截夠用的數字，都不會與別處取得的相同。如果光靠範例程式所能展示的性能，看不出這樣的結果。這也是我創作這套 ABC FORTH 的主因。
:

第(15)個程式是用來印出美國資訊交換標準碼(ASCII code)的範例程式。

它有什麼用？回顧第(2)個範例，我曾把大寫的 Y 碼修改成小寫的 y 碼後，程式才能正常執行出所有的結果。這就是它的典型應用例子，我曾用這個程式來印出 ASCII code　，進行比對，然後才仔細研究出原始程式的問題所在。

Forth 強調能夠執行即時控制(real time control)。單板電腦上裝 Forth ，開機啟動後，就能直接下命令執行對外控制。個人電腦上裝的 Forth ，也是啟動 Forth 後就能對事物進行控制。沒加裝特殊硬體時怎麼控制？就用鍵盤。所以會有 KEY , ?KEY 等指令可供即時控制需要之時使用。

鍵盤的存在，自古以來就很多變，但有個慣例不變，就是大家相約成俗的只使用 ASCII code 編碼。一般而言，初期發展 Forth 系統時，系統對大小寫嚴格區分。後來，大家覺得使用起來很不方便，於是，就將系統設計成大小寫不分，系統接受輸入後，就自行將大小寫的輸入都轉換成大寫的碼來處理。　Win32Forth , Lina64 都是，但只對執行比較、找到指令時的工作內容是這樣，做其他事情時，則還是保持對鍵盤輸入大小寫不同時，字碼嚴格區分。

至於後來出現文字編碼有所謂 BIG-5 或 UTF-8 編碼稱呼上的不同，放心，那只是用在電腦螢幕輸出顯示字型的編碼而已，對現行 Forth 系統的運作，不會有很大的影響，也就是鍵盤碼跟顯示字型，沒有直接的關係。如果您寫程式時，會在檔案中加註中文說明，那在轉換作業系統時，就會出問題。檔案的顯示碼必須轉換，才能看得到正確的文字顯示。克服這種問題的方法也不麻煩，只需一列作業系統操作指令就能解決，這裡就不討論。

這個範例程式設計成能將 ASCII code 對應碼印成表格，一個橫印，另個直印，鍵盤的控制碼與特殊鍵的碼，則不印出來。以後，必須學會如何以操作方式，獲得這些鍵的碼。以前，在國際網壇上，曾有貼文詢問如何獲得滑鼠鍵碼的問題，以便設計出能夠控制顯示 Google 地圖的應用程式，很久都沒有人回答，每當憶及這些苦況，心中常有不忍。後來，我在『行星齒輪模擬程式』那篇文章中就故意披露方法，大多數人可能不太會去注意。

一個長期使用 Forth 的人，會在電腦旁準備幾樣東西，令其垂手可得、利於參考。以我為例，我準備有 ： F79 , F83 與 ANSI 標準指令卡、 F79 與後來標準不同指令的對照表、能快速查得十進制與十六進制的 ASCII code 表、我自己創作之系統的使用說明.....等等等，這些東西，都是設計程式時的必備工具，經常要用。

還有，我也花了不少錢，印出所有發展後落定程式的源碼，才不怕破壞，也利於查核。中國有句俗話說：絹八百，紙一千。您就可以明白，無論電腦儲存媒體的發展有多高明，就是不如中國人講的紙一千，就是說紙張的資料能長保一千年的意思，千年之後再印成紙，又能有再一千年的保存能力。錄音帶、軟碟片、硬碟機、光碟片、USB 的 flush memory ，雲端庫存，全都沒有這種能力，戰亂一起，還都可能一切歸零。印成紙張，永恆的存活率最高。

這個程式在我們外出與好友談 Forth 技術時，如果實體表不在身邊，網上雜物一大堆，找到適當的好表也很浪費時間時，不妨就用這個程式直接印出對照表來使用。

我曾將幾百張 ANSI 標準指令表轉交給協會，只要願意踏入 Forth 行列者，都應該可以索討此表使用，永久留參。我搞出過這麼多東西後，知道自己根本不可能永遠記得住所有創作出來的使用規矩，所以，我自己也經常要參考自己出版的『使用手冊』，一輩子都要用。這個範例程式的意義也是這樣。

:

 
\ Example 15. 	ASCII Character Table

: Printable ( n -- n , convert non-printable characters to spaces )
        DUP 14 <                        ( 7-13 are special formatting )
        IF      DUP 6 >                 ( characters not displayable )
                IF DROP 32 THEN         ( substitute them by space )
        THEN
        ;

: HorizontalASCIItable ( -- )
        CR CR CR
        5 SPACES
        16 0 DO I 4 .R LOOP             ( show column header )
        CR
        16 0 DO                         ( do 16 rows )
                CR I 16 * 5 .R          ( print row header )
                16 0 DO                 ( print 16 characters in a row )
                        3 SPACES
                        J 16 * I +      ( current character value )
                        Printable
                        EMIT            ( print it )
                LOOP                    ( loop for next character )
        LOOP                            ( loop for next row )
        CR
        ;
: ht horizontalasciitable ;

: VerticalASCIItable ( -- )
        CR CR CR
        5 SPACES
        16 0 DO I 16 * 4 .R LOOP        ( show column headers )
        CR
        16 0 DO                         ( do 16 rows )
                CR I 5 .R               ( print row header )
                256 0 DO                ( do 16 columns )
                        3 SPACES
                        J I +           ( current character )
                        Printable EMIT
                16 +LOOP                ( skip 15 characters between columns)
        LOOP
        CR
        ;

: vt verticalasciitable ;

\s
Usage:

Type
	HorizontalASCIItable 
	VerticalASCIItable
or
ht
vt

第(16)個程式是根據輸入之年月日算出該日為星期幾的程式。

從現在開始，一百個範例，將進入一個新的階段，除非特別需要，後續範例程式，只要涉及數學計算，都將盡量以 ABC FORTH 的格式設計出程式。

這個程式本身的可讀性就已經相當高了，實在不需要特別解釋，但其中卻蘊含了高度的設計哲理。如果我不指點出來說明，您並不容易看得出來，因此，我只挑這方面的特點來強調。

首先，這個程式是所有利用到儒略日技術的總結，只有這種程式的寫法，能夠直接告訴大家儒略日的換算公式為何？程式就是答案，不必贅述。

其次，在正式走入 ABC FORTH 環境前，我最好把整件事情的重點，先行提出來列示於此，以利大家能夠快速進入狀況。

貼示在我的個人網頁 20180616 的網文：井田制度(# field farming system)中，曾經告訴全世界，我的系統可以精簡到只用這 11 個指令，便能寫出所有的數學計算程式。 :

 
1. LET     計算
2. FOR     就從
3. TO      做到
4. STEP    間隔
5. NEXT    配合
6. GOTO    跳到
7. IF      假如
8. THEN    前往
9. RUN     執行
10. PRINT  印出
11. END    結束

進一步講，只有 10 個，因為，最後的那一個 END 是強迫結束 BASIC 式環境回到 Forth 式環境所需要使用的指令。後續可以空白，也可以再添寫 FORTH 程式，如此而已。

只有 10 個要學、要記，教小學生都很容易了，所以有其特殊價值。

右邊對應列出來的中文名稱，只有 LET 被改稱為『計算』，與最原始的使用手冊內容不同，在此特別聲明。

我在長期研發系統時，考慮永久性的問題，也把將來會發展出來的語音輸入能力列入考慮，我拿 Google 的環境測試語音輸入，結果，就是這個 LET 指令，原用『令』為中文名稱時，成為無法解決的問題。中文的單個字，會有許多其他的單個字與其同音，但少有兩個字以上形成的詞會同音。所以，我就開始把 LET 改譯成『計算』，取其意義與這個指令執行的目的相關來定其名。第(16)個程式中同時列示了可執行的對應中文程式，可以看出，把『令』改成『計算』，不會影響整體精神，反而更符合程式的實際語意。試想，這個指令，若用在武器系統，讓操作者直接下口語命令要電腦執行，這時，就別再在那裡『令』、『另』難分，搞不出結果了，用『計算』就不會有此困擾。
:

 
\ 全英文版程式

5 INTEGERS julianDate WhatDay Year Month Day

: PrintWhatDay ( -- )
BASIC
10 IF WhatDay = 0 THEN 90
20 IF WhatDay = 1 THEN 110
30 IF WhatDay = 2 THEN 130
40 IF WhatDay = 3 THEN 150
50 IF WhatDay = 4 THEN 170
60 IF WhatDay = 5 THEN 190
70 IF WhatDay = 6 THEN 210
80 GOTO 1000
90 PRINT " This date is Sunday."
100 GOTO 1000
110 PRINT " This date is Monday."
120 GOTO 1000
130 PRINT " This date is Tuesday."
140 GOTO 1000
150 PRINT " This date is Wednesday."
160 GOTO 1000
170 PRINT " This date is Thursday."
180 GOTO 1000
190 PRINT " This date is Friday."
200 GOTO 1000
210 PRINT " This date is Saturday."
1000 END 1000
;

: MAIN   ( -- )
BASIC
10 PRINT " Enter the date as Year Month Day: " CR
20 INPUTI Year , Month , Day
30 LET julianDate
       = ( 367 * Year )
       + ( ( Month * 275 ) / 9 ) + Day + 1721029
       - ( ( ( ( Month + 9 ) / 12 ) + Year ) * 7 / 4 )
       - ( ( ( ( Year - ( ( Month + 9 ) / 12 ) ) / 100 ) + 1 ) * 3 / 4 )

40 LET WhatDay = ( JulianDate + 1 ) MOD 7
50 RUN PrintWhatDay
60 END
;

cr cr
.( Usage: main ) cr

\S
\ 全中文版程式

5 個整數變數 星期餘數 儒略日 年 月 日

: 印出星期幾   ( -- )
中文程式
10 假如 星期餘數 = 0 前往 90
20 假如 星期餘數 = 1 前往 110
30 假如 星期餘數 = 2 前往 130
40 假如 星期餘數 = 3 前往 150
50 假如 星期餘數 = 4 前往 170
60 假如 星期餘數 = 5 前往 190
70 假如 星期餘數 = 6 前往 210
80  跳到 1000
90  印出 " 這一天是星期天。"
100 跳到 1000
110 印出 " 這一天是星期一。"
120 跳到 1000
130 印出 " 這一天是星期二。"
140 跳到 1000
150 印出 " 這一天是星期三。"
160 跳到 1000
170 印出 " 這一天是星期四。"
180 跳到 1000
190 印出 " 這一天是星期五。"
200 跳到 1000
210 印出 " 這一天是星期六。"
1000 結束
;

: MAIN   ( -- )
中文程式
10 印出 " 請以 年 月 日 秩序輸入日期，數字之間必須空格 : " CR
20 輸入整數 年 , 月 , 日
30 計算 儒略日
        = ( 367 * 年 )
        + ( ( 月 * 275 ) / 9 ) + 日 + 1721029
        - ( ( ( ( 月 + 9 ) / 12 ) + 年 ) * 7 / 4 )
        - ( ( ( ( 年 - ( ( 月 + 9 ) / 12 ) ) / 100 ) + 1 ) * 3 / 4 )

40 計算 星期餘數 = ( 儒略日 + 1 ) MOD 7
50 執行 印出星期幾
60 結束
;

第(17)個程式是計算出任意階行列式值的程式。

數學問題一旦進展到矩陣的格式，就可算是將問題發展到了極致。矩陣常由聯立方程式形成，將問題聯立方程式化，等同於將問題進行系統性的綜合分析，所以說他算是極致。

處理矩陣數學的問題中，最基本的計算問題，就是算出矩陣的行列式值(determinant)，此值有很多地方要用到它。手算三階以下的行列式值，還算容易辦到，四階以上，就很難使用手算獲得，這個範例因此而安排一個很規則、但很簡單的四階行列式當輸入，來測示程式。但是，我測過，要算 100x100 的行列式也不會有問題，就算要應用到現在最熱門、最流行的雲端大數據處理問題，這個程式也能使用。

測試的方法有點講究，搞矩陣數學者都知道，相依聯立方程式無解，這種矩陣的行列式值等於 0。因此，必須找出線性獨立的聯立方程式來測試程式才有意義。我長期研究這方面的問題，拜網路之賜，我很快就發現了一種所謂的魔術方陣(magic matrix)之產生方法，確定這種矩陣絕對不會線性相依，於是先自行設計出一套它的產生程式，再將輸出結果交給這個範例程式測試，任何高階的魔術方陣都被證明確實是屢試不爽。可產生魔術方陣的程式，已經公佈在我的個人網頁中，請自行查找。所以，大家可以永久放心的使用這個範例。

有許多人喜歡拿程式語言有無自我叫用的(recursive)性能來強調語言的好壞，我則不以為然，能不用此性能時，絕不使用。 Forth 可以執行 recursive 式的程式，但參數的傳遞相當困擾，就算能夠執行，也會把使用者搞得昏頭轉向、不清不楚。另外，我深入探索過 Forth ，深知凡是執行 recursive 程式時，系統都必須將前一層次的相關參數，推放到一個堆疊上去，每解完一層後，才再取回上一層的參數，進行相同的遞迴式處理。堆疊容量都是有限的，耗用太多堆疊的程式，都不是正規的、也不是很好的使用方法，危險又混亂，所以我能不用 recursive 時就盡量不用。

不幸，後來發展出來的許多程式語言，都常只用 recursive 法求解行列式值，網上很多，我不想用，只好回頭找舊書，只找到了這篇可貴的資源：

P.N.Corkett and J.D.Tinsley, Practical Programming, Cambridge University Press, 1968, p.182

是 ALGOL 語言寫成的程式，有流程圖，很容易了解解法，但不容易寫成書中漂亮的程序，如果我不直接抄改原 ALGOL 程式，我承認，自己很難創作出這麼漂亮的程式。抄改過程中，我也才發現，我自己創作的 ABC FORTH 程式語法，特別適用於這種抄改古老程式的工作，凡是古聖先賢寫成的 Algol, PL1, FORTRAN-77, BASIC, PARSCAL 程式，幾乎都能直接抄改。

簡述計算出行列式值的算法為：

利用矩陣運算理論中的線性定理，將某一列之值同乘上一個常數再加到另一列去，矩陣本身的性質不變，這樣的轉換運算也適用於兩列之間的關係，這樣的特性就被稱為線性轉換關係。亦即與『所有的解乘上任意常數後相加，所得結果亦為其解』的意義相同(這也是現行最熱門流行之量子儲存擁有線性疊加性能理論的依據)。根據這個定理，我們就可以想辦法把任何一個矩陣，轉換成對應的上右三角、左下三角、右下三角、左上三角之等效矩陣，本質不變。結果，最後就能單只把對角線上的元素之值連乘起來，這就是行列式值。這個範例中轉換成的，是右上三角矩陣。

計算很冗長，別忘了它是在處裡矩陣，古聖先賢很精明，知道算法會弱化精確度，所以算前先找出最大元素，改放到對角線的最左上方，一切以他為準，來算出能消去半個三角區域內之元素所需乘上的理想常數，才進行計算，就會比較不失精確。此外，原作者也熟知電腦表現數字時有微量的不準度，該為零時可能會不為零，必定影響後續的計算處理，所以，當計算結果小於最低精確度時，必須強迫該元素為零。這些程序，這個程式中都做了，我已全用英文加好了註解，留供永久參考。

這麼漂亮、這麼複雜的程序，憑良心講，我自己是寫不出來的，所以程式很可貴。類似這種快要失傳的程式寫法，數量非常多也非常龐大，而且都很有用，不抄改寫進入我的 ABC FORTH 系統來庫存就覺得實在可惜。我常看舊書、看舊程式，發現有很多這種東西，目前還能獲得，日後就很難講。

我已不指望能有人共襄盛舉、一起做抄改寫的工作，因為要讀懂數值分析的書本，要自行認識數學問題，再找出可貴的古老創作，才能下手工作。這樣要求別人太難了，只能自己有空、有興趣時、自己做。

我放著許多這方面的資料無法看完，只能看一點做一點，選定書本後，持之以恆的有系統進行，產品已有一大堆，只是不宜放進百例，但這種工作，恆為有意義的工作，樂趣無窮。

 
\ 計算任意階行列式值程式2009-08-23

INTEGER N
4 4 MATRIX AA

: INPUT-DATA
  [[ N = 4 ]]
  {{ AA ( 1 1 ) = 1 }} {{ AA ( 1 2 ) = 2 }} {{ AA ( 1 3 ) = 3 }}
  {{ AA ( 1 4 ) = 4 }}
  {{ AA ( 2 1 ) = 2 }} {{ AA ( 2 2 ) = 3 }} {{ AA ( 2 3 ) = 4 }}
  {{ AA ( 2 4 ) = 1 }}
  {{ AA ( 3 1 ) = 3 }} {{ AA ( 3 2 ) = 4 }} {{ AA ( 3 3 ) = 1 }}
  {{ AA ( 3 4 ) = 2 }}
  {{ AA ( 4 1 ) = 4 }} {{ AA ( 4 2 ) = 1 }} {{ AA ( 4 3 ) = 2 }}
  {{ AA ( 4 4 ) = 3 }}    ;

INTEGER I INTEGER J INTEGER K
INTEGER P INTEGER Q
REAL A REAL D

: SHOW-AA    BASIC
10 FOR I = 1 TO 4
20 FOR J = 1 TO 4
30 PRINT { AA ( I  J ) }
40 NEXT J
50 NEXT I
60 END ;

: (DET) BASIC
\ A program to evaluate a determinant by the leading diagonal method,
\ using largest pivots.
\ 10 RUN INPUT-DATA
20 LET { D = 1 }
30 FOR K = 1 TO ( N - 1 )
40 LET P = K
50 LET Q = K
60 LET { A = ABS ( AA ( K K ) ) }
70 FOR I = K TO N
80 FOR J = K TO N
90 IF {  ABS ( AA ( I J ) ) > A } THEN 110
100 GOTO 140
110 LET P = I
120 LET Q = J
130 LET { A = ABS ( AA ( I J ) ) }
140 NEXT J
150 NEXT I
\ End of search for largest element.
160 IF P <> K THEN 180
170 GOTO 230
180 FOR J = K TO N
190 LET { A = AA ( P J ) }
200 LET { AA ( P J ) = AA ( K J ) }
210 LET { AA ( K J ) = NEGATE ( A ) }
220 NEXT J
\ End of interchange of rows P and K.
230 IF Q <> K THEN 300
240 GOTO 270
250 FOR I = K TO N
260 LET { A = AA ( I Q ) }
270 LET { AA ( I Q ) = AA ( I K ) }
280 LET { AA ( I K ) = NEGATE ( A ) }
290 NEXT I
\ End of interchange of columns Q and K.
\ Largest element is now the pivot.
300 LET { D = D * AA ( K K ) }
310 IF { ABS ( D ) < 1.0e-10 } THEN  430
320 FOR I = K + 1 TO N
330 FOR J = K + 1 TO N
340 LET { A = AA ( K J ) * AA ( I K ) / AA ( K K ) }
350 LET { AA ( I J ) = AA ( I J ) - A }
360 IF { ABS ( AA ( I J ) ) < ABS ( A ) * 1.0e-10 } THEN 380
370 GOTO 390
380 LET { AA ( I J ) = 0 }
390 NEXT J
400 NEXT I
\ End of reduction to upper triangular form.
\ 405 PAUSE
410 NEXT K
420 LET { D = D * AA ( N N ) }
430 IF {  ABS ( D ) > 1.0e-10 } THEN 450
440 LET { D = 0 }
450 LET { D = D }
\ Determinant is got in D.
\ 460 PRINT { D }
470 END ;

: DET BASIC
10 RUN INPUT-DATA
20 IF N = 1 THEN 40
30 GOTO 60
40 LET { D = AA ( 1 1 ) }
50 GOTO 110
60 IF N = 2 THEN 80
70 GOTO 100
80 LET { D = AA ( 1 1 ) * AA ( 2 2 ) - AA ( 1 2 ) * AA ( 2 1 ) }
90 GOTO 110
100 RUN (DET)
110 PRINT { D }
120 END ;


cr cr
.( Usage : ) cr
.( det ) cr

第(18)個程式是利用 DO ----- LOOP 指令印出以指定字元顯示的一個方塊。

這個程式沒有數學計算，所以不必強用 ABC Forth 格式寫程式。 使用 ABC Forth 時的基本精神是：使用時，並不應該因為有了 ABC Forth ，就得什麼程式都要用 BASIC 的格式寫程式。

我們在閱讀『FORTH入門』(Startiing FORTH)時，最先學會的程式寫法，就是使用 42 emit 在螢幕上輸出一個星號。那樣教 Forth ，也是在告訴大家，FORTH 的輸出能力非凡，比當時的任何程式語言都好，不只是控制螢幕輸出，連自動控制的輸出也都是這樣。因為，電腦的自動控制，只是在將一種位元花樣，放到指定記憶體的位址而已，其意義，與控制螢幕輸出是同一個道理。

我在收集範例程式時，考慮到螢幕輸出的範例，必定陸陸續續的會出現新發展與新構想，所以，只要與最原始的螢幕輸出相關的範例，我就收集，取其具有漸進式的教育效果，我才將其列為範例。後面還有更高級的螢幕輸出範例，可以輸出很複雜的字元圖，例如：第(85)個範例，可以以字元繪製出很複雜的 Mandelbrot 圖型。

有些範例，雖在前面已經用過，但發展到後來，我發現像寫文章一樣，把程式寫得更有規格，會更易於永久有用。於是，就再把前面已經用過的範例，重寫後再度編入百例。後面，您可以看到前面介紹的求行列式值程式的另種編排格式，又出現在第(92)個範例中，請不要覺得奇怪，新格式的範例程式，有其用途。

這個範例程式並不難看懂，不必深入解釋。但請實際執行一下，看看結果。您便可以發現螢幕輸出時，執行速度並不快，這一要點，倒是值得強調。

這個程式中的字元輸出，是逐個送出去顯示的，所以速度很慢。我在設計必須印出很多位數的數字時，如果採用每算出一個數字就執行一次印出，那麼，速度就會變慢。但是，如果一次就將所有位數的計算算完，待印之字元都暫時放在一個指定的緩衝區內暫存，印出時則一次印出，速度就會很快。

這是說明這個範例時值得特別指出性能差異的教材。
:

 
\ (28)ASCIIBox.f
\ 以指定之字元，印出指定寬與高之長方形，

: star 42 emit ;

: top    ( width -- )
  0 do star loop cr ;
: bottom ( widty -- )
  top ;
: middle ( height -- )
  star
  2 - 0 do space loop
  star cr ;

: box    ( width height -- )
  cr over
  top
  2 - 0 do dup middle loop
  bottom ;

page cr cr
.( Usage Example: ) cr
.( 40 10 box ) cr

第(19)個程式是一個求 pi 值的程式。

發展數位計算技術，都會拿求質數與求 pi 的問體考驗自己，我也不例外，題目的種類繁多，還時有新創意被發表出來，追蹤它可以忙得您不亦樂乎看都看不完。

這一範例，只能討論重點。

從大學起，我讀的是 1969 年第三版 MIT 教授 Thomas 的微積分教科書， chapter 17.6, p.591, 談到了 Computation of pi ，書中列有計算出 pi 的詳細公式，由於是運用 Taylor 級數展開式獲得的結果，我還能了解，並用程式來實現計算出無限位數的成果。許多別種無窮級數公式的計算法，也能寫成程式，但我也沒真正地做過。

算 pi 的方法有很多種，目的則只有一個，就是想獲得『快速收斂的計算』技術，以便在求取非常多的位數時派得上用場。這個範例使用一種稱之為水龍頭滴水(spigot)法的演算法求 pi 。他的基本公式定義為下式：
:

轉化成下式來寫程式。

pi = SUM(k = 0 --> OXO) (1/16^k) ( 4/8K+1 - 2/8k+4 - 1/8k+5 - 1/8k+6 )

上式中，我用 OXO 表示無限大，不深究公式。但是，式中有幾個要點，能解釋演算法的意義。所謂 spigot 算法，就是利用類似水龍頭滴水的方式，演算時，每算完一次，就只留最後的位數，供下次演算時使用，前次的結果就留在記憶體中供最後印出結果時使用。式中因為出現有 16, 8, 2 等係數，因此，演算時，可以將數字二進制、八進制、十六進制化，讓滴水的效果就是取進制後最後一位的餘數。這樣算，收斂速度會很快，至少比 Taylor 展開式的算法快，如此而已。

算法中有個根據設定位數宣告出來的大陣列，對陣列內容進行處理的指令，就是在作上列公式的 spigot 式計算。

這個範例的原設計者，是 Lina64 Forth 系統的設計人，荷蘭的高能物理學家 Albert ven der Host ，但程式被印度裔的加拿大高能物理學家 kForth 的設計人 Krishna Myneni 修改成他的 kForth 系統能執行得出來的程式。我取得之後，再度修改成我的 ABC Forth 可以執行得出結果的現在範例。我的本業也算是屬於高能物理，我跟 kForth 的設計人熟識，他們兩人也都曾在網上國際論壇，直接點名推介過好幾次我的程式。改出這個可用程式，似乎就是我們的共同喜好。

程式中列示的 http://shootout.alioth.debian.org/ 網頁，不用去逛了，因為撤網了，早就已經不存在。

程式中，凡是該列後面標有背斜線外加三個星號者，如 \ *** something ，都是被我修改過的地方，這樣， Win32Forth 系統才能跑得出結果。修改的技術，主要是呈現在如何解決局部變數(Local variables)的用法，在 Win332Forth 系統中，凡是指令之後以大括弧 { …… } 括起來的部份，就是宣告出局部變數的用法。說實在的，我很不喜歡使用局部變數設計程式，我有我自己的方法，能解決這方面的問題，此處暫且不表，為了能共容於別人，我只好這樣修改程式來適應大家。

另外，為了輸出上的美觀，我也利用了自己創作系統的『右側邊標』與『底部字元尺』，以每列形成五十、一百.....的方式印出結果，一般的 Forth 系統沒有這樣的功能。

這個範例算是 1995 年以後才出現的數學技術應用，很不幸，後來已有人證明了這個算法的不可用性。因為有人運用了上千台電腦連線跑過很多天，證明了算到 10 的 15 次方位數後(一千兆次以後），全部都為 0 ，證明了這個方法不可用。

大家看看我的修改、認識局部變數怎麼用、體會如何快速執行出結果.....等等的現象就夠了。我則耗費過時間，研究過這種冗長問題。
:

 
\ pi.f
\ 改寫自kforth之程式，只能run一次。故編譯完成後直接run。
\ 最後修正日期：20140114

\ Great Computer Language Shootout
\ http://shootout.alioth.debian.org/
\ contributed by Albert van der Horst, Ian Osgood

\ read NUM from last command line argument
\ 0. argc @ 1- arg >number 2drop drop constant NUM
\ 儲存計算結果的資料結構，靠這個常數來決定，編譯後就不能改，故用常數。
1300 constant NUM                                   \ ***

\
\ Arbitrary precision arithmetic
\ A p-number consists of a count plus count cells, 2-complement small-endian
\

\ Shorthand.
: p>size ( pn -- size ) POSTPONE @ ; IMMEDIATE
: p>last ( pn -- msb ) DUP p>size CELLS + ;
: [I] POSTPONE I POSTPONE CELLS POSTPONE + ; IMMEDIATE

\ Give sign of p
: p0< ( p -- flag ) p>last @ 0< ;

\ Copy a p-number to another buffer
: pcopy ( src dst -- ) OVER p>size 1+ CELLS MOVE ;

\ Check for overflow, extend the p-number if needed
: ?carry ( carry p -- ) 2DUP p0< <> IF 1 OVER +!  p>last ! ELSE 2DROP THEN ;

\ In-place multiply by an unsigned integer
: p* { n p -- }
  p p0<  0. ( sign dcarry )
  p p>size 1+ 1 DO
    p [I] @       ( digit )
    n UM* D+ SWAP ( carry digit )
    p [I] ! 0
  LOOP
  ROT n UM* D+ DROP  p ?carry ;

\ Ensure two p-numbers are the same size before adding
0 value signn                                  \ *** 取消二度使用局部變數
: extend { p n -- }                            
  p p0<  to signn \ { sign }                   \ *** sign to signn
  p p>size 1+  n p +!
  p p>size 1+ SWAP DO signn p [I] ! LOOP ;     \ ***
: ?extend ( p1 p2 -- p1 p2 )
  OVER p>size OVER p>size - ?DUP IF
    DUP 0< IF >R OVER R> NEGATE
    ELSE OVER SWAP
    THEN extend
  THEN ;

\ In-place addition of another p-number
: p+  ?extend { src p -- }
  src p0< p p0<  0. ( sign sign dcarry )
  p p>size 1+ 1 DO
    p   [I] @ 0 D+
    src [I] @ 0 D+ SWAP
    p   [I] ! 0
  LOOP
  DROP + + p ?carry ; \ add signs, check for overflow

\ In-place subtraction of another p-number
: p-  ?extend { src p -- }
  src p0< p p0<  0. ( sign sign dcarry )
  p p>size 1+ 1 DO
    p   [I] @ 0 D+
    src [I] @ 0 D- SWAP
    p   [I] ! S>D
  LOOP
  DROP + + p ?carry ; \ add signs, check for overflow

\
\ pi-spigot specific computation
\

\ approximate upper limit on size required (1000 -> 1166)
NUM 2* CELLS constant SIZE

\ Current z transformation
CREATE aq 1 , 1 , SIZE ALLOT
CREATE ar 1 , 0 , SIZE ALLOT
    \ "as" identical zero and remains so
CREATE at 1 , 1 , SIZE ALLOT

\ Generate non-zero parts of next matrix ( K 4K+2 2K+1 )
VARIABLE KK                                                       \ *** K to KK
: generate ( -- q r t ) 1 KK +!   KK @  DUP 2* 1+  DUP 2* SWAP ;  \ ***

\ HERE is used as a temporary p-number

\ Multiply z from the left
: compose< ( bq br bt -- )
  DUP at p*  ar p*  aq HERE pcopy  HERE p*  HERE ar p+  aq p* ;

\ Multiply z from the right
: compose> ( bt br bq -- )
  DUP aq p*  ar p*  at HERE pcopy  HERE p*  HERE ar p-  at p* ;

\ Calculate z at point 3, leaving integer part and fractional part.
\ Division is by multiple subtraction until the fractional part is
\ negative.
: z(3)  ( -- n pfract )HERE  aq OVER pcopy  3 OVER p*  ar OVER p+
  0 BEGIN SWAP at OVER p-  DUP p0< 0= WHILE SWAP 1+ REPEAT ;

\ Calculate z at point 4, based on the result for point 3
\ and decide whether the integer parts are the same.
: z(4)same? ( pfract -- flag ) aq OVER p+  p0< ;

: pidigit ( -- nextdigit)
    BEGIN z(3) z(4)same? 0= WHILE DROP generate compose< REPEAT
    1   OVER 10 *   10   compose> ;

: .digit ( -- ) pidigit [CHAR] 0 + EMIT ;

\ : .count ( n -- ) .\" \t:" 1 U.R CR ;    \ ***
: .count ( n -- )  3 spaces ." :" 1 U.R CR ;

\ Spigot n digits with formatting
: spigot ( digits -- )
  cr 0                                     \ ***
  BEGIN 50 +  2DUP > WHILE                 \ ***以下全原為10改成50
    50 0 DO .digit LOOP  DUP .count        \ ***
  REPEAT
  2DUP 50 - DO .digit LOOP  OVER - SPACES  .count ;   \ ***

((
: TypeCountingRuler ( -- )                 \ *** 自建底標尺
  cr ." =========1=========2=========3=========4=========5"
  cr ." 12345678901234567890123456789012345678901234567890"
  cr ." ========(c) 2014 Copyright, Counting Ruler========"
  cr ;
))

NUM spigot TypeCountingRuler \ bye         \ *** 取消bye

第(20)個程式是一個用來展示叫用系統功能以便獲得數據的的範例程式。

這個程式的來源是 yahoo 國際論壇網頁中的 Q＆Ａ 資料。這個網頁原本是一個活動尚稱熱絡的 Win32Forth 國際論壇網頁，後來因為貼文者亂來，不顧世人的反感，拿著別人的創作當自己的作品亂貼一通，結果引起全世界的公憤，這個國際論壇就被唾棄而沒落到快要消失了。在這個網頁中作亂的人，都在網上盜取別人的設計，猛講說：網上能夠公開獲得的東西都沒有版權，難怪所有 Forth 好手都要唾棄他們。網上公開的東西確實是可以引用，尊重創作者並不是難事，引用時告知讀者出處是該做的事情。

此範例之 Q&A 的回答者是個有修養的荷蘭人叫 Jos ，他在國際論壇中很活躍，他也有自己的個人網頁，是 Win32Forth 系統應用發展方面的長期貢獻者，對系統也很了解。我只在國際論壇上與他對過話而已。

舉這個範例，取其簡單，但不太有意義。提出問題者問：如何取得視窗的長、寬設定值？ Jos 隔了幾天才答覆說，可以利用叫用作業系統的功能程式來獲得。我試 run 後，顯示所言不假，但是，不是我很想要的東西，但能當簡單範例來留參。

我曾再度在 W10 中試 run ，結果仍正確，我獲得了 1366x768 的設定結果，輸出結果隨電腦的不同會不同， Asus 的 XP , hp 的 W7 , Acer 的 W10 ，答案都會不一樣，您執行後的結果，可能也會與我的結果不一樣。

第(19)、(20)兩個程式，都只能在 Windows’ 作業系統下使用，不能在 Linux 環境中使用。而這個程式使用專有的 call 指令來叫用作業系統的功能程式，一個編號是 0 ，另個是 1 ，如果您把 SM_ＣXSCREEN 直接輸入，再用 dot ( . ) 印出結果，您就可以看到是 0 放在堆疊上。也就是說： call 指令前面需要一個作業系統功能程式的編號，就能正確地執行出結果。結果則可以執行 GetSystemMatrics 指令而獲得。這就是一個典型最簡單的作業系統功能程式叫用規則。 Win32Forth 系統也仿效常用的 C 語言，把編號都賦予了慣用的名稱。 Win32Forth 系統內也因此而擁有成千上萬個常數的名稱。

叫用作業系統的功能程式，是 Forth 向外發展時不得不採用的性能擴充技術。在 DOS 的時代，是透過呼叫 INT xx 中斷程式的方式來完成此事， Windows’ 以後的作業系統，則改成了 call 呼叫編號的方式來完成。而在 XP 以前，參數都放在堆疊上傳遞，在 W7 以後，參數就改成放在 CPU 中的暫存器傳遞。在 Windows’ 中的放置秩序與在 Linux 中的放置秩序有所不同，使用時必須注意。詳細的細節，可以透過審視 Forth 系統中 call 指令的源程式了解規矩， Forth 的使用者就能看到這些可貴的資源，但是，被呼叫的對方，只願單純使用 Forth 的人就很難繼續深入了，被 call 那邊的程式資源，不是百分之百透通的。

以上是靜態連結程式的基本叫用法。 隨著作業系統的不同，另外發展出了一種直接透過文字視窗大肆叫用作業系統指令的方式完成的呼叫功能，使用一個稱為 system 的指令可以達到目的，這些都是屬於靜態連結技巧的使用方法。關於動態連結程式的叫用方法，也有一般使用上的通用協定，範例程式涉及時才討論。 :

\ HeightWidth.f
\ 20140210
\ 資料源自 Win32Forth 論壇

: HeightWidth ( -- )
SM_CYSCREEN call GetSystemMetrics   \ 768  ok (y for the height)
cr ." The height of screen is " .
SM_CXSCREEN call GetSystemMetrics   \ 1366  ok (x for the weight)
cr ." The wide of screen is " .
;

: main HeightWidth ;

cr cr
.( Usage: ) cr
.( main ) cr