2024年9月30日 星期一

一百個例題 (76 ~ 80)

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

張貼者： 於

沒有留言:

張貼留言

訂閱： 張貼留言 (Atom)