一百個例題 (96 ~ 100)
曾慶潭 Ching-Tang Tseng
ilikeforth@gmail.com
Hamilton, New Zealand
7 October 2024
第(96)個範例,是一個只能在 Lina64 中執行,可以將網址內原為中文格式碼的中文文字顯示出來的程式。
由於在 W10 環境中,文字視窗內中文文字的顯示方式與 Ubuntu 作業系統中者不同,所以,不要在 W10 環境中使用這個程式。
在 Win32Forth 續用機會逐漸接近尾聲的時候,為它重新開發這個主題的程式,意義不大,此處只能盡量解釋程式的精神所在。
程式後面附貼了 20241007 在 Ubuntu V20.04 啟動 Lina64 實際載入並執行這個程式後的顯示結果,讓大家體會程式的實際用處。
我每天使用電腦來做許多事情,瀏覽網路時,收集無數的訊息,影、音、圖、字通通都有。除了訊息,為了便於回查資料出處,一定要把網路的地址也同時留下來。又由於國際各國文字編碼格式的不同,如果網址中使用了非英文的文字,每當操作滑鼠進行反白、複製取得網址,再到文字檔案中去貼上,儲存這個網址時,就會出現所謂 UTF-8 編碼與 BIG-5 編碼之不同的問題,而令取得的網址出現很多以 『 % 』 號前引,後面跟著十六進制碼的東西,無法讓人直接看懂。而且,在 W10 環境擷取這種網址存進檔案,再轉往 Linux 環境打開檔案來看時,就只能顯示這類編碼,想看正確的結果,都得經過字型編碼的規格轉換。當然,系統中有可以自動把整個檔案內容翻譯出來的指令可用,翻譯完就不能再回頭給原作業系統使用了。我很少這樣做,每次,也只有單一個網址的需求問題。
這樣擷取的網址資料,經過 20 幾年之後,留下來的有用網址,已經是成千上萬個了。晚期的網址,才會有這樣的問題,早期的,都是只有英文。現在,每當想要回查網址訊息時,就常常會遭遇到因尾綴 UTF-8 字型編碼,而無法直接看出網頁實際名稱的狀況。
我用 Forth 寫程式,處理所有我會遇到的問題,上述狀況,也是在改用 64 位元之 Linux OS 環境後,才出現的困擾,所以就設計了這個程式。實際上,在 32 位元的 Win32Forth 系統中,也能設計得出具有同樣功能的程式,程式的難度,只屬於二級的程度而已。處理關鍵,只在設法進行刪除所有的 % 符號,轉換出來的結果,在純文字視窗上,就能顯示出原來的中文字了。 % 符號,等於只是一個前引的格式碼,沒有其他的意義。
關於字串處理方面的程式,我在百例中的介紹並不很多,原因就是它們都不屬於數學計算的範疇。我在設計應用程式時,仍然常需面對字串處裡問題。可是, Forth 標準指令中,並沒有許多現成的字串處理指令可用,例如字串的存取、裁剪、比較、增刪 ..... 等等,都沒有現成指令可用,必須自己設計。我並不喜歡討論這方面的技術,因為換個系統之後,一切都得重來,為了適應新環境,重新設計程式都很耗時間,除非是為了非做不可的事情,我才會去仔細研究。
我在 20181202 貼出過關於字串處裡的基礎教學文章:字串。隨後,國際論壇上大概也因此而大肆討論了許久,最後,好像就如我所說的,許多指令雖很有用,但公認不必列為標準,會用的人,兩三行就能自行設計得出來,不熟 Forth 者仍然不會。自此以後,也不再有熱心的人士,願意像我這樣,大耗篇幅寫例題程式了。
\ (96) url translator 20180710
0 value ptr
create url 32 cells allot
create web 32 cells allot
: reset$buf
url 32 cells erase
web 32 cells erase
cr ." Input url string $ " cr
input$ url $! ;
: setptr
url $@ drop to ptr ;
: url>web
begin
ptr c@ 0 <>
while
ptr c@ [char] % =
if HEX 0. ptr 1+ 2 >number 2drop d>s web $c+ DECIMAL 3 +to ptr
else ptr c@ web $c+ 1 +to ptr
then
repeat ;
: show
cr ." url = " url $@ type cr
cr ." web = " web $@ type cr ;
: main
reset$buf
setptr
url>web
show ;
main
\s
ching@ctt:~$ ./l
AMDX86 ciforth 5.3.0
fload 96.f
Input url string $
http://www.teepr.com/376386/lexychang/%E4%B8%B9%E9%BA%A5%E5%A4%A9%E4%BD%BF%E7%9A%84%E7%B4%A0%E9%A1%8F%E7%85%A7%E6%9C%83%E8%AE%93%E4%BD%A0%E7%9F%A5%E9%81%93%E4%BB%80%E9%BA%BC%E6%98%AF%E3%80%8C%E7%9C%9F%E6%AD%A3%E7%9A%84%E4%BB%99%E5%A5%B3/
url = http://www.teepr.com/376386/lexychang/%E4%B8%B9%E9%BA%A5%E5%A4%A9%E4%BD%BF%E7%9A%84%E7%B4%A0%E9%A1%8F%E7%85%A7%E6%9C%83%E8%AE%93%E4%BD%A0%E7%9F%A5%E9%81%93%E4%BB%80%E9%BA%BC%E6%98%AF%E3%80%8C%E7%9C%9F%E6%AD%A3%E7%9A%84%E4%BB%99%E5%A5%B3/
web = http://www.teepr.com/376386/lexychang/丹麥天使的素顏照會讓你知道什麼是「真正的仙女/
OK
第(97)個範例程式是用來測試 Lina64 系統中的矩陣及亂數函數應用的結果,兩項成果顯示都有很好的性能。
關於這兩大項的測試,我分開來解釋,並附上實際執行後的測試輸出結果,便於說明。
關於矩陣的測試,只偏重於測出帶 0 指標的單元,是否已被正常的納入系統?用不用都無所謂,輸出時就必須顯示它們的存在。這一特點,在原本的 Win32Forth ABC Forth 系統中沒有刻意安排。我是在使用許多大型數學計算套裝軟體後,才體會出必須令其存在的,否則會有很多現成的程式,在抄用時會遭遇很大的困擾,必須人工自行調整指標來對應,不如就在建系統時加入此一功能,會方便許多。當然,每次宣告陣列、矩陣來用時,就會損失一些可能用不上的記憶體。今日世界,這已不再是問題,那就不要省,對我而言,加設計也只是舉手之勞。這是這個範例測試出設計成功的一個項目。宣告後,沒用到的具有 0 指標之單元,內容全為 0 ,這就是成功的結果。
其次,關於亂數,分成三大項來看結果,我設計系統時沒有前人的文獻可以參考,自己設計後,使用上的『協定』也只好自己規定了。
所有的亂數都是炒作整數亂數後延生出來的,炒作方法就是把位元花樣向左、向右移動幾個質數次,再與原花樣進行 XOR 運算,保證與原花樣完全不一樣,這樣搞個幾次,就成了新的亂數,下一個照樣畫葫蘆。起始花樣,則藉 CPU 旁邊會跑個不停的計時器內容之下半段位元花樣來亂化,使用這種軟、硬通吃的技術,效果就會很好。
這樣產生的整數亂數,因用到全部的位元來生亂,最高位元當然也被炒到了,程式才容易設計,所以,直接使用時,會有負值出現,加一個組合語言指令 ABS 就能解決問題,亂數產生的速度非常快,是 giga 分之幾秒就能完成。但我認為負的亂數也很有用,功能就保留下來了,用時自己注意。
我的浮點數是全用整數四則運算完成的十進制浮點,拿整數亂數直接用,當然就更方便。只是用時要注意,一般的基本用法,都是先取得只介於 0 至 1 之間的浮點亂數,想放大多少,自己乘上去後才使用,就得了。
我的複數系統是一般數學運算套裝軟體系統上不太具有的東西,所以就不用找規格了,使用『協定』也只好自己創立,這個範例中有簡單的使用規格,叫用亂數函數時,使用者提供的參數就是新產生之亂數的上限範圍,實數的部分可以與虛數的部分不一樣,例如此範例中,實數上限是 100 ,虛數上限是 10 ,測試結果,顯示都很好,非常恰當。
\ (97 random number test
\ beware! In integer mode, need a NEGATE operation to get negative random number
: NTEST 10 0 DO CR 10 0 DO 100 RAND negate 4 .R LOOP LOOP ;
: FTEST 10 0 DO CR -1.0 E 2 FRAND F. LOOP ;
\ beware! upper bound input format ignore + and ii
: ZTEST 10 0 DO CR -1.0 E 2 1.0 E 1 ZRAND Z. LOOP ;
3 integers i j k
complex z0
3 3 3 [tensor] zzz
: init
basic
10 let z{ z0 = 0.0 E 0 + 0.0 E 0 ii }z
20 for i = 1 to 3
30 for j = 1 to 3
40 for k = 1 to 3
\ attention! this is the usage to create a complex random number
50 let z{ zzz ( i j k ) = RND ( 1.0 E 2 - 1.0 E 1 ii ) }z
60 next k
70 next j
80 next i
90 end ;
: MAIN
randomize
ntest cr ftest cr ztest cr
init
basic
10 for i = 0 to 3
20 for j = 0 to 3
30 for k = 0 to 3
40 print " zzz( " ; i ; j ; k ; " )= " ; z{ zzz ( i j k ) }z
50 let z{ z0 = z0 + zzz ( i j k ) }z
60 next k
70 next j
80 next i
90 print cr cr " z0 = " ; z{ z0 }z
100 end
cr ." --- The end of testing --- "
;
\s
ching@center:~$ cd lina64
ching@center:~/lina64$ ./f5124
AMDX86 ciforth 5.3.0
fload 97
I : ISN'T UNIQUE
J : ISN'T UNIQUE
OK
main
-10 -8 -72 -61 -31 -86 -14 -7 -86 -20
-42 -94 -92 -27 -86 -45 -57 -77 -34 -4
-98 -89 -61 -28 -86 -66 -36 -71 -68 -71
-35 -10 -10 -84 -85 -82 -23 -14 -89 -20
-34 -51 -10 -87 -58 -67 -21 -64 -43 -51
-26 -18 -63 -99 -30 -22 -34 -92 -25 -78
-58 -58 -98 -64 -93 -54 -23 -43 -99 -38
-83 -68 -50 -68 -11 -74 -44 -49 -85 -42
-52 -18 -23 -93 -39 -41 -37 -16 -49 -67
-40 -39 -66 -99 -70 -45 -99 -16 -17 -28
-1.913628913
-89.96769746
-90.80277685
-93.25782136
-12.5438047
-68.49962742
-8.071056155
-61.09243441
-47.21671679
-32.33685587
-3.989974468+ 4.9740950527ii
-77.58595045+ 0.7664840648ii
-42.15458337+ 0.0229446231ii
-5.083676585+ 6.5485581007ii
-62.74442857+ 5.2724854219ii
-70.1131553+ 6.2933488809ii
-37.53475328+ 6.8368182786ii
-38.22238984+ 1.2874766663ii
-96.85347849+ 4.7545419669ii
-3.281314505+ 5.6388045221ii
zzz( 0 0 0 )= 0
zzz( 0 0 1 )= 0
zzz( 0 0 2 )= 0
zzz( 0 0 3 )= 0
zzz( 0 1 0 )= 0
zzz( 0 1 1 )= 0
zzz( 0 1 2 )= 0
zzz( 0 1 3 )= 0
zzz( 0 2 0 )= 0
zzz( 0 2 1 )= 0
zzz( 0 2 2 )= 0
zzz( 0 2 3 )= 0
zzz( 0 3 0 )= 0
zzz( 0 3 1 )= 0
zzz( 0 3 2 )= 0
zzz( 0 3 3 )= 0
zzz( 1 0 0 )= 0
zzz( 1 0 1 )= 0
zzz( 1 0 2 )= 0
zzz( 1 0 3 )= 0
zzz( 1 1 0 )= 0
zzz( 1 1 1 )= 16.138346694- 1.6922562207ii
zzz( 1 1 2 )= 11.673369547- 0.0886917048ii
zzz( 1 1 3 )= 12.744789187- 3.7724680765ii
zzz( 1 2 0 )= 0
zzz( 1 2 1 )= 31.962691299- 1.3816571182ii
zzz( 1 2 2 )= 40.44457283- 1.5831628868ii
zzz( 1 2 3 )= 20.592126444- 4.1070317085ii
zzz( 1 3 0 )= 0
zzz( 1 3 1 )= 5.7211277157- 1.9009267784ii
zzz( 1 3 2 )= 79.010045666- 1.7076510055ii
zzz( 1 3 3 )= 36.175071964- 1.1219479774ii
zzz( 2 0 0 )= 0
zzz( 2 0 1 )= 0
zzz( 2 0 2 )= 0
zzz( 2 0 3 )= 0
zzz( 2 1 0 )= 0
zzz( 2 1 1 )= 15.89740986- 2.4554373865ii
zzz( 2 1 2 )= 82.495140485- 8.5362387692ii
zzz( 2 1 3 )= 17.291546835- 8.3884016077ii
zzz( 2 2 0 )= 0
zzz( 2 2 1 )= 15.663485123- 8.9183735609ii
zzz( 2 2 2 )= 87.857473069- 2.1572277506ii
zzz( 2 2 3 )= 8.0974462118- 4.8481376385ii
zzz( 2 3 0 )= 0
zzz( 2 3 1 )= 16.89172751- 5.8226671362ii
zzz( 2 3 2 )= 67.783749813- 1.849634518ii
zzz( 2 3 3 )= 57.98394664- 6.3832837116ii
zzz( 3 0 0 )= 0
zzz( 3 0 1 )= 0
zzz( 3 0 2 )= 0
zzz( 3 0 3 )= 0
zzz( 3 1 0 )= 0
zzz( 3 1 1 )= 5.9138480523- 6.632469444ii
zzz( 3 1 2 )= 92.038071469- 2.2556417577ii
zzz( 3 1 3 )= 18.729357865- 7.4502652215ii
zzz( 3 2 0 )= 0
zzz( 3 2 1 )= 67.299863245- 2.0522561644ii
zzz( 3 2 2 )= 80.622464446- 2.3218369173ii
zzz( 3 2 3 )= 74.024361093- 3.2010649041ii
zzz( 3 3 0 )= 0
zzz( 3 3 1 )= 95.399681736- 5.8407040896ii
zzz( 3 3 2 )= 59.354835252- 6.4971174443ii
zzz( 3 3 3 )= 12.681339537- 9.1802079012ii
z0 = 1130.4878896- 112.1467594ii
--- The end of testing --- OK
第(98)個範例程式是用來解決一個根據巴斯卡(Pascal)三角塔原理為基礎所演變出來的問題,讓程式可以進行矩陣運算,也讓結果易於被核驗。
高階二項式 (a+b)**n 的展開式,每一項的係數所形成的固定變化,曾被一位法國的數學家叫布萊斯‧巴斯卡(Blaise Pascal,1623年6月19日-1662年8月19日),妥善整理出來,是一個漂亮的三角塔式的係數數列,故得其名。但實際上,對這個問題有過研究的古老文獻,非常的多,顯然,發現人並不是巴斯卡,甚至於中國人也有,這方面,就不用去考究了,我們直稱其為巴斯卡三角塔便可。
巴斯卡三角塔並不容易被電腦程式漂亮的印出來,為了方便,只適合把三角塔印成上三角矩陣,或者是下三角矩陣。因此,這個程式的題目,就要求計算出把這兩種矩陣加起來後的結果到底為何?
我控留此程式的目的,主要是因為這樣的題目能磨練程式的寫法,另外,也是個能夠用來驗證矩陣的資料結構設計得是否正確之題目。算是一個優良的測試材料。關於磨練程式設計技術的部份,請自行看看程式本身,就能了解。關於驗證資料結構是否正確的問題,則還能進一步探討。
巴斯卡三角塔的整齊結構,易於只用肉眼觀察,就能看出對錯。兩個矩陣相加後,因為元素內容的花樣,係與對角線對稱,也易於只用肉眼觀察,就能看出對錯。一般而言,我們使用的電腦,宣告出 20X20 的矩陣,螢幕上尚能印出不需跳列的輸出,能測試到這樣的尺度,大約也足夠表示設計無誤了。
所謂的進一步探討,並不絕對需要,但對矩陣運算系統的設計與發展而言,就很有必要。大家可以看看程式的寫法,您能發現,就算要把矩陣印出來,也得在被印出程式之指令外面,先行宣告出矩陣的維度,才辦得到。後續那個矩陣的加法運算,亦然。這就表示,寫成的這套程式,代表著所有的矩陣運算程式都得這樣處理。
我在發展出許多數學計算系統後,自覺能夠改善矩陣體系運算系統的性能,能夠設計得出 M+ , M- , M* , M/ , Mx , M. , 1/M , Mt , M0 , Mneg ..... 等等,可以不管宣告維度是多少,就能直接使用的運算指令。這樣的性能,必須依靠矩陣在宣告產生時,就把維度也存在資料結構中才行,我辦到了這項設計,但沒有進行矩陣運算系統的發展,因我個人很少有接觸這方面問題的機會。
事實上,我還有許多的構想,可以設計出更多的數學運算體系,都很有意義。但自從進入 64 位元的環境後,面對的系統發展問題,不是更多的數學體系,而是更多的不同 Forth 基本平台。因此,發現研究的主要方向,自己必須嚴格把握住,才有發展的坦途可走。整個系統的基礎,若不專心設計得更好,發展再多的數學體系,只會弄亂自己的陣腳。於是,停止了矩陣體系運算系統的開發了。
Win32Forth 系統,可以順利的執行出這個程式。純粹只處理整數的數學問題,這類程式都能這樣被直接執行得出來。
\ (98) Pascal matrix generator
4 integers i j k n
20 20 (matrix) mm
: i-1 i 1- ;
: j-1 j 1- ;
: printmm
basic
100 for i = 1 to n
110 for j = 1 to n
120 let k = mm ( i j )
130 run k 8 .r space
140 next j
150 run cr
160 next i
200 end ;
: upm
basic
10 for j = 1 to n
20 let mm ( 1 j ) = 1
30 next j
40 for i = 2 to n
50 let mm ( i 1 ) = 0
60 for j = 2 to n
70 let mm ( i j ) = mm ( i-1 j-1 ) + mm ( i j-1 )
80 next j
90 next i
500 end ;
: lpm
basic
10 for i = 1 to n
20 let mm ( i 1 ) = 1
30 next i
40 for j = 2 to n
50 let mm ( 1 j ) = 0
60 for i = 2 to n
70 let mm ( i j ) = mm ( i-1 j-1 ) + mm ( i-1 j )
80 next i
90 next j
500 end ;
: spm
basic
10 for i = 1 to n
20 let mm ( i 1 ) = 1 :: mm ( 1 i ) = 1
30 next i
40 for j = 2 to n
50 for i = 2 to n
60 let mm ( i j ) = mm ( i j-1 ) + mm ( i-1 j )
70 next i
80 next j
90 end ;
: pascal ( n -- )
[[ n ]] !
basic
10 run cr upm
20 run cr printmm
30 run cr lpm
40 run cr printmm
50 run cr spm
60 run cr printmm
70 end ;
page 5 pascal
第(99)個範例是在『Euler數學問題網頁』中,條列之一個問題的電腦程式解法。
問題問:直角三角形的邊長為 1000 ,請問三邊邊長均為整數時的長度各為多少?
Euler’s 數學問題網頁,是一個最近兩年才聞名起來的網頁,所出問題,大部分是需要能夠運算大數目字的系統才能解答的題目,現已大約收集了六百多個問題,貼出在網上,讓人隨意發揮。剛開始時,網頁出入管制,沒有那麼嚴格,我有空時,就會去拿題目,試一試自己的系統,那時,大約我都能解,解題時,只覺得必須全靠電腦,很難有單憑思考就能解析出來的題目。後來,該網頁開始進行管制了,必須登記身份才能使用,我就放棄了追蹤。
話雖如此,那種數學題目,反而代表了現今網路環境熱門話題的演變趨勢,都屬於同一類型的問題。例如,要您從海量資料中,比對出合乎指定條件的網上訊息,您就得像處理無限位數的大數那樣,把一小節數字用電腦找出來。所以,這個數學問題網頁的建立,是一個時代的產物。以後,如果您有興趣也有本領解這種題目,就可以根據名稱,自行前往該網頁進行探索。我不參與,所以也沒留網址,但很容易搜索得到。
我留下這個程式的目的,除了題目優良外,還有個簡單意義,就是畢氏定理的運用,能讓我寫成的程式,比別種程式語言的語法更為清楚、簡單。本來,問題的表面只有ㄧ個條件,就是 a+b+c=1000 ,其中,有兩個獨立變數,一個因變數,條件只有一個時,沒辦法有特解,只能有也附帶條件式的通解。如果學過了數學上的畢氏定理,您就知道另有一個很好的輔助條件可以利用,也就是斜邊的平方必須等於另兩個直角邊之各自平方的和,亦即 a**2+b**2=c**2 。
這麼一來,寫程式時,語法好的系統,就很容易讓人看懂,執行起來,就不容易出差錯, BASIC 程式語言的語法,就具有這種優點。仔細看這不到十列的程式,使用兩個環路、一個邏輯判斷來完成,兩個環路就是給兩個獨立變數用的,範圍就指定在 1 到 1000 之間。邏輯判斷就只給畢氏定理使用,以取得因變數之值,合條件者就是答案,不合條件者就不是答案。這樣的電腦解法,讓普通人純用手算,就會成為不切實際的解題方法,用電腦算,又快又準,立刻就能得到答案。此程式,只處理整數,就可以在我設計的系統中執行得出來,我就不在此處提供解答,印出來的答案一定對,您可以自行核算。
本來, Forth 系統的本質,就是一種具有延伸性的程式語言,我沒附加 BASIC 式語法的執行能力前,早就有名家為 Win32Forth 系統添加了組合語言的功能,只是後來大家都不太使用了,所以罕有人再強調。但是,這套 Win32Forth 中的 Assembler 非常有名,許多專門介紹組合語言的網頁中,竟然也把這套組合語言列為一套正式的特產,性能既好也完整。
事實上,我自創添加進去的 BASIC 式語法,等同於 Forth 系統中的 Assembler ,就是 Forth 中的 BASIC 。
一旦一套 Forth 系統,同時具有以他種語言之語法執行出結果的功能時,這個 Forth 系統,就也具有能夠在語法中互相來去、交流穿插的功能。這也就是為什麼我能利用這種精神,在 BASIC 語法中創建出一個 RUN 指令的原因。我把雙方語言的特質都搞清楚了,賦予這種好處,就不是難事。這樣的發揮, Forth 系統的發明人並沒有善加利用。實際上,他有點痛恨 BASIC 在當年霸佔市場的強勢狀態,那時,他好像是只為了要藐視一下 BASIC ,才發表了他的那篇 Tiny BASIC Compiler 論文。我,則是取長補短,設計出了自己的東西。
另外, Charles H. Moore 終生都痛恨浮點數運算,所以他也從來不寫浮點數運算方面的程式。這一點,我個人認為不對,所以,我自創系統中所有的浮點數計算功能,如何實現?都是我自己安排、自己實現的設計。當然,仍屬於 Forth 浮點範疇的部份,只要能師法於現成的成果,我就取用,如果沒有,我就自行設計,現在也已很熟了。
我仍然是站在巨人的肩膀上創作系統,自己的成就只有一小部份。
(99)euler9
\ Subject: Project Euler Problem 9 : unique Pythagorean triple summing to 1000
\ a+b+c=1000
3 integers a b c
: main
basic
10 for a = 1 to 1000
20 for b = a to 1000
30 let c = 1000 - a - b
40 if ( a * a + b * b ) = ( c * c ) then 60
50 goto 70
60 print a , b , c
70 next b
80 next a
90 end
;
第(100)個範例,是一個能夠操作出讓單一個位元為 1 或為 0 的工具程式。
整台電腦中任一位址的內容,可以被取來放入 bitbuffer 變數內,然後才用這個程式中的操作指令,來令其中的單一個位元成為 1 或成為 0 。
操作的方法,在程式中有舉例,也只列出了一個顯示結果,要怎樣用才會更為方便?請自己加添設計出更簡單的指令,就能隨您所願。原本這個程式是在 64 位元環境中由 Bob Dickow 設計的,我試用後覺得比我自己設計的更好用,程式略作修改後,就被收集進這個範例了。
這種『位元操作』能力的展現,是一種 Forth 程式語言之特殊性能,別種程式語言都很難辦到這麼方便的設計格式,我用了 40 多年的 Forth ,只要想到要進行一些位元操作的工作時,深知非 Forth 莫屬。
我設計過一份以 Wina32 為平台,全由叫用硬體 co-processor 方式來完成所有浮點計算工作的系統。那個系統是標準 32 位元的設計,數據在 CPU 與 co-processor 間傳遞,執行所有的軟體指令都沒問題,硬體接線也與數據傳輸要求規格相符,所以,我很快便完成了設計,實際使用時,也很方便,很快就完成了 ABC Forth 的增建。
後來,我轉往 64 位元以 Lina64 為平台重新發展系統時,首先遭遇到的問題,就是從 CPU 內傳送 64 位元的單整數當 32 元的雙整數到 co-processor 後,因為硬體接線的格式與 32 位元者完全不同,我就不能再用 co-processor 設計浮點系統了。這問題不單純只是軟體指令的不能使用而已,硬體接線與數據傳輸時脈規格也完全不同,所以不能再寫組合語言實現同樣的夢想。除非能令系統回到 32 位元執行浮點運算,算完,再回到 64 位元系統顯示結果,這很麻煩,叫用 C 的庫存程式可以完成,但 gForth 系統就顯示了執行速度非常慢的問題,不值得這樣做。
我為什麼知道上述問題?就是因爲我使用了 Forth 的位元操作能力,仔細的測試過問題。把簡單位元花樣的數字,在 CPU 與 co-processor 間來回傳輸、簡單執行,再印出前後結果,一經比對,發現連一個簡單位元花樣的 1.2E0 浮點數,都沒辦法處裡得正確。無論以軟體設計了多少個委屈程式,都沒辦法獲得正確的結果。最後,只好放棄使用硬體,改成走自創純軟體設計之十進制浮點系統的發展方向,現已完成了設計。
這百例中,有許多範例都是我在發展系統期間所使用的工具程式,都有著深刻的印象與使用經驗,這些範例都不是幾天之內就能完成的研究,它們能用時,都是經過長久研發與各種測試步驟後,才得到的結果。
我貼完了百例,也對自己過去所做過的工作,進行了一次完整的複習,對我自己而言,很有益處。
\ (100)FORTH BIT Arrays , by Bob Dickow dickow@turbonet.com
\ version 1.0.3 Mar 2, 2018
\ Notes:
\ endian agnostic, address unit agnostic (untested)
\ Vocabulary (internal dependencies) ^2MASKS bmask _BIT@
\ Vocabulary: 2^ BIT@ BIT! BOOL@ BOOL!
\ BOOL@ and BOOL! are convenience words to avoid converting FALSE to 1
\ 2^ may be handy outside the array access words. Included for compatibility with earlier versions
\ Buffers for a bit array should be a minimum of 1 CELL in size
DECIMAL
CELL 8 * CONSTANT AddressUnit \ 8 is bits in a byte
\ calulates the power of 2 ( n -- 2^n )
: 2^ ( n -- 2^n )
1 swap lshift
;
\ compiling word, creates an array of CELL * 8 powers of 2 for masking use in the array fetch and store words.
\ does: ( nindex -- nvalue ) \\ nindex is 0, 1, 2 ... giving 2^31, 2^30, 2^29 ... 2^0
: ^2MASKS ( -- )
CREATE
0 AddressUnit 1- DO
i 2^
, -1 +LOOP
DOES> SWAP CELLS + @
;
^2MASKS bmask \ lookup table of powers of 2 for the BIT@ and BOOL@ words
\ factored out, internal use by BIT@, BOOL@
: _BIT@ ( caddr n-index -- n )
AddressUnit /MOD ( caddr n n )
CELLS ( caddr n n )
ROT + ( n caddr )
@ SWAP ( nval n )
bmask AND ( n )
;
\ given a memory buffer address caddr and index n representing at bit-wise
\ index, return the value of the bit at that index as 1 | 0.
: BIT@ ( caddr n-index -- n )
_BIT@ ( n )
0= NOT
ABS
;
\ store a value TRUE|FALSE or 1|0, given the value (TRUE|FALSE or 1|0), an address, and a bit-wise index
\ (synonym for BOOL! below)
: BIT! ( n-value caddr n-index -- )
ROT ABS 1 AND ( caddr n n ) \ bounds check
IF ['] OR ELSE ['] XOR THEN >R \ set masking op
AddressUnit /MOD ( caddr n n )
CELLS ( caddr n n )
ROT + DUP ( n caddr caddr )
@ ROT ( caddr caddr nval )
bmask R> EXECUTE ( caddr n )
SWAP !
;
\ given a memory buffer address caddr and index n representing at bit-wise
\ offset, return or store the value of the bit as TRUE|FALSE values.
: BOOL@ ( caddr n-index -- flg )
_BIT@ ( n )
0= NOT ( n )
;
\ Synonym for BIT!
\ store a value TRUE|FALSE ( or 1|0 ), given the value n, a base address of the array, and a bit-wise index into it.
: BOOL! ( n caddr n-index -- )
BIT!
;
DECIMAL
\ create the buffer:
CREATE bitbuffer 2 CELLS ALLOT
\ initialize buffer to all on bits:
\ 4294967295 DUP bitbuffer DUP ROT ROT ! CELL+ !
-1 DUP bitbuffer DUP ROT ROT ! CELL+ !
\ store a bit at position 34
0 bitbuffer 34 BIT!
0 bitbuffer 67 bit!
0 bitbuffer 127 bit!
\ for lina64
: bottomruler130 cr
." 01234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890 " cr
." 0=========1=========2=========3=========4=========5=========6=========7=========8=========9=========0=========1=========2=========3 " cr ;
: binary 2 base ! ;
: (du.) <# #s #> ;
: du. (du.) type bottomruler130 ;
\ test it:
bitbuffer 2@ binary du. decimal \ displays: 1111111111111111111111111111111111011111111111111111111111111111
一百個例題 (91 ~95)
Ching-Tang Tseng
Hamilton, New Zealand
6 October 2024
第(91)個範例是專門探討以二分法(bisection)求解方程式實數根的技術。
請注意前面的第(76)個範例,也是利用二分法技術,用來算出浮點數的平方根。
這個方法很好用,也絕對可靠,只要對此技術認識得夠清楚,就能用來解決所有方程式求根的問題。
我個人認為,二分法是當今全世界從事於電腦數學運算應用者都應該熟悉的技術。
這個範例的全套英文解說,曾經貼示於我的個人網頁,是 20180416 所刊出的文章: Bisection Method To Solve Equation 。(該文因網頁改版,原貼內容已遭破壞,我會找時間進行修改。)
範例所解的題目,是經過刻意挑選的,您無法手解出它所有的三個實根。題目中的方程式,是代表某種化學反應趨勢的實測數據描述出來的表示式。當初,我大女兒念生化系時,拿這個方程式回家來問我,如何才能把正的兩個實根給算出來?那兩個根,要用來定出激烈化學反應的界限。於是,我發展出了這套程式。
我用英文貼文,大部分台灣人都看不懂,教大家看英文資料,不是我的責任,我貼文的目的,是讓全世界知道我的系統能解這種問題。很多人看得懂英文,但也只看出了表面意義。唯一看過此文,並在國際論壇回應網貼意見者,只有 ciForth 的作者,他才意識到系統能解此題的一個很重要的特點,他貼文表示此特點必須為後世所師法。
我若不刻意指出這個特點,大家可能也很難注意到,但如果從函數圖來看問題,就能看出問題所在了。函數圖中在 x = 80 的位置,是一個奇點,它的函數值是正的無限大,過去,所有的數學計算系統,在演算過程中,若碰到這種計算,會立刻執行例外處理、發出警告、停止執行。可是,我的系統就故意從這一個奇點開始計算,還進行了兩次,系統碰到了,照樣繼續算,直到產生結果,不需要發出任何警告。我的系統可以接受無限大的數字,而且,這個數字恆就是個無限大,使用者若想印出來看,系統就印出是帶有 + 或 - 符號的 infinit 。
二分法的基本精神,前文已經介紹過,這個範例則專挑學校老師也無法手解的問題來示範二分法的功能,它的發揮,主要是必須先靠繪出函數圖來輔助,您才能確認函數的實根大約會落在那裡?現行電腦的功能,已經能夠非常快速輕鬆的解決函數繪圖問題,求解實根時所需提供給程式使用的上下限二值,根據圖形就很容易取得,剩下來的冗長計算,全部交給電腦。
我在發展到 Lina64 系統之後,有感於這種可用的寶貴程式值得庫存。因此,開始把程式寫得非常整齊,並添加了許多必要的註解,讓自己在日後想要再度使用時,一目了然:
新給的方程式定義,就只需修改 FX 的定義內容,而且,寫法就是日常書寫的格式。
新給的上下限二個值,就只需修改 root1 , root2 , root3 ..... 等等就夠了。
仍然只需執行 main 指令,就能得到所有的答案。
我在 Lina64 中也發展了語音輸出及單音樂譜的執行能力,所以在貼文時,故意展示出這兩項功能。
語音播出是『莫斯科廣播電台的台號音樂與呼叫』,人聲與樂聲均能使用程式輸出。
繪圖的輸出,是藉 Forth 系統控制 Python 程式去操控 GnuPlot ,是個論秒就能繪出函數圖的設計,那個環境,可以直接操作圖形歸檔。
所有執行程序均在 Lina64 中一氣呵成。
浮點運算部分與 ABC Forth 的功能,是我個人的創作,此系統未曾釋出,也不想釋出,因為大家可以注意看所有的浮點數輸入方式,跟此前全世界的任何系統都不一樣,假數部分與指數部分是被中間插入一個 e 或 E 來完成的,執行速度很快,因為數字都不需要轉換,就能直接算出浮點運算的結果。我不想擾亂世界的規矩,自己使用就無所謂,所以不釋出系統。
系統的所有功能, Lina64 的原作者大約也很難完成,所以他仔細讀過我的貼文後刻意回應。他說,他就是希望他的系統能夠辦到這些事情。他講的一點也沒錯,能實現這些功能的基本要素,是他的系統內已經設計好了跟作業系統所需要的介面工具指令,有好幾個,各個都能正確執行出對的結果。只是這位作者可能沒有像我一樣去普讀許多系統資料,供自己套進 Lina64 使用。這位作者的出身是高能物理的博士,我跟他一樣,都懂核子科學。我未曾與他聯繫過,只在國際論壇上神交,但我很感激他公開所有設計的系統。
:
\ Bisection method to find real roots of equation in Lina64 ABC forth
\ Equation: y=ln|140/(80-x)|-((60+x)/53.0516477)
\ Guessing range: root1[-100, -50], root2[50, 80], root3[80, 100].
\ ===============================================
\ (1). Variables declare heap
variable Counter
integer Steps
10 reals X Y XA XB XP FA FB FP EPS TEMP
\ ===============================================
\ (2). Input heap
\ Equation is defined here.
: FX
{{ Y = LN ( ABS ( 140.0 e 0 / ( 80.0 e 0 - X ) ) ) - ( ( 60.0 e 0 + X ) / 53.0516477 e 0 ) }} ;
\ guessing range of roots are defined here.
: root1
{{ XA = -100.0 e 0 }} {{ XB = -50.0 e 0 }} ;
: root2
{{ XA = 50.0 e 0 }} {{ XB = 80.0 e 0 }} ;
: root3
{{ XA = 80.0 e 0 }} {{ XB = 100.0 e 0 }} ;
\ All parts above are changeable setting.
\ ===============================================
\ ===============================================
\ (3). Main code heap
\ All parts below are fixed.
: Fbisection
0 COUNTER !
BEGIN
{{ X = XA }} FX {{ FA = Y }}
{{ X = XB }} FX {{ FB = Y }}
{{ EPS = ABS ( XA - XB ) - 1.0 E -12 }}
EPS F0>
WHILE
{{ XP = ( XA + XB ) / 2.0 e 0 }}
{{ X = XP }} FX {{ FP = Y }}
{{ TEMP = FA * FP }}
TEMP F0>
IF {{ XA = XP }} ELSE {{ XB = XP }} THEN
1 COUNTER +!
REPEAT
CR ." Counts = " COUNTER @ .
{{ TEMP = ( XA + XB ) / 2.0 e 0 }}
CR ." Solution = " TEMP F. cr
;
: Bbisection
BASIC
10 let steps = 0
20 LET { X = XA }
30 RUN FX
40 LET { FA = Y }
50 LET { X = XB }
60 RUN FX
70 LET { FB = Y }
80 LET { EPS = ABS ( XA - XB ) - 1.0 E -12 }
90 IF { EPS <= 0.0 e 0 } THEN 210
100 LET { XP = ( XA + XB ) / 2.0 e 0 }
110 LET { X = XP }
120 RUN FX
130 LET { FP = Y }
140 LET { TEMP = FA * FP }
150 IF { TEMP > 0.0 e 0 } THEN 180
160 LET { XB = XP }
170 GOTO 190
180 LET { XA = XP }
190 LET Steps = Steps + 1
200 GOTO -20
210 PRINT " Steps = " ; steps
220 LET { TEMP = ( XA + XB ) / 2.0 e 0 }
230 run cr ." Solution = " TEMP f. cr
240 END
;
\ ===============================================
\ ===============================================
\ (4). Manipulation heap
: mainF
BASIC
10 run root1
20 run Fbisection
30 run root2
40 run Fbisection
50 run root3
60 run Fbisection
100 end
;
: mainB
root1 Bbisection
root2 Bbisection
root3 Bbisection ;
\ Russia: Moscow Radio Station signed voice
: Signal
MSO T/2 MSO T/2 HDO 3T/2 MSI T/2 MLA T/2 MSI T/2 HDO T/2 HRE T/2 HDO 1T MSO 2T ;
: signed
CTONE ADAGIO Signal 500 ms Signal
s" DeLaswedge! This, Is, Moscow, Radio! Station! " girlsay ;
: main
cr ." ===== FORTH style program outcomes ===== " cr
s" FORTH style program outcomes " girlsay
mainF
cr cr ." ===== BASIC style program outcomes ===== " cr
s" BASIC style program outcomes " boysay
mainB
cr ." ============ the end =================== " cr
s" The end " boysay
signed
;
\ ===================================================
\s
main
===== FORTH style program outcomes =====
Counts = 46
Solution = -59.99999999
Counts = 45
Solution = 67.291782187
Counts = 45
Solution = 88.517724806
===== BASIC style program outcomes =====
Steps = 46
Solution = -59.99999999
Steps = 45
Solution = 67.291782187
Steps = 45
Solution = 88.517724806
============ the end ===================
第(92)個範例與第(17)個範例的內容相同,都是用來求解整數元素之方形矩陣的行列式值程式。
這樣的安排,對一般人而言,不太具有什麼意義,對我自己而言,則有深層的意義。
在第(17)範例說明中,已經介紹過程式設計的最重要精神,是不使用叫用自用副程式的方法來完成設計,而是根據能夠處理矩陣的數學原理,採用調整出等效上三角矩陣的方式解題,最後,把對角線上的所有元素之值乘起來,就是行列式值了。
我在設計範例時,覺得在校的學生,所能遇到的題目,都是只給整數元素,然後要求解題的,因此,在設計程式時,也採用這樣的方式演進程式。在實務工程中,很少有這麼簡單的題目,例如我在面對3D繪圖座標的工程計算時,都是從浮點數座標作為元素開始,反而是在演算完畢後、想繪點時,才需要轉換成整數座標來繪圖,實務與教材程序恰恰相反。
為了讓程式永久有用,這個範例中加寫了許多的說明與註釋,以便讓程式在未來無論是遇到整數元素或遇到浮點數元素時,都能使用。
還有,請注意這一套解題方法,在程式中,無論如何,都得採用浮點數運算,就算最後只想得到整數的輸出結果,演算過程中,也仍得使用浮點運算。請注意,我在 Lina64 中加裝的浮點系統是我自創的,一個能在傳統浮點系統中算得正確的問題與程式,在我自創的浮點系統中,當然也得獲得同樣的正確結果。我在 2008 年寫成的舊程式,在十年之後的 2018 年,再度拿出來使用,也必須有同樣的正確結果,才能表示我創作的理念完全正確,這個範例證明了這樣的成就。
除此之外,我並不滿足於程式僅只能測試這麼小規模的尺度。後面第(95)個範例,您將可以見到,能夠全自動提供保證是線性獨立、龐大無比的矩陣元素之獨特技術。(92)與(95)兩套範例,結合起來使用,我的系統就可以測試出幾百乘幾百的方形矩陣,顯示出照樣能得到正確的結果。這種事情,只能我自己做,自己驗證了,這是另一種成就,我辦到過。
準備被永久使用的設計,將程式分列成幾個區塊(heap),日後再度使用時,只須在數據輸入區塊內更換數據,就可以執行出整個程式得到結果。這個程式我已用過無數次。
:
\ A procedure to evaluate a determinant by the leading diagonal method, using largest pivots. 2008-08-23
\ Lina64 version 20180508
\ ******************************
\ (1).Data structure declared heap
6 integers N I J K P Q
2 reals a d
20 20 (MATRIX) IA
20 20 MATRIX AA
\ ******************************
\ (2).Data input heap
: INPUT-IDATA
[[ N = 4 ]]
[[ IA ( 1 1 ) = 1 ]] [[ IA ( 1 2 ) = 2 ]] [[ IA ( 1 3 ) = 3 ]] [[ IA ( 1 4 ) = 4 ]]
[[ IA ( 2 1 ) = 2 ]] [[ IA ( 2 2 ) = 3 ]] [[ IA ( 2 3 ) = 4 ]] [[ IA ( 2 4 ) = 1 ]]
[[ IA ( 3 1 ) = 3 ]] [[ IA ( 3 2 ) = 4 ]] [[ IA ( 3 3 ) = 1 ]] [[ IA ( 3 4 ) = 2 ]]
[[ IA ( 4 1 ) = 4 ]] [[ IA ( 4 2 ) = 1 ]] [[ IA ( 4 3 ) = 2 ]] [[ IA ( 4 4 ) = 3 ]]
;
: INPUT-DATA
INPUT-IDATA
basic
10 for i = 1 to N
20 for j = 1 to N
30 let { AA ( i j ) = i>r ( IA ( i j ) ) }
40 next j
50 next i
60 end ;
\ *******************************
\ (3).Auxiliary instructions heap
: SHOW-AA
BASIC
10 run cr
20 FOR I = 1 TO N
30 FOR J = 1 TO N
40 let { a = AA ( I J ) }
50 run a fs. 3 spaces
60 NEXT J
70 run cr
80 NEXT I
90 END ;
\ ********************************
\ (4).Main code heap
1.0 e -12 fconstant feps
: (DET)
BASIC
\ 10 RUN INPUT-DATA
20 LET { D = f1.0e0 }
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 250
240 GOTO 300
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 ) < feps } 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 ) * feps } THEN 380
370 GOTO 390
380 LET { AA ( I J ) = f0.0e0 }
390 NEXT J
400 NEXT I
\ End of reduction to upper triangular form.
410 NEXT K
420 LET { D = D * AA ( N N ) }
430 IF { ABS ( D ) > feps } THEN 450
440 LET { D = f0.0e0 }
450 LET { D = D }
\ Determinant is now in D
\ 460 PRINT { D }
470 END ;
\ *********************************
\ (5).Manipulation heap
: round ( f1 -- f2 )
over abs nlog 12 > if fround then ;
: DET
BASIC
10 RUN INPUT-DATA
20 run cr ." matrix " n . ." X " n . ." : " cr show-aa cr
30 IF N = 1 THEN 50
40 GOTO 70
50 LET { D = AA ( 1 1 ) }
60 GOTO 120
70 IF N = 2 THEN 90
80 GOTO 110
90 LET { D = AA ( 1 1 ) * AA ( 2 2 ) - AA ( 1 2 ) * AA ( 2 1 ) }
100 GOTO 120
110 RUN (DET)
120 PRINT { " determinant = " ; D }
130 let { d = d round }
140 run cr ." determinant = " d f. cr
150 END ;
det
\ ********** the end ***************
\ det = 160
第(93)個範例是一個傅立葉(Fourier)級數轉換的數值解法,請注意!不是快速傅立葉級數轉換(FFT)的課題。
這個範例程式的詳細解說,包括引用公式,請參考我個人網頁 20180502 的貼文: Coefficients of Fourier series 。該文是全用英文寫的,文中也有實作成果的驗證貼圖,足以用來解釋整個程式的基本精神。
當初,我在發展語音輸出技術時,接觸到套裝軟體的『輸入規格』。規格中要求給予將近 20 個 sin 或 cos 函數的係數,作為語音輸出時的指定係數,套裝軟體就能發出各種聲音。
究其原因,就是各種聲音的波譜函數,均能重新表示成傅立葉級數展開式,所以,只要取得級數各項的係數,代入程式,就能令語音合成晶片發出同樣波譜的聲音。
音調的高低,只需設定單個頻率參數就能改變。聲音的大小,只是音量放大的設定。持續的時間,也只是設定共振週期的長度便可辦到。套裝軟體的功能一應俱全,說明資料也很完整,軟體還有許多套,套套都很容易使用,規格也都相當近似。
這麼簡單的語音輸出道理,引起我的興趣,於是就寫了這個程式。為了解釋效果,我也再度利用圖示,顯示理論效果。發現,以此方式完成的工作,簡直就跟原來的波譜完全相似,效果非凡。
換句話說,如果我有心從事於這項工作,絕對能把任何歌星的聲音重現,而且有把握讓聲音以清唱的方式,從語音輸出晶片發送出去。
這一生可以研究的主題很多,不缺這種課題,因為我寫出這個範例程式、看過效果之後,覺得學問並不高深,但要做這方面的事情時,必須花費非常多的時間才能完成。
語音的花樣跟數學題目一樣,是無限多個。有用的一整套規模,至少也得發展出一套中文語音輸出所需的數量,才有值得稱道的地方。我在網文上表示過,這應該是國家政策上的計畫,不是我個人該做的事情。想找齊標準男生與女生的音源,就不是憑我個人能力所能及的事情。還要系統性的造出許多必要的波譜,才能給我這種程式把係數分離出來,提供套裝軟體使用。後續的控制,就更有得設計的了,最後才能成為有用的產品。可是,整個工作,需要用到電腦計算,取得係數的地方,只需像我寫的這個範例程式,就夠用了。別人也許也有更特別的方法,最後,也不外乎取得係數。
這麼說來,請看看我的貼圖,體會一下數學計算系統的能力,把一個很尖銳的週期性三角形波譜,全用傅立葉級數展開式重新表示時的差別,也就夠了。
我們在大學裡學過的傅立葉級數,是能夠直接手解解出三角波的理論係數的,這個範例程式沒有直接取用手解解答所得係數來繪圖,是真正使用電腦的數位解法所得之數值係數繪出的圖形,只取五項,我只在單一個週期內繪圖,就已經能夠有這種效果了。
程式的輸入數據,可以是任意的波譜讀數。點數也沒有限制,點取得多了,解析度就高。所用到的理論式子,大二的工程數學中都有。我因長期從事過信號處理的軟體研發工作,熟悉這方面數學與軟體結合起來的應用技術。
:
\ Calculating Fourier coefficients in Lina64 ABC forth
\ ===============================================
\ (1) Data structure declare heap
12 value km
km 2 / value nm
2 integers k n
km array f
km array t
km array u
nm array a
nm array b
5 reals c d q r m
\ ===============================================
\ ===============================================
\ (2) Data input heap
: input
{{ m = i>r ( nm ) }}
{{ f ( 0 ) = 1.0 e 0 }}
{{ f ( 1 ) = 6.666 e -1 }}
{{ f ( 2 ) = 3.333 e -1 }}
{{ f ( 3 ) = 0.0 e 0 }}
{{ f ( 4 ) = -3.333 e -1 }}
{{ f ( 5 ) = -6.666 e -1 }}
{{ f ( 6 ) = -1.0 e 0 }}
{{ f ( 7 ) = -6.666 e -1 }}
{{ f ( 8 ) = -3.333 e -1 }}
{{ f ( 9 ) = 0.0 e 0 }}
{{ f ( 10 ) = 3.333 e -1 }}
{{ f ( 11 ) = 6.666 e -1 }}
{{ f ( 12 ) = 1.0 e 0 }}
;
\ All parts above are changeable setting.
\ ===============================================
\ ===============================================
\ All parts below are fixed.
\ (3) Main code heap
: main
basic
10 run input
20 for n = 1 to nm
30 let { q = 0.0 e 0 }
40 let { r = 0.0 e 0 }
50 for k = 1 to km
60 let { c = i>r ( k ) }
70 let { d = i>r ( n ) }
80 let { t ( k ) = ( f ( k ) * cos ( ( d * c * f(pi) ) / m ) ) / m }
90 let { u ( k ) = ( f ( k ) * sin ( ( d * c * f(pi) ) / m ) ) / m }
100 let { q = q + t ( k ) }
110 let { r = r + u ( k ) }
120 next k
130 let { a ( n ) = q }
140 let { b ( n ) = r }
150 print " format 1================== "
160 print " ( " ; n ; " ):: " , { q , r }
200 print " format 2=================== "
210 print " a ( " ; n ; " ) = " ; { q }
220 print " b ( " ; n ; " ) = " ; { r }
300 print " format 3=================== "
310 run cr ." a ( " n . ." ) = " q f.
320 run cr ." b ( " n . ." ) = " r f. cr cr
400 next n
410 end ;
\ ===============================================
\ ===============================================
\ (4) Manipulation heap
: check
main
basic
20 for k = 0 to 360
30 let { c = i>r ( k ) }
40 let { d = f(2pi) * c / 360.0 e 0 }
50 let { m = a ( 1 ) * cos ( d ) + a ( 3 ) * cos ( 3.0 e 0 * d ) + a ( 5 ) * cos ( 5.0 e 0 * d ) }
60 run cr d f. 5 spaces m f.
100 next k
110 end ;
\ ================================================
第(94)個範例,也是談積分技術的範例程式。
積分技術的細部討論,已經在第(74)範例中討論過了,不再贅述。但這個範例與(74)比較,有所不同,因此,只談不同的地方。
原(74)範例,是一個純用 Forth 程式語法寫成的程式,而(74-1)則改為全用 BASIC 式的語法寫成,(94)範例則亦為 BASIC 式語法寫成。
前曾述及,積分結果,可以精確到與用來提供正式使用時的函數值完全一樣。這個性質,凸顯了積分的精確性。前面舉例,提到 ln(x) 的微分是 1/x ,所以,如果您不知道自然對數的函數值為多少時,您可以透過對 1/x 的函數進行分割成很小間隔的 dx 來積分,就能取得很精確的 ln(x) 函數值。在(74-1)中,我設計了那樣的程式,很容易重新設定,進行驗證。
在這個(94)範例中,我進一步提供另一個類似問題的函數,就是: arc tan 函數的微分是 1/(1+x**2) 。以其驗證積分可以獲得很精確答案的論述。
(94)範例是我在 2016 年,發展完 Lina64 系統後,才用它來測試系統的一個程式。它的效應,令我對自己完成的系統設計工作充滿了信心,這種測試是不可能有人願意幫忙做的事情,通常都是系統設計人得自己做,我也不例外。因此,那時,每當我測試系統告一段落之後,就在我的個人網頁上貼文,向全世界公告帶有影音、圖示以及詳細說明的結果。
這樣做,是個博士鄰居建議我積極進行的,他對技術專利很敏感,因為他自己開公司,以他自己特有的技術,設計產品賣全世界,他告訴我說,如果是自己苦心經營出來的創作,就算不賣,先行在網上刊出使用成果,就能保護自己的創作,留下親筆筆記,是更好的做法,萬一將來被人告上法庭,非打官司不可時,絕對有用。所以,我就這樣做了。這個範例的詳細說明,也能在我的個人網頁中找到,我不想花時間重新探尋是那一天刊出的了,因為每次自己去找刊出日期,就會令那個網頁的訪客統計記錄被自己的點閱而扭曲。我個人覺得,就算都沒人想看我的個人網頁了,我也不必自我陶醉,自己點開來增加點閱的統計人數。網頁的價值該是怎樣,就讓它怎樣,我自己不去扭曲應有的事實。
從微積分課本的導函數表中,可以找出許多具有上述性質的函數來進行同樣的測試,但是,導函數必須是不為本函數者,才有這樣測試的意義。例如:下列舉出幾個本函數及其導函數的關係,前者是本函數,次者是其導函數。
x**a -----> a*x**(a-1)
exp(x) ---> exp(x)
a**x -----> (a**x)*ln(a)
ln(x) ----> 1/x
sin(x) ---> cos(x)
cos(x) ---> -sin(x)
tan(x) ---> sec(x)**2
arc tan --> 1/(1+x**2)
因此,您從上列導函數表中可以看得出來,並不是各個函數之值都適合以導函數之積分值測試的。設計系統時,基本函數可能是從另一個更基本的函數導出來的,所以,如果導函數仍為本函數者,如 : exp(x) 函數,就不適合用來從事本範例中強調的測試,否則將會陷入所謂的『循環驗證的誤謬』,我這樣介紹測試,必須先把原則搞清楚。
Lina64 ABC Forth 系統的測試結果列示於後:
請從間隔量為一千點,與間隔量提升為一萬點後的積分結果中,體會所得積分之值,與實際本函數之函數值比較,驗證其精確度是否提高了好幾位數?
:
\ Integrations in lina64 ABC FORTH
\ Integration.f
\ Author: Ching-Tang Tseng
\ 20161107, Hamilton, NZ
9 reals a b h sum1 sum2 x f(x) f(a) f(b)
2 integers n i
1000 value approximations
: 1/x ( -- )
basic
10 let { f(x) = 1. E 0 / x }
20 end ;
: 1/(1+x**2) ( -- )
basic
10 let { f(x) = 1. E 0 / ( 1.0 e 0 + x * x ) }
20 end ;
defer function
' 1/x is function
: range1/x ( -- )
basic
10 let n = approximations
20 let { a = f1.0e0 }
30 let { b = 100.0 e 0 }
40 end ;
: range1/(1+x**2) ( -- )
basic
10 let n = approximations
20 let { a = f0.0e0 }
30 let { b = f2.0e0 }
40 end ;
defer setting
' range1/x is setting
\ All parts above are changeable setting.
\ All parts below are fixed.
: SimpsonInit ( -- )
basic
10 run Setting
20 let { h = ( b - a ) / I>R ( n ) }
30 let { x = a + h / 2. E 0 }
40 run function
50 let { sum1 = f(x) }
60 let { sum2 = 0. E 0 }
70 let { x = a }
80 run function
90 let { f(a) = f(x) }
100 let { x = b }
110 run function
120 let { f(b) = f(x) }
130 end ;
: Simpson ( -- )
basic
10 run SimpsonInit
20 for i = 1 to n - 1
30 let { x = a + h * I>R ( i ) + h / 2. E 0 }
40 run function
50 let { sum1 = sum1 + f(x) }
60 let { x = a + h * I>R ( i ) }
70 run function
80 let { sum2 = sum2 + f(x) }
90 next i
100 print " Simpson integration = " ;
{ ( h / 6. E 0 ) * ( f(a) + f(b) + 4. E 0 * sum1 + 2. E 0 * sum2 ) }
110 end ;
: L-Rectangle ( -- )
basic
10 run Setting
20 let { h = ( b - a ) / I>R ( n ) }
30 let { x = a }
40 let { sum1 = 0. E 0 }
50 for i = 1 to n
60 run function
70 let { sum1 = sum1 + h * f(x) }
80 let { x = x + h }
90 next i
100 print " L-Rectangle integration = " ; { sum1 }
110 end ;
: R-Rectangle ( -- )
basic
10 run Setting
20 let { h = ( b - a ) / I>R ( n ) }
30 let { x = a }
40 let { sum1 = 0. E 0 }
50 for i = 1 to n
60 let { x = x + h }
70 run function
80 let { sum1 = sum1 + h * f(x) }
90 next i
100 print " R-Rectangle integration = " ; { sum1 }
110 end ;
: M-Rectangle ( -- )
basic
10 run Setting
20 let { h = ( b - a ) / I>R ( n ) }
30 let { x = a }
40 let { sum1 = 0. E 0 }
50 for i = 1 to n
60 let { x = x + h / 2. E 0 }
70 run function
80 let { sum1 = sum1 + h * f(x) }
90 let { x = x + h / 2. E 0 }
100 next i
110 print " M-Rectangle integration = " ; { sum1 }
120 end ;
: Trapezoidal ( -- )
basic
10 run Setting
20 let { h = ( b - a ) / I>R ( n ) }
30 let { sum1 = 0. E 0 }
40 let { x = b }
50 run function
60 let { f(b) = f(x) }
70 let { x = a }
80 run function
90 let { f(a) = f(x) }
100 let { sum1 = h * ( f(b) + f(a) ) / 2. E 0 }
110 for i = 1 to n - 1
120 let { x = x + h }
130 run function
140 let { sum1 = sum1 + h * f(x) }
150 next i
160 print " Trapezoidal integration = " ; { sum1 }
170 end ;
\ Output description depends on input setting.
: integrations ( -- )
L-Rectangle
R-Rectangle
M-Rectangle
Trapezoidal
Simpson
;
: function1
cr cr
." Integration of 1/x from 1 to 100 = ln(100) "
cr
integrations
cr cr 18 spaces
." ln(100) = " {{ ln ( 100.0 e 0 ) }}
fs. cr ;
: function2 ( -- )
cr cr
." Integration of 1/(1+x^2) from 0 to 2 = atan(2) "
cr
integrations
cr cr 18 spaces
." atan(2) = " {{ atan ( f2.0e0 ) }}
fs. cr ;
: (main)
['] 1/x is function ['] range1/x is setting function1
['] 1/(1+x**2) is function ['] range1/(1+x**2) is setting function2
;
: main
cr ." =====Integration===== "
cr ." approximations = 1000 "
1000 to approximations (main)
cr ." ========================================= "
cr ." approximations = 10000 "
10000 to approximations (main)
cr ." =======the end======= "
;
main
\s
ching@center:~$ cd lina64
ching@center:~/lina64$ ./f5124
AMDX86 ciforth 5.3.0
fload 94
A : ISN'T UNIQUE
B : ISN'T UNIQUE
I : ISN'T UNIQUE
=====Integration=====
approximations = 1000
Integration of 1/x from 1 to 100 = ln(100)
L-Rectangle integration = 4.6549910575 E 0
R-Rectangle integration = 4.5569810575 E 0
M-Rectangle integration = 4.6047625486 E 0
Trapezoidal integration = 4.6059860575 E 0
Simpson integration = 4.6051703849 E 0
ln(100) = 4.6051701859 E 0
Integration of 1/(1+x^2) from 0 to 2 = atan(2)
L-Rectangle integration = 1.1079486644 E 0
R-Rectangle integration = 1.1063486644 E 0
M-Rectangle integration = 1.1071487444 E 0
Trapezoidal integration = 1.1071486644 E 0
Simpson integration = 1.1071487177 E 0
atan(2) = 1.1071487177 E 0
=========================================
approximations = 10000
Integration of 1/x from 1 to 100 = ln(100)
L-Rectangle integration = 4.6100788525 E 0
R-Rectangle integration = 4.6002778525 E 0
M-Rectangle integration = 4.6051661027 E 0
Trapezoidal integration = 4.6051783525 E 0
Simpson integration = 4.605170186 E 0
ln(100) = 4.6051701859 E 0
Integration of 1/(1+x^2) from 0 to 2 = atan(2)
L-Rectangle integration = 1.1072287172 E 0
R-Rectangle integration = 1.1070687172 E 0
M-Rectangle integration = 1.107148718 E 0
Trapezoidal integration = 1.1071487172 E 0
Simpson integration = 1.1071487177 E 0
atan(2) = 1.1071487177 E 0
=======the end======= OK
第(95)個範例是能解出一組魔術方陣的程式,我有個人用途。
關於魔術方陣,網上可以找到許多參考文獻,魔術方陣是中國人最先發現的,三階魔術方陣又稱『洛書』,已有兩、三千年的歷史。魔術方陣的定義是把由 1 到 n 平方的整數,排成一個方陣(稱為 n 階方陣),使每一行、每一列,及兩條對角線上 n 個數字之和都相同。
為了展現中國古文化與現代科技的關聯性,我特地利用我家的圍棋,擺出三張所謂『洛書』的棋譜照相後貼文,這就代表了魔術方陣的解,恆『存在』但都不是『唯一』的。每一維度的方陣之解,都可以有很多個,解法也有很多種,而且可以用排列、組合的數學分析技術,算出到底可能有多少種解?這方面的問題,不是本文打算討論的重點。
我在 20180602 貼出的個人網文中,原本只貼程式執行後的輸出結果,後來才加補貼我自己設計出來的這套源程式。我在我的個人網頁上,經常對貼文修改內容,用來測試跟隨網文的系統會不會變化?後來確定,我每改一次網文內容,電腦系統確實是會重新掃描該網頁,可能是為了要修正 Google 公司的關鍵字資料庫系統之內容吧。
網路上也有很多能解出魔術方陣的程式,但很少有人能寫出隨便給一個正整數,就都能有答案的程式。我對問題分析後,看出只需把整數分成三大類,從 3 開始,所有 3 以上的整數,就都能憑程式解出一套解答來,就自創了這套程式。
從純數學的立場看這個主題,魔術方陣算是一種研究數字理論的課題。從三千年前中國人的立場看這個主題,魔術方陣是一種配合甲骨文來搞算命卜卦的工具。從現行許多媒體上經常刊載的魔術方格填字遊戲看這個主題,魔術方陣是一種電玩的項目。
從我個人的觀點呢?我把它當作用來產生矩陣程式所需之測試樣品的絕佳工具。大家請自己看程式執行後的輸出結果,每一列、每一行的元素,全都不一樣,這代表了什麼意義?這種矩陣就絕對是一種『線性獨立』形式的矩陣,不是嗎?拿來測試矩陣運算系統最好用了。否則,您要自己亂填出元素之值,又不能確保可能會『線性相依』,很麻煩。
至於程式本身,涉及我自己創作的部分,有個要點必須提醒大家:
特別注意!我所設計的 ABC Forth 系統中,矩陣或陣列的指標,不可以為運算式。凡運算式都得重新以高階定義方式特別指定,如本程式內所列出的三個指標: x+sqd2 , y+sqd2 , l+1 便是。
這個範例是個很特別的設計,要讀過很多資料後,才設計得出來。如果覺得沒有用,請自己印個方陣拆掉一些數字,當個電玩遊戲玩一玩就算了。
:
\ magic square matrics generator
\ Author: Ching-Tang Tseng
\ 20150528, Hamilton, New Zealand
\ ******************************
\ (1).Data structure declared heap
\ Attention! max is a reserved word of mathematical function.
10 integers x y l n q sqd2 q2 nr maxx msum
100 100 (matrix) sq
\ ***********************************
\ (2).Data input heap
: 4n+2init
basic
\ 10 let n = 6
20 let msum = n * ( n ** 2 + 1 ) / 2
30 let sqd2 = n / 2 :: q2 = sqd2 ** 2
40 let l = ( n - 2 ) / 4
50 let x = sqd2 / 2 + 1 :: y = 1
60 let nr = 1
70 end ;
: 2n+1init
basic
\ 10 let n = 7
20 let nr = 1
30 let x = n - ( n / 2 )
40 let y = 1
50 let msum = n * ( n ** 2 + 1 ) / 2
50 let maxx = n * n
60 end ;
: 4ninit
basic
\ 10 let n = 8
20 let msum = n * ( n ** 2 + 1 ) / 2
30 let q = n / 4
40 let nr = 1
50 end ;
\ ********************************
\ (3).Auxiliary instructions heap
: check-columms-and-rows
basic
10 For y = 1 To n
20 let nr = 0 :: q = 0
30 For x = 1 To n
40 let nr = nr + sq ( x y )
50 let q = q + sq ( y x )
60 Next x
70 If nr <> msum Or q <> msum Then 90
80 goto 100
90 run cr x 4 .r y 4 .r cr -1 abort" Error1: value <> magic sum " cr
100 Next y
110 end ;
: check-diagonals
basic
10 let nr = 0 :: q = 0
20 For x = 1 To n
30 let nr = nr + sq ( x x )
40 let q = q + sq ( x ( n - x + 1 ) )
50 Next x
60 If nr <> msum Or q <> msum Then 80
70 goto 90
80 run cr nr . 3 spaces q . 3 spaces cr -1 abort" Error2: value <> magic sum " cr
90 end ;
: EraseSq
basic
10 for y = 1 to 100
20 for x = 1 to 100
30 let sq ( x y ) = 0
40 next x
50 next y
60 end ;
: PrintSq
basic
10 Print " Magic square size : " ; n ; " * " ; n
20 Print " The magic sum = " ; msum
30 run cr
40 For y = 1 To n
50 For x = 1 To n
60 let q = sq ( x y )
70 run q 4 .R
80 Next x
90 run cr
100 Next y
110 end ;
\ *******************************
\ (4).Main code heap
: 2n+1
basic
10 if sq ( x y ) = 0 then 30
20 goto 100
30 let sq ( x y ) = nr
40 if nr mod n = 0 then 70
50 let x = x + 1
55 let y = y - 1
60 goto 80
70 let y = y + 1
80 let nr = nr + 1
100 if x > n then 120
110 goto 200
120 let x = 1
130 if sq ( x y ) <> 0 then 150
140 goto 200
150 let x = x + 1
200 if y < 1 then 220
210 goto 300
220 let y = n
230 if sq ( x y ) <> 0 then 250
240 goto 300
250 let y = y - 1
300 if nr > maxx then 400
310 goto -10
400 end ;
: 4n
basic
10 for y = 1 to n
20 for x = q + 1 to n - q
30 let sq ( x y ) = 1
40 next x
50 next y
60 For x = 1 To n
70 For y = q + 1 To n - q
80 let sq ( x y ) = sq ( x y ) xor 1
90 Next y
100 Next x
110 let q = n * n + 1
120 for y = 1 To n
130 for x = 1 To n
140 if sq ( x y ) = 0 Then 170
150 let sq ( x y ) = nr
160 goto 180
170 let sq ( x y ) = q - nr
180 let nr = nr + 1
190 Next x
200 Next y
210 end ;
\ Attention! in ABC FORTH, index can not to be an algebraic experssion
: x+sqd2 x sqd2 + ;
: y+sqd2 y sqd2 + ;
: l+1 l 1+ ;
: 4n+2
basic
10 if sq ( x y ) = 0 then 30
20 goto 130
30 let sq ( x y ) = nr
40 let sq ( x+sqd2 y+sqd2 ) = nr + q2
50 let sq ( x+sqd2 y ) = nr + q2 * 2
60 let sq ( x y+sqd2 ) = nr + q2 * 3
70 if nr mod sqd2 = 0 then 90
80 goto 110
90 let y = y + 1
100 goto 120
110 let x = x + 1 :: y = y - 1
120 let nr = nr + 1
130 if x > sqd2 then 150
140 goto 200
150 let x = 1
160 if sq ( x y ) <> 0 then 180
170 goto 200
180 let x = x + 1
200 if y < 1 then 220
210 goto 260
220 let y = sqd2
230 if sq ( x y ) <> 0 then 250
240 goto 260
250 let y = y - 1
260 if nr > q2 then 300
270 goto -10
300 for y = 1 to sqd2
310 for x = 1 to l
320 let q = sq ( x y )
330 let sq ( x y ) = sq ( x y+sqd2 )
340 let sq ( x y+sqd2 ) = q
350 next x
360 next y
400 let y = ( sqd2 / 2 ) + 1
410 let q = sq ( 1 y )
420 let sq ( 1 y ) = sq ( 1 y+sqd2 )
430 let sq ( 1 y+sqd2 ) = q
440 let q = sq ( l+1 y )
450 let sq ( l+1 y ) = sq ( l+1 y+sqd2 )
460 let sq ( l+1 y+sqd2 ) = q
500 for y = 1 to sqd2
510 for x = n - l + 2 to n
520 let q = sq ( x y )
530 let sq ( x y ) = sq ( x y+sqd2 )
540 let sq ( x y+sqd2 ) = q
550 next x
560 next y
1000 end ;
\ **********************
\ (5).Manipulation heap
: magics ( n -- )
[[ n ]] !
EraseSq
basic
10 if n mod 2 = 1 then 100
20 if n mod 4 = 2 then 200
30 goto 300
100 run 2n+1Init
110 run 2n+1
120 goto 400
200 run 4n+2Init
210 run 4n+2
220 goto 400
300 run 4nInit
310 run 4n
400 run PrintSq
500 run check-columms-and-rows
510 run check-diagonals
1000 end ;
: main ( -- )
21 3
do i magics loop ;
main
\ ********** the end ***************
一百個例題 (86 ~ 90)
Ching-Tang Tseng
Hamilton, New Zealand
5 October 2024
第(86)個範例是仿造其他程式語言所強調之巨集指令組(macro set)觀念而設計的一個類似範例。
實際上, Forth 本身的程式設計方式,根本就各個都是巨集指令,都在設計完成後,就能給後續指令當巨集指令使用。所以,這個範例實際上是個多餘的累贅研究。
用過別種程式語言者,經常在闖入 Forth 國際論壇網頁時,胡亂發言,拿一般程式語言中使用過的術語來考驗 Forth 的熟手,通常都不太有人願意理會,但挑釁久了,也難免會有高手或老手出面修理來犯者,於是,就有這種無聊程式可被收集,這個範例就是這麼來的。
不過,要想設計出這個範例中的 macro 指令,沒有三級 Forth 程度的使用能力,也並不容易一下子就寫得出程式。程式中使用了延遲(postpone)系統執行時態(tense)的技巧解決問題。
所謂的時態,指的是系統經常只能處於執行時態(execution tense)或編譯時態(compiling tense),否則就是當機時態(exeception tense)。如果系統處於當時的時態,而程式內容能令系統改變時態時,就叫做延遲系統的執行時態了。通常的用法,都是在程式中使用一個叫做 postpone 的指令,要系統在執行到它的後續指令時,不要執行,必須改在下一個執行時態時才執行。
這樣的哲理解釋,第一次碰到的人,很難理解它的意義,因為我沒拿實際例子來解釋。
所有相關於這一類指令的詳細說明,請參考鄒靖寧老師編譯的『FORTH入門』第十一章,這本『Starting FORTH』教材中的最後一章才有解說,那時的指令名稱叫 COMPILE , [COMPILE] ,現在叫 POSTPONE 。由於這方面的哲理深奧難懂,才會被編寫到最後一章。也由於指令功能效果非凡,所以才會令 Forth 程式語言與眾不同、功能強大無比,高手也才能拿它來發揮、再創新高。
我不諱言,我自創的系統,很多地方都涉及這方面的技術應用,才可以令指令被執行的時態強行改變,也才能改變程式的語法,我才寫得出整套的 BASIC ,可以不像傳統的 BASIC 作者那樣,非得寫 BASIC 的剖析器(parser)不可,我自創的系統中直接編譯 BASIC 程式,執行起來才會快速,使用者能感覺得出程式是被直接執行的。
這個巨集指令的設計方法,是將隨後需要用到的現成指令,全當字串來看待,設計時都不會被執行,等到下一次巨集指令再被用到時,就把這些字串當指令來編進系統。執行最後完成的 foo 測試指令時,巨集效果就產生了。把字串整串執行掉的指令叫 evaluate ,執行時立即安放字串的指令叫 sliteral 。
我所謂的累贅是這樣解釋的:如果不設計這個 macro ,您就用 『 : square dup * ; 』規規矩矩地寫完 Forth 標準設計,而在此處的 foo 測試程式,仍然是寫法不變,就能完成。那麼,在 Forth 中要 macro 做什麼?根本就可以完全不用。
我在 Win32Forth 系統的低階設計部分,發現後來的版本,前人很好的源程式被大量改寫了,改為使用許多內容為低階組合語言的巨集指令,令程式更難閱讀,有錯就更為難改,搞得我後來連追蹤、查證的興趣都沒有了。
:
\ (86)macro.f This macro coding style is effect in Lina64 system.20171201
: macro
: char parse postpone sliteral postpone evaluate
postpone ; immediate
;
\ Usage is e.g.
macro square " dup *" ok
: foo 5 square . ; ok
foo 25 ok
第(87)個範例是一個用來印出 16 進制與 2 進制雙整數的工具指令程式。
這是一個我在發展 Lina64 系統期間所開發的發展工具,用來檢驗雙整數存在於系統中的情況。
因為在發展期間,最方便的系統觀察指令,就是 DUMP ,他能按照所給的指定位址,把後續的多少個位元組內的內容,以 16 進制的方式傾印出來。
我若在程式中輸入了雙整數,想要看看被編譯進系統後的結果,以 Dump 指令來看,只能是 16 進制的顯示,想直接比對,使用這個 DH. 與 DB. 指令來操作就很方便。
這兩個指令若想在 Lina64 中運轉得起來,有一個名叫 CO 的指令需要用到,因為要操作到文字視窗畫面來顯示輸出,它會與一般 Forth 系統不同,所以該系統中增列了這個叫做 CO 的指令,我從系統說明書中的 8.3.5 節,把說明資料也收集在範例程式的最後面,以利參考。
這個範例程式很簡單,用 Forth 系統的人,必須能直接讀懂,裡面用到的設計觀念,只是強調在改變進制時,借助於回返堆疊能夠暫存數據的能力來設計程式。臨時放一下,用完就取回來,恢復系統原來的進制狀態,如此而已。
所有的 Forth 系統在轉換記憶體內容成可印出數字時,都仍使用自古以來就用的 <# , #s , #> 等指令,只是從 1994 年的 ANSI 標準出現以後,系統就只會處理正的數字,不管負的數字了,使用時要自己注意。
反過來解釋我為什麼要設計這樣的工具?如果不用這個 DH. ,那麼,我想要檢查輸入雙整數 1.23 後的花樣時,就得執行下列操作才看得到:
1.23
Hex d. Decimal
然後對每一個輸入的雙整數,都得重複上列的操作。我嫌老打那些更換進制的命令太麻煩,所以設計了這兩個不做進制操作就能顯示結果的指令。不怕麻煩的人,可以不用理會這個設計。
:
\ (87)DH. DB. in Lina64 with special word: CO.20171202
\ Switch to hex for the duration of the definition.
: HEX: R> BASE @ >R >R HEX CO R> BASE ! ;
: BINARY: R> BASE @ >R >R 2 BASE ! CO R> BASE ! ;
\ An example
: (DH.) HEX: <# #S #> ;
: (DB.) BINARY: <# #S #> ;
: DH. (DH.) type ;
: DB. (DB.) TYPE ;
\ 1.2 (DH.) TYPE
\ C OK
\ 1.2 DH.
\ C OK
\ 8.3.5 CO
\ Name: CO
\ No stackeffect
\ Attributes:
\ Description: Return to the caller, suspending interpretation of the current definition, such that
\ when the caller exits, this definition is resumed. The return stack must not be engaged, such as
\ between >R and R> , or DO and LOOP .
第(88)個範例是一個解常數係數微分方程(Ordinary differential equations, ODE)的典型程式。
前曾述及,電腦的積分算得很準,微分算得不太準,今天就來看看解微分方程的程式技術是怎麼設計的?
範例(88-1)的程式規格才是勉強能夠符合 Win32Forth 系統的程式寫法,但標號 110 , 140 的程式中,浮點數的方次符號 ** 必須改回『 ^ 』。
範例(88)是 Lina64 中 ABC FORTH 的寫法,跟(88-1)最大的不同,就是所有浮點數的輸入方式有很明顯的差異。
凡是 1.23e4 的格式,都得以空格分開成 1.23 e 4 的方式輸入。
從這樣的差異中,大家就能體會出任何數字都只有兩種格式,一個是帶小數點的數字,另一個就是不帶小數點的整數。浮點數就可以用這種方式存在,我創作的這種浮點系統,辦到了這樣的設計,好用的很。如果不這樣設計,您可以仔細追蹤 Win32Forth 系統是如何接受浮點數的,把源程式印出來看,需要四張 A4 用紙,不是高手就很難看得懂。我在那個所有的浮點數都需要用到的 e 或 E 上動了手腳,一列程式就設計出來了,我自創的系統就能這麼方便的輸入浮點數。
早期的發展,我沒顧及後來別種程式語言都用 ** 代表 ^ 了,至少在我放棄 Win32Forth 之 ABC Forth 前,我都仍照用 ^ 。
簡單介紹微分方程電腦程式解法的大要:
教科書中的常微方程,都是經過很理想的安排,寫得整整齊齊,才能那麼容易用來教學生手解,這樣的解法,在學術上叫解析解。我踏入社會後,如果碰得到微分方程問題,就都不可能是理想的方程式了,很多 ODE 都還是根本無法手解的東西,必需簡化掉許多項後才能手解,解完,才又考慮再帶有另一項時會是怎樣?再來修正,有時解答使用的函數型式要靠猜的,才能勉強含有意義。一般手解解法的步驟,是首先把 ODE 的通解先解出來,解答便帶有未定係數,然後將起始條件或邊界條件帶入通解,再把未定係數解出來,這樣,就得到了明確的特解,算是最後結果。
用電腦程式解微分方程時,無論是 ODE 或者是偏微方程(PDE) 都一樣,解法就跟上述步驟完全不同。電腦程式解法是一開始就把無論是能不能手解的微分方程式,全都改寫成數學家憑理論推導出來的差分方式表示之等效近似逼近式,然後,直接就從起始條件或邊界條件作為演算推進的起點,一點一點的把整個指定領域中所有的對應點之數值給算出來,所以,解答都是一長串的數據組,例如一個 x 座標值對應到一個 y 的函數值。
從這麼簡單的說明,大家可以明白,電腦解法列印出來的數據是函數的趨勢,兩點之間的事情不列入考慮。因此,根本沒有值得強調有幾位數精確度的意義。把這些數據畫出圖來,才能看得清楚解答的變化趨勢。我因為長期搞電腦的數學運算,也長期使用電腦畫函數曲線圖,所以熟知這些效果,講得出其中的道理。
我保證上訴文字,您此前在那裡都沒見過。
那麼,電腦解微分方程有何好處?有的,還真不少。前面說過,真實世界中的微分方程,大部分都是無法手解、很不理想的方程式,用電腦來解這種工程問題時,沒有困擾,不需要背記那一種方程式該用那一種有點碰運氣之猜得的三角或指數函數來匹配,可以直接寫逼近式,有時需要預作微分推導,會比較麻煩而已。微分方程的數值分析解法技術也已有很多種,也可以不用進行微分處裡就能解題的應用技術,這個範例只是眾多技術中的一種,叫作泰勒級數四次微分後的匹配解法。
解的過程在起始步驟中就得先給予指定範圍與每一步驟的推進差量,當然,差量越小就可以得到越精確的解答,只是電腦計算會多一些,印出來的數據也當然多一些,如此而已。最後,電腦繪圖才最有用,趨勢用圖一看就能一目暸然。
我曾經接觸到用電腦程式解微分方程的第一個問題,是計算出從太空返回地球時,物體自由落體狀態下的極限速度。一開始,只有自己設計的固定點浮點系統可用,但就算解得的輸出數據只有三位數精確度可言,拿來繪圖,卻也看得出很理想的變化趨勢了。
這幾十年我解過無數次的微分方程式,自己覺得就算拿來搞武器的線上直接控制,也能有用。掌握趨勢才是解微分方程的重點,小數點後面三位數,就有千分之一的精確度,打什麼目標需要這麼準?我搞過打軍艦的魚雷,解析度才 12 度左右就夠了,把我設計的系統裝上武器,可以直接使用。
:
\ (88)Taylor series method to solve ODE(order 4 ) in lina64 ./f5106
\ data M,h,t,x/200,0.01,-1.0,3.0/
2 integers k m
7 reals h t x x1 x2 x3 x4
: taylor basic
10 print " * "
20 print " *, Taylor series method (order 4) "
30 print " *, Section 8.1, Kincaid-Cheney "
40 print " * "
50 let m = 200
60 let { h = 0.01 e 0 }
70 let { t = -1.0 e 0 }
80 let { x = 3.0 e 0 }
90 print 0 , { t , x }
100 for k = 1 to m
110 let { x1 = cos ( t ) - sin ( x ) + t ** 2.0 e 0 }
120 let { x2 = negate ( sin ( t ) ) - x1 * cos ( x ) + 2.0 e 0 * t }
130 let { x3 = negate ( cos ( t ) ) - x2 * cos ( x ) + ( x1 * x1 ) * sin ( x ) + 2.0 e 0 }
140 let { x4 = sin ( t ) + ( x3 *8 3.0 e 0 - x3 ) * cos ( x ) + 3.0 e 0 * x1 * x2 * sin ( x ) }
150 let { x = x + h * ( x1 + ( h / 2.0 e 0 ) * ( x2 + ( h / 6.0 e 0 ) * ( x3 + ( h / 24.0 e 0 ) * x4 ) ) ) }
160 let { t = t + h }
170 print k , { t , x }
200 next k
210 end ;
\ (88-1)Taylor series method to solve ODE(order 4) in gforth64 abc forth
2 integers k m
7 reals h t x x1 x2 x3 x4
: taylor ( -- )
basic
10 print " * "
20 print " * Taylor series method (order 4) "
30 print " * Section 8.1, Kincaid-Cheney "
40 print " * "
50 let m = 200
60 let { h = 0.01e0 }
70 let { t = -1.0e0 }
80 let { x = 3.0e0 }
90 print 0 , { t , x }
100 for k = 1 to m
110 let { x1 = cos ( t ) - sin ( x ) + t ** 2.0e0 }
120 let { x2 = negate ( sin ( t ) ) - x1 * cos ( x ) + 2.0e0 * t }
130 let { x3 = negate ( cos ( t ) ) - x2 * cos ( x ) + ( x1 * x1 ) * sin ( x ) + 2.0e0 }
140 let { x4 = sin ( t ) + ( x3 ** 3.0e0 - x3 ) * cos ( x ) + 3.0e * x1 * x2 * sin ( x ) }
150 let { x = x + h * ( x1 + ( h / 2.0e0 ) * ( x2 + ( h / 6.0e0 ) * ( x3 + ( h / 24.0e0 ) * x4 ) ) ) }
160 let { t = t + h }
170 print k , { t , x }
200 next k
210 end ;
\s
\ Reference:
\ c
\ c Second Edition
\ c Numerical Analysis:
\ c The Mathematics of Scientific Computing
\ c D.R. Kincaid & E.W. Cheney
\ c Brooks/Cole Publ., 1996
\ c ISBN 0-534-33892-5
\ c COPYRIGHT (c) 1996
\ c
\ c Section 8.2
\ c
\ c Solving the initial value problem using Taylor Series
\ c
\ c
\ c file: taylor.f
\ c
\ data M,h,t,x/200,0.01,-1.0,3.0/
\ c
\ print *
\ print *,' Taylor series method (order 4) '
\ print *,' Section 8.1, Kincaid-Cheney'
\ print *
\ print 3,'k','t','x'
\ print 4,0,t,x
\ c
\ do 2 k=1,M
\ x1 = cos(t) - sin(x) + t**2.0
\ x2 = -sin(t) - x1*cos(x) + 2.0*t
\ x3 = -cos(t) - x2*cos(x) + (x1**2.0)*sin(x) + 2.0
\ x4 = sin(t) + ((x3)**3.0 -x3)*cos(x) + 3.0*x1*x2*sin(x)
\ x = x + h*(x1 + (h/2.)*(x2 + (h/6.)*(x3 + (h/24.)*x4)))
\ t = t + h
\ print 4,k,t,x
\ 2 continue
\ c
\ 3 format(a6,a9,a15)
\ 4 format(1x,i5,2x,e13.6,2x,e13.6)
\ stop
\ end
include taylor redefined k ok
taylor
*
* Taylor series method (order 4)
* Section 8.1, Kincaid-Cheney
*
0 -1.00000000000000E0 3.00000000000000E0
1 -9.90000000000000E-1 3.01400331867272E0
2 -9.80000000000000E-1 3.02803125373840E0
3 -9.70000000000000E-1 3.04208574304848E0
4 -9.60000000000000E-1 3.05616870908915E0
5 -9.50000000000000E-1 3.07028205863460E0
6 -9.40000000000000E-1 3.08442768236222E0
7 -9.30000000000000E-1 3.09860745442929E0
8 -9.20000000000000E-1 3.11282323201108E0
9 -9.10000000000000E-1 3.12707685479982E0
10 -9.00000000000000E-1 3.14137014446445E0
11 -8.90000000000000E-1 3.15570490407070E0
12 -8.80000000000000E-1 3.17008291746145E0
13 -8.70000000000000E-1 3.18450594859697E0
14 -8.60000000000000E-1 3.19897574085496E0
15 -8.50000000000000E-1 3.21349401629017E0
16 -8.40000000000000E-1 3.22806247485362E0
17 -8.30000000000000E-1 3.24268279357115E0
18 -8.20000000000000E-1 3.25735662568157E0
19 -8.10000000000000E-1 3.27208559973412E0
20 -8.00000000000000E-1 3.28687131864576E0
21 -7.90000000000000E-1 3.30171535871815E0
22 -7.80000000000000E-1 3.31661926861472E0
23 -7.70000000000000E-1 3.33158456829826E0
24 -7.60000000000000E-1 3.34661274792935E0
25 -7.50000000000000E-1 3.36170526672629E0
26 -7.40000000000000E-1 3.37686355178704E0
27 -7.30000000000000E-1 3.39208899687415E0
28 -7.20000000000000E-1 3.40738296116338E0
29 -7.10000000000000E-1 3.42274676795710E0
30 -7.00000000000000E-1 3.43818170336371E0
31 -6.90000000000000E-1 3.45368901494417E0
32 -6.80000000000000E-1 3.46926991032735E0
33 -6.70000000000000E-1 3.48492555579552E0
34 -6.60000000000000E-1 3.50065707484203E0
35 -6.50000000000000E-1 3.51646554670289E0
36 -6.40000000000000E-1 3.53235200486450E0
37 -6.30000000000000E-1 3.54831743554973E0
38 -6.20000000000000E-1 3.56436277618494E0
39 -6.10000000000000E-1 3.58048891385056E0
40 -6.00000000000000E-1 3.59669668371799E0
41 -5.90000000000000E-1 3.61298686747618E0
42 -5.80000000000000E-1 3.62936019175081E0
43 -5.70000000000000E-1 3.64581732651983E0
44 -5.60000000000000E-1 3.66235888352892E0
45 -5.50000000000000E-1 3.67898541471076E0
46 -5.40000000000000E-1 3.69569741061222E0
47 -5.30000000000000E-1 3.71249529883380E0
48 -5.20000000000000E-1 3.72937944248577E0
49 -5.10000000000000E-1 3.74635013866571E0
50 -5.00000000000000E-1 3.76340761696235E0
51 -4.90000000000000E-1 3.78055203799070E0
52 -4.80000000000000E-1 3.79778349196384E0
53 -4.70000000000000E-1 3.81510199730654E0
54 -4.60000000000000E-1 3.83250749931645E0
55 -4.49999999999999E-1 3.84999986887835E0
56 -4.40000000000000E-1 3.86757890123739E0
57 -4.30000000000000E-1 3.88524431483701E0
58 -4.19999999999999E-1 3.90299575022768E0
59 -4.09999999999999E-1 3.92083276905230E0
60 -3.99999999999999E-1 3.93875485311441E0
61 -3.89999999999999E-1 3.95676140353516E0
62 -3.79999999999999E-1 3.97485174000507E0
63 -3.69999999999999E-1 3.99302510013663E0
64 -3.59999999999999E-1 4.01128063892351E0
65 -3.49999999999999E-1 4.02961742831219E0
66 -3.39999999999999E-1 4.04803445689171E0
67 -3.29999999999999E-1 4.06653062970685E0
68 -3.19999999999999E-1 4.08510476820015E0
69 -3.09999999999999E-1 4.10375561028756E0
70 -2.99999999999999E-1 4.12248181057268E0
71 -2.89999999999999E-1 4.14128194070366E0
72 -2.79999999999999E-1 4.16015448987716E0
73 -2.69999999999999E-1 4.17909786549267E0
74 -2.59999999999999E-1 4.19811039396072E0
75 -2.49999999999999E-1 4.21719032166743E0
76 -2.39999999999999E-1 4.23633581609803E0
77 -2.29999999999999E-1 4.25554496712079E0
78 -2.19999999999999E-1 4.27481578843282E0
79 -2.09999999999999E-1 4.29414621916815E0
80 -1.99999999999999E-1 4.31353412566839E0
81 -1.89999999999999E-1 4.33297730341515E0
82 -1.79999999999999E-1 4.35247347912321E0
83 -1.69999999999999E-1 4.37202031299245E0
84 -1.59999999999999E-1 4.39161540111606E0
85 -1.49999999999999E-1 4.41125627804200E0
86 -1.39999999999999E-1 4.43094041948356E0
87 -1.29999999999999E-1 4.45066524517488E0
88 -1.19999999999999E-1 4.47042812186596E0
89 -1.09999999999999E-1 4.49022636645160E0
90 -9.99999999999992E-2 4.51005724922771E0
91 -8.99999999999993E-2 4.52991799726795E0
92 -7.99999999999993E-2 4.54980579791316E0
93 -6.99999999999993E-2 4.56971780236541E0
94 -5.99999999999993E-2 4.58965112937790E0
95 -4.99999999999993E-2 4.60960286903159E0
96 -3.99999999999992E-2 4.62957008658907E0
97 -2.99999999999992E-2 4.64954982641547E0
98 -1.99999999999992E-2 4.66953911595623E0
99 -9.99999999999925E-3 4.68953496976106E0
100 7.52869988573934E-16 4.70953439354318E0
101 1.00000000000008E-2 4.72953438826286E0
102 2.00000000000008E-2 4.74953195422398E0
103 3.00000000000008E-2 4.76952409517247E0
104 4.00000000000008E-2 4.78950782238533E0
105 5.00000000000008E-2 4.80948015873896E0
106 6.00000000000008E-2 4.82943814274579E0
107 7.00000000000008E-2 4.84937883254824E0
108 8.00000000000007E-2 4.86929930985933E0
109 9.00000000000007E-2 4.88919668383937E0
110 1.00000000000001E-1 4.90906809489885E0
111 1.10000000000001E-1 4.92891071841753E0
112 1.20000000000001E-1 4.94872176837072E0
113 1.30000000000001E-1 4.96849850085358E0
114 1.40000000000001E-1 4.98823821749537E0
115 1.50000000000001E-1 5.00793826875578E0
116 1.60000000000001E-1 5.02759605709597E0
117 1.70000000000001E-1 5.04720904001792E0
118 1.80000000000001E-1 5.06677473296597E0
119 1.90000000000001E-1 5.08629071208515E0
120 2.00000000000001E-1 5.10575461683161E0
121 2.10000000000001E-1 5.12516415243122E0
122 2.20000000000001E-1 5.14451709218274E0
123 2.30000000000001E-1 5.16381127960298E0
124 2.40000000000001E-1 5.18304463041180E0
125 2.50000000000001E-1 5.20221513435567E0
126 2.60000000000001E-1 5.22132085686877E0
127 2.70000000000001E-1 5.24035994057177E0
128 2.80000000000001E-1 5.25933060660847E0
129 2.90000000000001E-1 5.27823115582151E0
130 3.00000000000001E-1 5.29705996976865E0
131 3.10000000000001E-1 5.31581551158175E0
132 3.20000000000001E-1 5.33449632667108E0
133 3.30000000000001E-1 5.35310104327801E0
134 3.40000000000001E-1 5.37162837287958E0
135 3.50000000000001E-1 5.39007711044894E0
136 3.60000000000001E-1 5.40844613457587E0
137 3.70000000000001E-1 5.42673440745194E0
138 3.80000000000001E-1 5.44494097472537E0
139 3.90000000000001E-1 5.46306496523059E0
140 4.00000000000001E-1 5.48110559059800E0
141 4.10000000000001E-1 5.49906214474942E0
142 4.20000000000001E-1 5.51693400328501E0
143 4.30000000000001E-1 5.53472062276751E0
144 4.40000000000001E-1 5.55242153990973E0
145 4.50000000000001E-1 5.57003637067127E0
146 4.60000000000001E-1 5.58756480927059E0
147 4.70000000000001E-1 5.60500662711828E0
148 4.80000000000001E-1 5.62236167167787E0
149 4.90000000000001E-1 5.63962986525971E0
150 5.00000000000001E-1 5.65681120375425E0
151 5.10000000000001E-1 5.67390575531022E0
152 5.20000000000001E-1 5.69091365896359E0
153 5.30000000000001E-1 5.70783512322274E0
154 5.40000000000001E-1 5.72467042461543E0
155 5.50000000000001E-1 5.74141990620275E0
156 5.60000000000001E-1 5.75808397606512E0
157 5.70000000000001E-1 5.77466310576543E0
158 5.80000000000001E-1 5.79115782879393E0
159 5.90000000000001E-1 5.80756873899947E0
160 6.00000000000001E-1 5.82389648901151E0
161 6.10000000000001E-1 5.84014178865698E0
162 6.20000000000001E-1 5.85630540337598E0
163 6.30000000000001E-1 5.87238815264011E0
164 6.40000000000001E-1 5.88839090837694E0
165 6.50000000000001E-1 5.90431459340389E0
166 6.60000000000001E-1 5.92016017987487E0
167 6.70000000000001E-1 5.93592868774237E0
168 6.80000000000001E-1 5.95162118323788E0
169 6.90000000000001E-1 5.96723877737319E0
170 7.00000000000001E-1 5.98278262446484E0
171 7.10000000000001E-1 5.99825392068402E0
172 7.20000000000001E-1 6.01365390263375E0
173 7.30000000000001E-1 6.02898384595535E0
174 7.40000000000001E-1 6.04424506396569E0
175 7.50000000000001E-1 6.05943890632684E0
176 7.60000000000001E-1 6.07456675774938E0
177 7.70000000000001E-1 6.08963003673064E0
178 7.80000000000001E-1 6.10463019432890E0
179 7.90000000000001E-1 6.11956871297446E0
180 8.00000000000001E-1 6.13444710531846E0
181 8.10000000000001E-1 6.14926691312004E0
182 8.20000000000001E-1 6.16402970617256E0
183 8.30000000000001E-1 6.17873708126922E0
184 8.40000000000001E-1 6.19339066120860E0
185 8.50000000000001E-1 6.20799209384029E0
186 8.60000000000001E-1 6.22254305115098E0
187 8.70000000000001E-1 6.23704522839102E0
188 8.80000000000001E-1 6.25150034324157E0
189 8.90000000000001E-1 6.26591013502245E0
190 9.00000000000001E-1 6.28027636394051E0
191 9.10000000000001E-1 6.29460081037847E0
192 9.20000000000001E-1 6.30888527422417E0
193 9.30000000000001E-1 6.32313157424002E0
194 9.40000000000001E-1 6.33734154747233E0
195 9.50000000000001E-1 6.35151704870052E0
196 9.60000000000001E-1 6.36565994992567E0
197 9.70000000000001E-1 6.37977213989839E0
198 9.80000000000001E-1 6.39385552368558E0
199 9.90000000000001E-1 6.40791202227567E0
200 1.00000000000000E0 6.42194357222237E0 ok
第(89)個範例是兩套設計出常用但非標準指令 within ( n1 n2 n3 — flag ) 的方法,設計方式帶有技巧。
使用 Forth 系統應該養成一點寫程式時的好習慣,以免日後回頭追蹤時,給自己帶來困擾。這個習慣,就是每當設計任何新指令時,要記得在指令名稱的同一列內,加添表示這個指令執行前後,堆疊數據的變化狀況,例如:此例為 ( n1 n2 n3 — flag )。
within 的功能是判斷出 n1 是不是介於 [n2, n3) 之間?這是標準嚴謹的階梯函數之數學表示法,寫成數學式則為:
n2 =< n1 < n3
亦即,n1 屬於 [n2, n3) 之間的所有數字,包括 n2 但不包括 n3 ,這樣的規格,可以從微積分課本第一章中獲得。
程式執行結果,則會得到一個旗號,為真時得 -1 ,為假時則得 0 。
一般而言,常人在得到上述問題時的直接反應,會是去採用邏輯比較法來設計程式,也就是採用假如怎樣就會怎樣的思考方式來設計程式,程式中就會出現許多 IF ..... ELSE ..... THEN 的結構,這個範例程式中卻完全沒有。
我經常瀏覽國際論壇上的討論,主要目的就是去發掘一些名家的設計,可以發現許多採用非常規思考方式所設計出來的漂亮結果,這個範例就是在這種情況下所獲得的。
設計者完全不用邏輯比較結構,卻能設計得出這個必須至少得用兩次比較才能產生正確輸出結果的程式。至少得用兩次的比較,被採用一次減法運算,與另一個無號數的比較指令 u< ,實現出了等效的結果。
能設計出這種技術程式者,算是高手,簡化了程式量,也讓程式執行速度大量加速,很值得欽佩。
這幾十年來,我的個人筆記已經累積了 20 幾本,本本都記錄了這類的成就,筆記中除了自己的思考記錄,其餘的,就都是名家設計的精華。我每次為新系統發展自創部份時,都會拿出來參考,好幾本筆記本都已用到脫頁了,表示它們長期貢獻過設計工作,我很珍惜。
這兩套設計方法都很漂亮, within1 採用將過渡期的暫態數據暫存於回返堆疊的方法來設計程式,我後來覺得有些新式的 Forth 系統,把回返堆疊設計得很不像樣,還猛用,總有一天會出問題,所以,我還是按照 Forth 的傳統精神,改寫出規規矩矩全只使用堆疊指令的方式,設計出了 within2 ,這就是這兩個功效一樣的程式設計上之差異。
想要測試指令也很簡單,這種事我常做,在此提醒大家,處理標準數學問題時,千萬要記得『數字有正負,範圍有邊界』,測試程式時,它們都是經常會被忽略的問題,務必永遠記得,此範例都能測試得到的ㄧ序列測試如下:
-7 -5 5 within .
-5 -5 5 within .
0 -5 5 within .
5 -5 5 within .
7 -5 5 within .
如果您將上列測試執行完畢後,所得 0 或 -1 的旗號結果,都符合要求,那麼,這樣的設計,就算是絕對沒有問題了。
:
\ (89)high level within
: within1 ( n1 n2 n3 -- flag )
over - >r - r> u< ;
: within2
OVER - -ROT - SWAP U< ;
第(90)個範例是有關陣列(array)、矩陣(matrix)及張量(tensor)資料結構的測試程式。
由於這三種資料結構又可分為三種不同的內容:整數、實數、複數。結構內容也會不同,所以都必須測試,測試的環境,是 64 位元的 gForth 系統,測試的方法,是使用指定數字與亂數置入之後,再按秩序整齊的列印出來,印出的同時,就根據各種維度的指標,來換算出應該取得的資料結構之位址,取出內容,直接印出。
這個程式不能在 Win32Forth 系統中使用,主要的影響,是整數指標不能在實數與複數的運作環境中使用。
這是一個很麻煩而難以處理的問題,要使用令系統在處理當時能忽略掉指標為整數而不是實數或複數的特殊技術。該忽略幾次?隨維度為多少而變化。還有,有些 Forth 系統裡,恆把整數 -1 , 0 , 1 , 2 都設計成常數,以便提高執行速度。在我自創的 ABC Forth 系統中,這些問題都造成很大的困擾,所以設計出來的應用規則,就會有很大的差異。
上述涉及的問題,在不同的系統中,解決問題的方法都不一樣,我沒有辦法把它們一一列示於此。關於這方面的技術,我就不談了。
我在發展有維度的資料結構設計時,進展到矩陣階段,發現原在 Win32Forth 中的設計方式不太理想,原始的設計格式是採用了 Forth 鼻祖 Charles H. Moore 的規格,他很省記憶體,又喜歡凡事都倒過來用,結果,對矩陣二維度以上的東西,用起來就會有麻煩。
矩陣在數學上有其獨特的運算,如點乘(dot production)或叉乘(cross production),運算時,都得憑維度指標來取得元素位址,才能有正確的數字可用。如果維度在資料結構中存放不足或秩序顛倒,設計起來,就有很大的麻煩。我從 32 位元環境走進 64 位元環境後,就完全放棄了原來的設計規格,改採自己在受盡折磨後才體會出來的方法,修正設計。這也是為什麼我在設計出 gForth64 的創作後,仍須重新測試這批資料結構的原因。
:
\ (90)A M T test in gforth64 ABC.20171202
4 (array) aa
3 4 (matrix) mm
2 2 2 (tensor) tt
3 integers i j k
: initaa basic
10 for i = 0 to 4
20 let aa ( i ) = i
30 next i
40 end ;
: initmm basic
10 for i = 0 to 3
20 for j = 0 to 4
30 let mm ( i j ) = i * j
40 next j
50 next i
60 end ;
: inittt basic
10 for i = 0 to 2
20 for j = 0 to 2
30 for k = 0 to 2
40 let tt ( i j k ) = i * j * k
50 next k
60 next j
70 next i
80 end ;
: .aa basic
10 for i = 0 to 4
20 print " aa( " ; i ; " )= " , aa ( i )
30 next i
40 end ;
: .mm basic
10 for i = 0 to 3
20 for j = 0 to 4
30 print " mm( " ; i ; j ; " )= " , mm ( i j )
40 next j
50 next i
60 end ;
: .tt basic
10 for i = 0 to 2
20 for j = 0 to 2
30 for k = 0 to 2
40 print " tt(" ; i ; j ; k ; " ) = " , tt ( i j k )
50 next k
60 next j
70 next i
80 end ;
4 array faa
3 4 matrix fmm
2 2 2 tensor ftt
: initfaa basic
10 for i = 0 to 4
20 let { faa ( i ) = I>R ( i ) }
30 next i
40 end ;
: initfmm basic
10 for i = 0 to 3
20 for j = 0 to 4
30 let { fmm ( i j ) = I>R ( i ) * I>R ( j ) }
40 next j
50 next i
60 end ;
: initftt basic
10 for i = 0 to 2
20 for j = 0 to 2
30 for k = 0 to 2
40 let { ftt ( i j k ) = I>R ( i ) * I>R ( j ) * I>R ( k ) }
50 next k
60 next j
70 next i
80 end ;
: .faa basic
10 for i = 0 to 4
20 print " faa( " ; i ; " )= " , { faa ( i ) }
30 next i
40 end ;
: .fmm basic
10 for i = 0 to 3
20 for j = 0 to 4
30 print " fmm( " ; i ; j ; " )= " , { fmm ( i j ) }
40 next j
50 next i
60 end ;
: .ftt basic
10 for i = 0 to 2
20 for j = 0 to 2
30 for k = 0 to 2
40 print " ftt( " ; i ; j ; k ; " ) = " , { ftt ( i j k ) }
50 next k
60 next j
70 next i
80 end ;
4 [array] zaa
3 4 [matrix] zmm
2 2 2 [tensor] ztt
: initzaa basic
10 for i = 0 to 4
20 let Z{ zaa ( i ) = 1.0e0 - 2.0e0 ii }Z
30 next i
40 end ;
: initzmm basic
10 for i = 0 to 3
20 for j = 0 to 4
30 let Z{ zmm ( i j ) = ZRAND }Z
40 next j
50 next i
60 end ;
: initztt basic
10 for i = 0 to 2
20 for j = 0 to 2
30 for k = 0 to 2
40 let z{ ztt ( i j k ) = ZRAND }Z
50 next k
60 next j
70 next i
80 end ;
: .zaa basic
10 for i = 0 to 4
20 print " zaa( " ; i ; " )= " , Z{ zaa ( i ) }Z
30 next i
40 end ;
: .zmm basic
10 for i = 0 to 3
20 for j = 0 to 4
30 print " zmm( " ; i ; j ; " )= " , Z{ zmm ( i j ) }Z
40 next j
50 next i
60 end ;
: .ztt basic
10 for i = 0 to 2
20 for j = 0 to 2
30 for k = 0 to 2
40 print " ztt( " ; i ; j ; k ; " ) = " , Z{ ztt ( i j k ) }Z
50 next k
60 next j
70 next i
80 end ;
: main basic
10 run initaa
20 run initmm
40 run inittt
50 run .aa cr
60 run .mm cr
70 run .tt cr
110 run initfaa
120 run initfmm
130 run initftt
140 run .faa cr
150 run .fmm cr
160 run .ftt cr
210 run initzaa
220 run initzmm
230 run initztt
240 run .zaa cr
250 run .zmm cr
260 run .ztt cr
300 end ;
main
\s
ching@center:~$ cd gforth
ching@center:~/gforth$ ./abc06
Gforth 0.7.2, Copyright (C) 1995-2008 Free Software Foundation, Inc.
Gforth comes with ABSOLUTELY NO WARRANTY; for details type `license'
Type `bye' to exit
include testamt redefined i redefined j redefined k ok
main
aa( 0 )= 0
aa( 1 )= 1
aa( 2 )= 2
aa( 3 )= 3
aa( 4 )= 4
mm( 0 0 )= 0
mm( 0 1 )= 0
mm( 0 2 )= 0
mm( 0 3 )= 0
mm( 0 4 )= 0
mm( 1 0 )= 0
mm( 1 1 )= 1
mm( 1 2 )= 2
mm( 1 3 )= 3
mm( 1 4 )= 4
mm( 2 0 )= 0
mm( 2 1 )= 2
mm( 2 2 )= 4
mm( 2 3 )= 6
mm( 2 4 )= 8
mm( 3 0 )= 0
mm( 3 1 )= 3
mm( 3 2 )= 6
mm( 3 3 )= 9
mm( 3 4 )= 12
tt(0 0 0 ) = 0
tt(0 0 1 ) = 0
tt(0 0 2 ) = 0
tt(0 1 0 ) = 0
tt(0 1 1 ) = 0
tt(0 1 2 ) = 0
tt(0 2 0 ) = 0
tt(0 2 1 ) = 0
tt(0 2 2 ) = 0
tt(1 0 0 ) = 0
tt(1 0 1 ) = 0
tt(1 0 2 ) = 0
tt(1 1 0 ) = 0
tt(1 1 1 ) = 1
tt(1 1 2 ) = 2
tt(1 2 0 ) = 0
tt(1 2 1 ) = 2
tt(1 2 2 ) = 4
tt(2 0 0 ) = 0
tt(2 0 1 ) = 0
tt(2 0 2 ) = 0
tt(2 1 0 ) = 0
tt(2 1 1 ) = 2
tt(2 1 2 ) = 4
tt(2 2 0 ) = 0
tt(2 2 1 ) = 4
tt(2 2 2 ) = 8
faa( 0 )= 0.00000000000000E0
faa( 1 )= 1.00000000000000E0
faa( 2 )= 2.00000000000000E0
faa( 3 )= 3.00000000000000E0
faa( 4 )= 4.00000000000000E0
fmm( 0 0 )= 0.00000000000000E0
fmm( 0 1 )= 0.00000000000000E0
fmm( 0 2 )= 0.00000000000000E0
fmm( 0 3 )= 0.00000000000000E0
fmm( 0 4 )= 0.00000000000000E0
fmm( 1 0 )= 0.00000000000000E0
fmm( 1 1 )= 1.00000000000000E0
fmm( 1 2 )= 2.00000000000000E0
fmm( 1 3 )= 3.00000000000000E0
fmm( 1 4 )= 4.00000000000000E0
fmm( 2 0 )= 0.00000000000000E0
fmm( 2 1 )= 2.00000000000000E0
fmm( 2 2 )= 4.00000000000000E0
fmm( 2 3 )= 6.00000000000000E0
fmm( 2 4 )= 8.00000000000000E0
fmm( 3 0 )= 0.00000000000000E0
fmm( 3 1 )= 3.00000000000000E0
fmm( 3 2 )= 6.00000000000000E0
fmm( 3 3 )= 9.00000000000000E0
fmm( 3 4 )= 1.20000000000000E1
ftt( 0 0 0 ) = 0.00000000000000E0
ftt( 0 0 1 ) = 0.00000000000000E0
ftt( 0 0 2 ) = 0.00000000000000E0
ftt( 0 1 0 ) = 0.00000000000000E0
ftt( 0 1 1 ) = 0.00000000000000E0
ftt( 0 1 2 ) = 0.00000000000000E0
ftt( 0 2 0 ) = 0.00000000000000E0
ftt( 0 2 1 ) = 0.00000000000000E0
ftt( 0 2 2 ) = 0.00000000000000E0
ftt( 1 0 0 ) = 0.00000000000000E0
ftt( 1 0 1 ) = 0.00000000000000E0
ftt( 1 0 2 ) = 0.00000000000000E0
ftt( 1 1 0 ) = 0.00000000000000E0
ftt( 1 1 1 ) = 1.00000000000000E0
ftt( 1 1 2 ) = 2.00000000000000E0
ftt( 1 2 0 ) = 0.00000000000000E0
ftt( 1 2 1 ) = 2.00000000000000E0
ftt( 1 2 2 ) = 4.00000000000000E0
ftt( 2 0 0 ) = 0.00000000000000E0
ftt( 2 0 1 ) = 0.00000000000000E0
ftt( 2 0 2 ) = 0.00000000000000E0
ftt( 2 1 0 ) = 0.00000000000000E0
ftt( 2 1 1 ) = 2.00000000000000E0
ftt( 2 1 2 ) = 4.00000000000000E0
ftt( 2 2 0 ) = 0.00000000000000E0
ftt( 2 2 1 ) = 4.00000000000000E0
ftt( 2 2 2 ) = 8.00000000000000E0
zaa( 0 )= 1.00000000000000E0 - 2.00000000000000E0 ii
zaa( 1 )= 1.00000000000000E0 - 2.00000000000000E0 ii
zaa( 2 )= 1.00000000000000E0 - 2.00000000000000E0 ii
zaa( 3 )= 1.00000000000000E0 - 2.00000000000000E0 ii
zaa( 4 )= 1.00000000000000E0 - 2.00000000000000E0 ii
zmm( 0 0 )= 2.18418296993905E-1 + 9.56317576559408E-1 ii
zmm( 0 1 )= 8.29509233976486E-1 + 5.61695442796543E-1 ii
zmm( 0 2 )= 4.15307081497883E-1 + 6.61187349195214E-2 ii
zmm( 0 3 )= 2.57577792395641E-1 + 1.09956793538275E-1 ii
zmm( 0 4 )= 4.38289977814206E-2 + 6.33965712335876E-1 ii
zmm( 1 0 )= 6.17272290688600E-2 + 4.49538960330905E-1 ii
zmm( 1 1 )= 4.01306281518799E-1 + 7.54673486461245E-1 ii
zmm( 1 2 )= 7.97286954148340E-1 + 1.83837115850224E-3 ii
zmm( 1 3 )= 8.97504060947105E-1 + 3.50752337999061E-1 ii
zmm( 1 4 )= 9.45447502166707E-2 + 1.36168915841807E-2 ii
zmm( 2 0 )= 8.59096855325204E-1 + 8.40847450700052E-1 ii
zmm( 2 1 )= 1.23103915771052E-1 + 7.51236407436075E-3 ii
zmm( 2 2 )= 2.60302997781105E-1 + 9.12483707029598E-1 ii
zmm( 2 3 )= 1.13664046448499E-1 + 3.51628659922457E-1 ii
zmm( 2 4 )= 8.22887316729355E-1 + 2.67132270274280E-1 ii
zmm( 3 0 )= 6.92066499820010E-1 + 5.61662474908709E-1 ii
zmm( 3 1 )= 8.61215790669069E-1 + 4.53793775035904E-1 ii
zmm( 3 2 )= 9.11977028433223E-1 + 5.97916877175643E-1 ii
zmm( 3 3 )= 1.88954691024942E-1 + 7.61492056195388E-1 ii
zmm( 3 4 )= 3.96988475880115E-1 + 1.85314117085801E-1 ii
ztt( 0 0 0 ) = 5.74365861050024E-1 + 3.67026667747193E-1 ii
ztt( 0 0 1 ) = 6.17204827078248E-1 + 3.61528704111245E-1 ii
ztt( 0 0 2 ) = 2.12929997692318E-1 + 7.14471214783597E-1 ii
ztt( 0 1 0 ) = 1.17706867921030E-1 + 2.99329148744852E-1 ii
ztt( 0 1 1 ) = 8.25002954725643E-1 + 8.24660073884139E-1 ii
ztt( 0 1 2 ) = 6.18617707220194E-2 + 7.10780524979709E-1 ii
ztt( 0 2 0 ) = 8.82833339685031E-2 + 7.77994008631443E-1 ii
ztt( 0 2 1 ) = 7.45303068657081E-1 + 3.08674919562729E-1 ii
ztt( 0 2 2 ) = 8.99373090779117E-1 + 7.63536724617489E-1 ii
ztt( 1 0 0 ) = 7.61730646137954E-1 + 4.06969640593496E-1 ii
ztt( 1 0 1 ) = 9.38749454887002E-1 + 5.62088285834570E-1 ii
ztt( 1 0 2 ) = 1.78200216115546E-2 + 5.01103225397460E-1 ii
ztt( 1 1 0 ) = 4.19092551068912E-2 + 3.68850581519702E-1 ii
ztt( 1 1 1 ) = 2.71723601627966E-1 + 8.58572561227983E-1 ii
ztt( 1 1 2 ) = 2.90365587123840E-2 + 1.74422790377598E-2 ii
ztt( 1 2 0 ) = 1.52383787628442E-1 + 1.14318671223856E-1 ii
ztt( 1 2 1 ) = 3.53907259345943E-1 + 1.19307827260023E-1 ii
ztt( 1 2 2 ) = 2.06652759204923E-1 + 2.12923957134096E-1 ii
ztt( 2 0 0 ) = 6.12947552750328E-1 + 8.09519074768535E-1 ii
ztt( 2 0 1 ) = 5.87089634773829E-1 + 2.15491643741490E-1 ii
ztt( 2 0 2 ) = 7.68056363225010E-1 + 7.23296722734020E-1 ii
ztt( 2 1 0 ) = 4.48018990665683E-1 + 8.55176118135069E-1 ii
ztt( 2 1 1 ) = 9.45017496098307E-1 + 9.09056924241156E-1 ii
ztt( 2 1 2 ) = 5.19725721105806E-1 + 3.01946252725062E-2 ii
ztt( 2 2 0 ) = 4.81066955011835E-1 + 2.92312883908075E-1 ii
ztt( 2 2 1 ) = 9.02639843012504E-1 + 6.67841511158199E-1 ii
ztt( 2 2 2 ) = 4.12278035847600E-1 + 1.56948490607063E-1 ii
ok