一百個例題 (21~24)

一百個例題 (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.   
  ;