2024年12月15日 星期日

知難行易

知難行易


曾慶潭 Ching-Tang Tseng
ilikeforth@gmail.com
Hamilton, New Zealand
16 December 2024


這是    國父孫中山先生的學說,我們小的時候使用過的教科書中,有過這麼一課,課文的標題就是:『知難行易』。

說實在話,當時讀到這一課時,我並不能深入了解知難行易的真正意義。因為課文的內容,舉的例子很簡單,是描述    國父旅居海外時,有一天找了水匠來家裡修理水龍頭,只見水匠拿個工具扳手扭一扭,就把水龍頭的問題解決了。課文用這樣的例子說明:要想知道水龍頭問題的道理並不容易,不是常人輕易就能辦到的事情。可是 國父親眼目睹水匠只用工具扭了一扭就解決了問題,意思是實際的修理工作,做起來就非常容易。

我移民海外也將近三十年了,知道紐西蘭的水匠是一種辛苦、骯髒的高收入職業,他們就算在某家公司或行號上班,整車的工具也都是自備的。請水匠登門修理東西,通常都得付出高昂的代價,從出車開始計時、計費,工時按每半小時算作一個計費單元來算,耗用的工料最多只佔總費用的三分之一,修理工錢大約也要三分之一。換句話說,他們能賺得的利潤,就是剩下的三分之一。還有,紐西蘭的水匠與電匠是分開的,不像台灣的水電工,水與電的修理工作都可以同時做。

我在這裡度過了將近三十年的居家生活,但從未找過水匠來家裡修過東西,發揚了紐西蘭人凡事喜歡 DIY(Do It Yourself) 的精神,家庭用水設施與裝置上的任何問題我都自行解決。家裡水龍頭漏水、滴水、關不死是經常發生的事情,不尋求水匠來修理,那就得自己親自動手解決問題。因此可以體會知難行易課文的意義。

我住在紐西蘭北島中央一個不濱海的 Hamilton 市區內,本市民生供水所需要的經費,只從收取的年度房屋稅中分配來使用。居民用水不必另行計費的狀況下,家裡的水龍頭漏水、滴水時,大家就都不太當它是一回事,通常都不急著要找水匠來及時修理。我自己有工具,又熟悉各種水系統的運行原理,能自行解決問題。

知難行易所舉之例,是一般住家內最簡單的水龍頭修理工作,可是    國父就是不懂、不知道如何解決,表示常人大多數也不太理解這方面的問題。我在單純的修理這種水龍頭問題時,體會出一些必須經常身體力行才會注意得到的要點,如果在修理水龍頭時不注意這些事情,只是隨便轉緊所有的螺絲,有可能修過的水龍頭都照樣還會繼續漏水或滴水。

網上有許多影片教您如何修理漏水的水龍頭,但我沒見過有人能將安裝水龍頭時的工作要點強調得出來的影片,討論知難行易的道理,講一講水龍頭,就是有意義的事情。

讓我們先從水龍頭的內部結構開始談起。我所借用的圖示,取材自一個收費網頁,因為可以免費下載 15 張例圖,這是我下載的第一張圖,所以沒有計費,我們使用它時,仍應該表示感謝。網頁是:

https://www.dreamstime.com/stock-illustration-cross-section-brass-faucet-isolated-white-background-d-illustration-image84036492

我取得圖示後,自行添加了必要配件的中文名稱,圖示如下:


水龍頭轉桿的向上運動,會將壓緊在閥座上的桿底墊片向上提昇,讓自來水流過閥座噴出水龍頭的噴口。轉桿向下運動到底,則將桿底墊片壓緊在閥座上,水龍頭便停止供水。

假設水龍頭所有的配件都沒有壞,安裝這樣的水龍頭時,操作步驟有個要領,照著步驟做,才能裝得成能夠理想工作的水龍頭。這個要領就是安裝時,首要重視轉桿所處的位置,轉桿應該是轉到下死點的位置,但是絕對不可以在下死點處轉得太死,理想的位置應該是輕鬆轉到下死點後再往回轉回半圈,讓轉桿成為鬆而不緊的狀態。這個預處理效果,是為了避免在安裝固定螺絲時,整體會因轉桿同時被壓死在閥座上而妨礙了讓固定螺絲繼續緊定於必須絕對封死的位置。這個回轉半圈的動作,在水龍頭正常使用時,也是桿底墊片彈性工作間隙密閉時所需要的壓緊空間。

上述要領,也就是修理水龍頭更換小配件後,重新安裝時必須注意的要領。固定螺絲要轉得越緊越好,才能避免漏水,為了轉緊,工具可能滑脫,會傷了水龍頭拋光後的漂亮表面,老手都會另外安排不易滑脫的套筒,例如套上紙管或布來保護,但這都不是要點,只影響美觀,或避免手受傷害。其他的小螺絲,例如:在最下方緊定桿底墊片的小螺絲,在最上方緊定轉柄的小螺絲,或者是為了美觀而在最外面另外加裝的匣蓋盒子螺絲,只要裝緊了便可,都不強調必須轉得很緊,以防下次修理或年久銹蝕之後難以拆卸。

同樣的道理也經常被使用於原子爐系統的操作中,大型閥門關死時,當然要關得夠緊。大型閥門被開啟時,就要避免在到達全開位置後,閥門轉盤繼續被強轉到非常緊鎖的狀態,否則下次進行關閉閥門的操作時,會發生轉盤難以轉動的問題。正確的操作方式,是確定轉到全開位置後,轉盤回轉半圈,讓下次要關閉操作轉盤時,容易轉得動。

駕駛汽車打方向盤時也是這樣,要避免把方向盤打死,左轉或右轉確定方向盤轉到底了,方向盤就應該立即回轉半圈來控制轉向系統,否則方向機的油壓系統容易受損,高壓油可能會擠爆油封而開始漏油。

這些轉到極限點立刻回轉半圈的訣竅,能確保被操作之裝置使用起來比較安全的道理,都是其來有自的。

所有的水龍頭,無論它後來的發展出現了任何造型,只要是旋轉式的,其內部結構之基本運作機制,跟上列的圖示都是一樣的。這張展示圖最為精簡、清楚,弄懂了這個知難行易所舉之例的難懂東西,有利於自己動手修理水龍頭。碰到新的產品,也可以舉一反三的方式照樣推理以解決問題。


我們使用現代的任何 Forth 系統,也會遇到許多類似知難行易的情況,舉個簡單例子來說,例如:想把數字印出來的指令『 . 』(唸作 dot ) ,在一般流行性程式語言中常被定名為 print ,其使用方式對任何程式語言來說都是有點複雜的。早期電腦還是低位元的時代,我們一直強調 Forth 的透通性,那個時候,想要追蹤與硬體有關的訊息輸出指令之內容,需要了解的資料並不太多,直讀系統程式的源碼,很快就能明白訊息都是如何被印出來並被顯示在螢幕上的。後來,環境變了,作業系統成為電腦中的重要角色,程式語言都得依它而建,Forth 系統百分之百的透通性就無法再被強調,這也對我們使用 Forth 系統時,形成了『知難』 的障礙。不過,我們還是有辦法將『知難』化解出『行易』的,就像搞懂水龍頭內部的關鍵結構後,當然可以憑經驗培養出修理水龍頭時該有的工作素養,修水龍頭是這樣,順利使用各種 Forth 系統也是這樣。

回頭檢視 20081129 最原始的網文:FORTH-83 標準指令。你可以發現標準指令中只有兩個指令直接執行數字輸出,一個就是 . ,另一個就是 u. 。數字有很多格式,除了整數、無號數,還有雙整數、無號雙整數、大數、浮點數、複數 ..... 等等,浮點數又可再分為帶 10 的幾次方的科學數字表示格式、不帶方次的格式、工程上愛用的一百以內之兩位整數表現格式 ..... 等等等,格式雖然複雜,用起來則必須採用『行易』的觀念,由兩個基本數字印出指令而延伸出來之其他的數字輸出指令,都不再被列為標準指令,因為他們的使用規格都是一致的,不需要另行解釋。

此外,我個人有些系統設計經驗,體會出所有的數字輸出指令,都應該先行備妥其核心結構指令後再用比較好。數字輸出指令係直接將堆疊上的數字輸出,先行準備之相應的核心結構指令則將堆疊上的數字先行轉換到某一記憶體指定區,然後轉為提供記憶體指定區的起始位址(addr),另外也提供待印數字的字長(len),這兩個參數可以讓後續打算設計的數字輸出指令更具有使用上的彈性。尤其是在當數字輸出不僅只是送往螢幕這麼簡單的地方而已時,例如:是一個物件或一個檔案時,後續的設計就容易用它來解決問題。這類核心指令,按照公稱的使用習慣,我們都以其相應的指令外加一組小括號來表示,對應關系列示如下:


 
它們之間的關係則簡單到形同這個樣子:
: . ( n -- ) (.) type space ;

整數 . ( n -- ) (.) ( n -- addr len )
無號數 u. ( u -- ) (u.) ( u -- addr len )
雙整數 d. ( d -- ) (d.) ( d -- addr len )
無號雙整數 ud. ( ud -- ) (ud.) ( ud -- addr len )
浮點數 f. ( f -- ) (f.) ( f -- addr len )
科學表示的浮點數 fe. ( f -- ) (fe.) ( f -- addr len )
複數 z. ( f1 f2 -- ) (z.) ( f1 f2 -- addr len )


根據經驗培養出來的工作素養,深入剖析知難行易的道理,告訴我修理水龍頭時應該有什麼重點要求。
同樣的,設計系統很久之後,深入剖析知難行易的道理,也告訴我在設計各種數字輸出指令時,有什麼重點要求該預先設計。

2024年12月1日 星期日

盍各言爾志

盍各言爾志


曾慶潭 Ching-Tang Tseng
ilikeforth@gmail.com
Hamilton, New Zealand
1 December 2024


1960 年代,是我們這群 70 幾歲的老人開始上初中的年代。最近,我剛好碰上了一個外孫女初中畢業,要準備上高中了。她經常來我這裡陪我們生活,有過許多的隔代交流,每次回家時,我女婿就會問她今天又跟阿公學了什麼?
我能教她什麼?都是隨興的,不拘形式。

我記得她剛上初中時,我告訴過她三項上初中時會與上小學時不同的心情感受:
一是初中生會覺得老師在生活方面不再怎麼細微的照顧學生了。
二是你上中學時不再那麼容易的得到各種鼓勵性質的獎狀。
三是你會發現從別的小學畢業而來的同學,出現了好幾個能力都比你優秀的同學。
她必需經過比較後了解這三種情況,上學時,才不容易產生挫折感。
現在,她初中畢業了,一切狀況都還很優秀,表示歷經小學轉上初中的學業情況,與學習上的心理狀況都保持的很不錯。

我當然希望她再接再厲,繼續以平衡的心態接受高中教育。這個時候,我也還能再給她一點鼓勵,我能準備的說明題材就是教她如何立定志向?順便告訴她一點中國古文化在教育中被使用、推演的方式。
我告訴她,學生上了初中、高中,是應該立定自己志向的時候了。但我不強調對未來職業的選擇該是那一種,而是要強調在個人的修養與品行上應該立定什麼樣的志向。順便也告訴她,我當年上中學時的情況,跟現在的學生並沒有太大的差別。只不過,我們是幸運的,那時的台灣,有著許多剛從中國大陸保留來台的精英學者,他們認認真真的成立過中華文化復興委員會,教育部頒發的國語、國文課本,內容都是經過這些精英挑選過的,我們受教時立刻受益。現在的台灣已經不是這樣了,政治掛帥的結果是教材內容亂七八糟,沒有人恭維。

初中也是我們剛開始要正式接受中國古文化教育的時期,一定要開始接觸古書中的文言文,也就是文章言論方式的文字,有別於常人說話時講話方式的語體文文字。最典型之文言文的第一課,就是孔子與學生對談交流『志向』的課文。這一課,學生都要背誦,終生牢記,經常引用,身體力行。課文的內容在『論語』書中第五章第 26 條,是這樣言簡意賅的記錄著:

顏淵、季路侍。子曰:「盍各言爾志?」子路曰:「願車馬、衣輕裘,與朋友共,蔽之而無憾。」顏淵曰:「願無伐善,無施勞。」子路曰:「願聞子之志。」子曰:「老者安之,朋友信之,少者懷之。」

剛從小學畢業讀過國語、只強調學會認字的學生,並不能直接看懂這樣的文言文,拿它當教材就是一篇進行中國古文教育時很好的啟蒙材料。

我當講故事般的方式講給外孫女聽。2500 年前,孔子與兩位傑出的學生談論志向,一個是學生中最有錢的叫子路,另一個是學生中最窮的叫顏回。

孔子對他們說:『何不各自談一談你們的志向呢?』子路先回答說:『我希望我能作到我的車、馬、昂貴的輕質衣裳,與朋友分享使用時,用壞了也不會覺得遺憾而抱怨。』顏回則說:『我希望我不會吹噓自己的能力,也不會誇大自己貢獻過工作後的功勞。』隨後子路恭敬的說:『請問老師的志向為何?』孔子回答說:『我希望老年人能安享天年,朋友都信任我的為人,少小的人都能受到關懷照顧。』

課文中,盍各言爾志、蔽之而無憾、願聞子之志.....等等句子都是標準的文言文。句子中帶有『之』 字的用法者,就是標準文言虛字『之』字的用法,能美化句子,加強語氣,強調重點。

文章凸顯了 2500 年前生活富裕的現象,跟現在一樣,有馬、有車、服飾輕柔者代表的是有錢人生活時使用的物質。中國的絲綢早在 4800 年前就已有文字記載,孔子的時代也正是西方國家的羅馬時代,那時中國的絲綢,經絲路貿易銷往羅馬時,羅馬人要用黃金來購買,可見穿絲綢之『輕裘』者是多麼的富有。

顏回雖窮,他 32 歲就去世了,卻是孔子最常表揚的學生,說他好學、不遷怒於無辜者、從不犯下第二次的錯誤。他家雖窮,住在貧民區內、日常生活的飲食都有困難,但他仍然不改變安貧樂道的生活態度,算得上是個賢能的人。

我們在接受論語教育時,並不需要每章、每節、逐條的學習,許多書上記錄之孔子的言行,都是我們在讀完高中、有了文言文閱讀能力的基礎後,自行去找整本論語來自己看懂的,所以能夠知道孔子是如何批評顏回是個賢人的依據。

每個人年輕時都該有個自己的志向,你可以稱它是個人的座右銘,拿它當個人發展時的方向指標。從孔子教育學生的言論中,你也可以看出一個偉人立定志向時寬大的心胸,兩個傑出的學生志向,偏重於修養自己,都沒有野心。孔子就不同了,他的志向考慮到社會群體,令老人、朋友、小孩都能受惠,他的志向與眾不同。世人就憑他對整個世界的胸懷,確認他是聖人,是萬世師表。

『盍各言爾志』這一課課文教會了我們如何步入中國古文化的文明境界。

Forth 程式語言對比於中國文章,可以說 Forth 程式都算是文言文,大部份都很難直接看懂,需要借助於說明。讀文言文時,我們都借助於古聖先賢對文言文所做過的註解來深入了解文章的意義,所以,程式的說明很重要。除此之外,Forth 程式語言的標準指令還隨著時代演進,主要的精髓是不變的,只在定義上有點差異,研究 Forth 時就不能不注意。下列是一份整理出來的 Fig 指令標準與 ANS 指令標準的對照參考表,你讀 Forth 的古程式時需要用到。


 
Fig Forth vs. ANS Forth
xxx represents any Forth word.
"..." replaces a section of Forth code

State                   Fig                     ANS
                        " This is a string"     S" This is a string"
                        (non-Fig)
compiling:              ' xxx                   ['] xxx >BODY
executing:              ' xxx                   ' xxx >BODY
compiling:              ' xxx CFA               ['] xxx
executing:              ' xxx CFA               ' xxx
                        (KEY)(non-Fig)          KEY
                        (NUMBER)                >NUMBER(not exactly the same)
                        XOR                     INVERT
                        DUP                     ?DUP
compiling:              ."                      ."
executing:              ."                      .(
                        0. 2VARIABLE xxx        2VARIABLE xxx
                        0 VARIABLE xxx          VARIABLE xxx
                        3 PICK                  2 PICK
                        4 ROLL                  3 ROLL
                        ;S                      EXIT
                        <BUILD ... DOES>        CREATE ... DOES>
                        ?1DATA(non-Fig)         KEY?
                        ?TERMINAL               KEY?(any key, not just ctrl-c)
                        ABORT                   START(not ANS)
                        ALLOT                   RAM ALLOT ROM
                        BELL                    $07 EMIT
                        BETWEEN                 1+ WITHIN
                        BLANKS                  BLANK
                        CFA                     BODY>(non ANS)
                        COMPILE                 POSTPONE
                        CREATE xxx SMUDGE       CREATE xxx
                        D*/ (non-Fig)           M*/
                        DLITERAL                2LITERAL
                        DMINUS                  DNEGATE
                        DOWN                    $0A EMIT
                        END                     UNTIL
                        ENDIF                   THEN
                        HERE                    ROM HERE
                        HOME                    $0B EMIT
                        ID.                     .ID
                        IN                      >IN
                        LEFT                    $08 EMIT
                        LFA                     BODY> >NAME N>LINK
                        MINUS                   NEGATE
                        NFA                     BODY> >NAME
                        NOT                     0=
                        PFA                     N>LINK LINK> >BODY
                        R                       R@
                        RIGHT                   $1C EMIT
                        RTE, SMUDGE             RTE,
                        RTS, SMUDGE             RTS,
                        S->D                    S>D
                        U*                      UM*
                        U/                      UM/MOD
                        UP                      $1F EMIT
                        VLIST                   WORDS
                        WORD HERE ALIGNED       WORD
                        [COMPILE]               POSTPONE(usually but see ANS TDS)


2024年11月14日 星期四

我們以數學技術為這個世界扭製出美麗的圖譜

我們以數學技術為這個世界扭製出美麗的圖譜


曾慶潭 Ching-Tang Tseng
ilikeforth@gmail.com
Hamilton, New Zealand
15 November 2024


1. 前言

這本是一篇於 20130416 曾經刊登於 https://how-are-we.webnode.tw/ 網頁的網文。

本人珍惜自己的創作,多年之後的今天,特將此文彙整、集中再刊出於此網頁,以利查找。

文內敘述是標準之繪出平面極座標函數與平面參數方程式函數圖形的方法,程式設計時就根據此方法繪圖。

原文的圖示以淺黃色繪製座標圈,視覺上很不明顯。程式略施修改,改採淺綠色繪製座標圈,效果就明顯的不同。

本文與原文的不同處,是此處刊出了繪圖的源程式。這裡才談 Forth ,程式只適合在這裡刊出。這種大型程式,很難在他處獲得,想詳細了解程式內容則需要下點功夫,程式不容易設計、也不容易用來教會新手使用,繪圖原理不變,才是主要重點。今早貼文前,我試 run 了這個 11 年前設計的程式,照樣能一鍵一圖,效果比當年只能在 XP 作業系統中執行時的效果好太多了。請注意!程式中欲繪之圖的數學函數,只須一列標準的數學表示式便能完成設計,它是 ABC Forth 的特色。程式中一口氣提供了 20 個簡單的範例函數式,想用就這樣用。成品圖包含了徑長的刻度數字,程式展示了圖文並茂的繪製方法。

我所使用的背景系統是這樣的:

Ubuntu 20.04 OS,
加裝視窗模擬軟體 wine 後,
直接執行 ABC660forthV61505.exe,
載入此繪圖程式後,壓下 F2 鍵,
就得到了繪製結果。

Linux 作業系統能夠立即擷圖,只須壓下 shift-print 雙鍵,然後操作滑鼠,框住擷圖區域,就能得到下列圖形:


2. 繪圖原理

描述變數間關係的數學式子稱為函數。


2-1. 平面極座標函數(Polar Function)

極座標函數是一種專門用來描述徑長(r)隨角度(θ)變化的數學式子。人們常因打字不便,而將希臘字母『θ』(讀作theta)改以 t 來取代。

一個最簡單的極座標函數,可以是 r = t 。

進一步,以常用三角函數變化,可以是 r = sin ( t ) ,我們就稱其為:以數學技術扭製出了一個新的數學函數。

更進一步,再以四則運算配合其它各種數學基本函數,更可以扭製出千千萬萬種 r 隨 t 而變化的數學式子。

我們利用扭製出來的數學式子,在以極座標表示的平面上繪圖,就能繪製出各種圖譜,其中,有許多圖譜具有奇特或美麗的花樣。

通常,極座標製圖時,都以右向的水平橫軸,作為角度為 0 的起始基準角度,以其計算出對應的函數徑長值後,在座標圖上繪製出起始之點。隨後,以逆時鐘旋轉方向,依序每次增加一個定量之角度,再根據同一數學函數,依次計算出各個對應的徑長值。將這一序列的極座標數據組,以打點方式標示,或以連線串接起來,直到計算並繪製至完成一圈 360 度時為止,就能獲得相應於極座標函數之完整的圖譜。

圖譜是千變萬化的,所以,您能創作,您能拿產品去小畫家中自行彩繪。

繪製完成的視窗圖中:

(1) 輔以座標及其刻度
(2) 列示單一刻度的對應數值
(3) 列示繪圖所根據的數學函數式子

適合供作學術性方面使用。我們也能提供不帶視窗窗框的無邊顯示圖,更適合在學術論文中使用。另有自定需求者,可以將圖檔輕易載入小畫家系統中,自行補白或美化。

直角座標與極座標線,分別以不同的顏色標示於圖中時,可以增加讀圖時對函數曲線變化的趨勢感,是否顯示這一部份?可以自由選定。

函數曲線比座標線來得重要,因此,我們強調以採用『多隻繪圖筆原理』的繪圖方式製圖。

函數曲線的顯示,可以選擇自己希望的顏色,函數曲線的繪製,也可選用以點狀、細連線、或粗彩線三種不同方式來繪製。


2-2. 平面參數方程式函數(Parametric Equation Function)

參數方程式函數的數學表示方式,與極座標函數不同,製圖方法也有所不同。

參數方程式函數曲線,使用直角座標來繪圖。

參數方程式函數,直接以下列方式表示 x 及 y 所有座標點的變化趨勢:

x = f(t)
y = g(t)

其中 x 與 y 均只單純的隨著另外一個 t 參數而變化,而由 t 所形成的函數,也是無限多種,因此,可繪之圖就能千變萬化。

參數方程式函數中的 t 變數,不適合將其想像成角度從 0 度變化到 360 度的對應關係,比較適合將其想像成是 x 與 y 的座標值是隨時間 t(time) 而變化的趨勢, t 可以是過去(負值)、現在(0值)、與未來(正值)的時間變化。

相對於極座標函數繪圖而言,參數方程式函數繪圖方式的主要不同之處,為函數繪圖時必須另行採用定義域來製圖。

極座標製圖時,函數的定義域為固定的,也就是 t(theta) 從 0 度變化到 360 度,便足以明確顯示函數曲線之變化趨勢。

以參數方程式函數繪圖時,其 t(time) 的定義域,則必須在繪出圖形之前先行選定,選定的範圍沒有限制,可以從很大的負值到很大的正值,選定之 t 的定義域範圍,會決定圖示的品質。

3. Win32Forth 附加 ABC Forth 功能後設計出來之源程式


 
\ Copyright (C) 2013 Ching-Tang Tseng

include winplot.f

\ ---------------------------------------------------------------------

\ 設定整個畫面底色的方法20140325

: SetBkColor ( cilor-obj -- )
  >R
  0 0 screen-width screen-height R> FillArea: demo-dc ;
\ Polar plot

: SetWhiteBG ( -- )
  White SetBkColor
;

\ ---------------------------------------------------------------------

500 to screen-width
560 to screen-height

\ ---------------------------------------------------------------------

\ 令WINPLOT主程式中具有八個單按鍵功能的向量式設計。

DEFER PLOT2
DEFER PLOT3
DEFER PLOT4
DEFER PLOT5
DEFER PLOT6
DEFER PLOT7
DEFER PLOT8
DEFER PLOT9

\ 先讓function key的功能生效,故預先進行下列設定,且不影響後續程式的規劃。

' noop IS PLOT2
' noop IS PLOT3
' noop IS PLOT4
' noop IS PLOT5
' noop IS PLOT6
' noop IS PLOT7
' noop IS PLOT8
' noop IS PLOT9

: Frame
  LTBLUE LINE-COLOR
              0             0 MOVETO
              0 screen-height LINETO
   screen-width screen-height LINETO
   screen-width             0 LINETO
              0             0 LINETO
;

' Frame IS PLOT9

\ -------------------------------------------------------------------

3 INTEGERS I XI YI
6 REALS r Rmax degree t X Y

10000 VALUE N

N ARRAY RR
N ARRAY XX
N ARRAY YY

: r(t) ( -- )
\ {{ r = 1.2345E23 * degree * cos ( 6 * t ) }} ; \ (1)花瓣逐漸縮小圖
\  {{ r = cos ( t ) }} ;                         \ (2)右圓圈圖
\  {{ r = exp ( -1.2 * t ) }} ;                  \ (3)飛機翅膀橫截面圖
\  {{ r = t }} ;                                 \ (4)螺旋線圖
  {{ r = sin ( 4 * t ) * cos ( 8 * t ) }} ;      \ (5)花辦圖
\  {{ r = 2 * ( 1 - cos ( t ) ) }} ;             \ (6)左心形圖
\  {{ r = 2 * ( 1 + cos ( t ) ) }} ;             \ (7)右心形圖
\  {{ r = SQRT ( 10 * cos ( 2 * t ) ) }} ;       \ (8)雙紐線圖
\  {{ r = 10 }} ;                                \ (9)正圓圈
\  {{ r = ABS ( sin ( t ) ) }} ;                 \ (10)上下雙圓圖
\  {{ r = 1 + cos ( t ) * sin ( t ) ^ 2 }} ;     \ (11)三凸輪圖
\  {{ r = 2 - sin ( 7 * t ) - ( cos ( 30 * t ) / 2 ) }} ; \ (12)苔蘚圖
\  {{ r = ABS ( sin ( t ) ) }} ;                 \ (13)up and down circles
\  {{ r = ABS ( 2 * ( 1 + sin ( t ) ) ) }} ;     \ (14)up cardioid
\  {{ r = ABS ( 3 * cos ( 3 * t ) ) }} ;         \ (15)6 leaved Ross
\  {{ r = ABS ( 2 + 3 * cos ( t ) ) }} ;         \ (16)limacon
\  {{ r = 3 * t }} ;                             \ (17)archimedian Spiral
\  {{ r = exp ( 3 * t ) }} ;                     \ (18)logarithmic spiral
\  {{ r = 2 / ( 1 - 0.8 * cos ( t ) ) }} ;       \ (19)橢圓,關鍵在0.8-->0.99
\  {{ r = cos ( 5 * t ) + ( 4.5 ) * cos ( t ) }} ; \ (20)heart to bell
\  {{ r = cos ( 7 * t / 2 ) }} ;

: rEvaluate BASIC
10 REM r evaluate
20 FOR I = 0 TO N
30 LET { degree = I>R ( I ) * 360 / I>R ( N ) }
40 LET { t = ( degree / 180 ) * FPI }
50 RUN r(t)
60 LET { RR ( I ) = r }
70 NEXT I
80 END ;

: RmaxFind BASIC
10 REM Rmax find out
20 LET { Rmax = RR ( 0 ) }
30 FOR I = 0 TO N
40 LET { Rmax = MAX ( Rmax RR ( I ) ) }
50 NEXT I
60 END ;

: PolarToCartesianConvert BASIC
10 REM polar to cartesian convert
20 FOR I = 0 TO N
30 LET { degree = I>R ( I ) * 360 / I>R ( N ) }
    :: { t = ( degree / 180 ) * FPI }
    :: { XX ( I ) = 250 + ( RR ( I ) * cos ( t ) * 200 / Rmax ) }
    :: { YY ( I ) = 250 - ( RR ( I ) * sin ( t ) * 200 / Rmax ) }
40 NEXT I
50 END ;

: DataPreparation ( -- )
  rEvaluate
  RmaxFind
  PolarToCartesianConvert ;

: PolarPlot
BASIC

100 REM circles plot
110 run 1 SetBkMode: demo-dc
120 run black SetTextColor: demo-dc
130 RUN LTgreen    LINECOLOR: demo-DC
140 RUN 250 250  50    CIRCLE: demo-DC
150 RUN 250 250 100    CIRCLE: demo-DC
160 RUN 250 250 150    CIRCLE: demo-DC
170 RUN 250 250 200    CIRCLE: demo-DC

200 REM radiation scale lines plot
210 FOR I = 0 TO 360 STEP 10
220 LET { degree = I>R ( I ) }
     :: { t = ( degree / 180 ) * FPI }
     :: { X = 250 + 200 * cos ( t ) }
     :: { Y = 250 - 200 * sin ( t ) }
230 LET XI = INT ( X ) :: YI = INT ( Y )
240 RUN  XI  YI MOVETO: demo-DC
250 RUN 250 250 LINETO: demo-DC
260 NEXT I

300 REM text print
310 run   4   4 S" Fig. Polar Plot" textout: demo-dc
320 RUN   4  21 S" r=sin(4*t)*cos(8*t)" TEXTOUT: demo-DC

330 RUN   4 460 S" R(max)    = " TEXTOUT: demo-DC
340 RUN   4 480 S" scale/div = " TEXTOUT: demo-DC

350 RUN 490 251 S" X" TEXTOUT: demo-DC
360 RUN 252   1 S" Y" TEXTOUT: demo-DC

400 REM rmax, scale/division floating point values print
410 RUN 98 460 Rmax (FG.) TEXTOUT: demo-DC
420 LET { r = Rmax / 4 }
430 RUN 98 480    r (FG.) TEXTOUT: demo-DC

500 REM x axis and y axis plot
510 RUN ltblue LINECOLOR: demo-DC
520 RUN 250   2  MOVETO: demo-DC
530 RUN 250 504  LINETO: demo-DC
540 RUN   2 250  MOVETO: demo-DC
550 RUN 498 250  LINETO: demo-DC

600 REM axis scale marks plot
610 FOR I = 50 TO 450 STEP 50
620 RUN 250  I MOVETO: demo-DC
630 RUN 260  I LINETO: demo-DC
640 RUN  I 250 MOVETO: demo-DC
650 RUN  I 240 LINETO: demo-DC
660 NEXT I

900 REM function curve plot
910 RUN LTRED LINECOLOR: demo-DC
920 LET XI = INT ( XX ( 0 ) ) :: YI = INT ( YY ( 0 ) )
930 RUN XI YI MOVETO: demo-DC
940 FOR I = 1 TO N
950 LET XI = INT ( XX ( I ) ) :: YI = INT ( YY ( I ) )
960 IF ( ( XI - 250 ) ^ 2 + ( YI - 250 ) ^ 2 ) > 40400 THEN 980
\ 970 RUN XI YI LINETO: demo-DC
\ 970 RUN XI YI BLACK SETPIXEL: demo-DC
970 RUN LTRED BRUSHCOLOR: demo-DC
972 RUN XI YI 2 FILLCIRCLE: demo-DC
980 NEXT I

990 run frame
1000 END
;
' PolarPlot is plot2

\ ---------------------------------------------------------------
\ Top Level program starts here
\ ---------------------------------------------------------------

: WINPLOT       ( -- )
                Start: DEMOW
                StartPos: DEMOW 50 + swap 50 + swap message-origin
                blue line-color
                RANDOM-INIT             \ initialize random number generator
                erase-demo  SetWhiteBG  \ 20140325 ***
                begin   Refresh: DEMOW
                        key             \ handle keyboard interpretation
                        case
                        'O' +k_control  of  open-demo-bitmap            endof
                        '1'             of   1 save-bitmap              endof
                        '2'             of   4 save-bitmap              endof
                        '3'             of   8 save-bitmap              endof
                        '4'             of  16 save-bitmap              endof
                        '5'             of  24 save-bitmap              endof
                        '6'             of  32 save-bitmap              endof
                        'S' +k_control  of  16 save-bitmap              endof
                        'V' +k_control  of  paste-demo-bitmap           endof
                        'P' +k_control  of  print-demo                  endof
                        'Q' +k_control  of  16 print-demo-bmp           endof
\ -----------------------------------------------------------------------------
                        k_F1            of  about-demo                  endof
\ -----------------------------------------------------------------------------
\ 增加可用功能鍵直接繪出圖形功能,讓各種圖形可以疊加存在。

                        K_F2            OF  PLOT2                       ENDOF
                        K_F3            OF  PLOT3                       ENDOF
                        K_F4            OF  PLOT4                       ENDOF
                        K_F5            OF  PLOT5                       ENDOF
                        K_F6            OF  PLOT6                       ENDOF
                        K_F7            OF  PLOT7                       ENDOF
                        K_F8            OF  PLOT8                       ENDOF
                        K_F9            OF  PLOT9                       ENDOF
\ -----------------------------------------------------------------------------
\ 停止但保存所繪之圖,可以回到系統操作畫面,並透過相同操作繼續繪製想疊加之圖。

                        K_F10           OF  EXIT                        ENDOF
\ -----------------------------------------------------------------------------
                        k_cr            of  run-demo                    endof
                        k_cr +k_control of  line-walk                   endof
                        k_esc           of  erase-demo  SetWhiteBG      endof          \ ***
                       k_esc +k_control of                              endof
                'P' +k_control +k_shift of  GetHandle: DEMOW
                                            Setup: ThePrinter           endof
                'C' +k_control          of  false copy-demo-bitmap      endof
                'C' +k_control +k_shift of  true  copy-demo-bitmap      endof
                'X' +k_control          of  false copy-demo-bitmap
                                            k_esc pushkey               endof
                        endcase
                again   ;


\ ---------------------------------------------------------------

: main
  DataPreparation
  winplot
;

page
cr cr
.( Usage: ) cr
.( main --> K_F2 )
main

2024年11月1日 星期五

以運轉原子爐的技術疏通排水管

以運轉原子爐的技術疏通排水管


曾慶潭 Ching-Tang Tseng
ilikeforth@gmail.com
Hamilton, New Zealand
2 November 2024


『排水管堵塞』是生活中經常會遇到的問題,卻很少人知道誰都可以不費力氣、不花金錢、不必求人自行解決問題。

原子爐水系統中的管路堵塞問題,跟我們日常生活中的排水管堵塞問題,本質上問題是同樣的,差別只在於堵塞物有所不同而已。原子爐內的水系統都帶有輻射物質,解決管路堵塞問題時比較麻煩,必須採取輻射防護措施。日常生活中排水管的堵塞物,最糟也只是污穢的不得了、又髒又臭,沒有輻射問題。

我運轉過原子爐十幾年,我們是怎樣解決管路堵塞問題的? 職業關係,我們熟知水系統內各處的水壓,知道水壓的威力。因此,每逢堵塞問題發生時,我們都採取藉力使力的方式解決問題。也就是說,最好的解決辦法就是借助水的壓力把堵塞物沖洗出來。當然,工作前得先知道能接通獲得多高壓力的水?沖洗出來的東西會出現在那裡?沖洗水的蓄積壓力能否建立得起來?運轉工程師熟知這些情況,看著系統藍圖就能指揮工作,而且,十幾年來履試不爽,從未碰到過解決不了的管路堵塞問題。有的時候,管線半堵不堵,水壓蓄積不起來,我們還會故意塞進無害物質,強把半堵塞情況搞到全堵死的程度,然後才再採用蓄積水壓把所有的堵塞物全沖出來。

同樣的,只要你不怕髒,力行孔子的教誨:『多能鄙事』,日常生活中的任何排水管堵塞問題都能自行解決。你只要準備一條能接通到水龍頭的供水軟管,另備一塊抹布,就可以開始工作了。如果情況許可,工作前,請找好被沖出來的堵塞物會出現在那裡的地點,疏通後到那裡去找到成果,能保證完成工作後感到心滿意足。沒有出口點可以觀賞堵塞物時也沒關係,只要堵塞物能被沖到更寬大的排水管道去,也就完成了工作。詳細的解釋,看圖說明會容易些,請看懂下圖:





這個排除管路堵塞物的方法,在我貼文以前,你無法從網路上獲得這種資訊。我的首貼是珍貴的,更貴在讓您學會後自己力行實踐,想轉載者,請尊重出處。

我以專業知識講解工作原理:

一般而言,常人能用口吹氣球所使出的壓力,大約只有 5 ~ 10 psi 。psi 是壓力單位:磅/平方英吋。
搭高樓升降梯、客機著路前的增壓,壓力變化都只低於能吹氣的壓力而已,很多人就已經會產生痛感。沒有表達能力的嬰兒當然會因此而哭泣,同機旅客要諒解。
公稱的大氣壓力為 14.7 psi 。
潛入水中,每下沉十公尺就會增加一個大氣壓的壓力。
家庭常用瓦斯減壓閥的瓦斯輸出壓力為 20 ~ 25 psi 。
自用轎車輪胎壓力為 28 ~ 32 psi 。

那麼,家用自來水的供水壓力,一般為 50 psi 。我要用它 !

這樣的經驗技術數據告訴大家,藉力使力的根據,就在於直接引用自來水的供水壓力疏通排水管。
工作的要領,則在使用雙手捂住溼抹布時,要確保堵塞物之前的管路內能夠蓄積起自來水系統供應的 50 psi 壓力。
另一個要點是:你必須確定堵塞管路這一段內沒有其他歧管,水壓不會因歧管的存在而無法蓄積起來,如果有歧管,歧管的端口也必須完全塞死才行。
堵塞物被疏通的瞬間,手上捂緊的溼抹布會瞬間失去壓力,沖通時壓力的失去會很容易的感覺得出來。

我住紐西蘭,我家廚房排水管的長度超過十公尺,我很體諒我太太每天在廚房裡工作的辛苦,鑑於自己有能力解決任何的管路堵塞問題,我從不限制她讓任何東西流進排水管,反正我就是有辦法、不花錢解決任何管路堵塞問題。移民近 30 年,我親手做過的疏通排水管事件不計其數,包括室外排水瓦管亦然,在南島的但尼丁市做過,從未失手。不怕髒,疏通抽水馬桶堵塞問題也做過,網路上教過 Forth 同好鄧淵源解決他家的問題,沖出了黃小鴨玩具。也指揮過多人合作解決有歧管問題的案例,被指揮的一群人還是好幾個印度人。半堵不堵我強塞廢報紙弄成全堵的情況更是常有,十月份我家就剛做過。採用的都是孔子教誨的: 『吾道一以貫之』的方法,當然,忠恕之道也讓我教會了許多朋友,解決過同樣的問題。

排水管壁通常都是光滑直通到底的,否則管路更容易堵塞,採用這套辦法解決問題時,不要擔心沖不出堵塞物。這種沖法,是靜壓力蓄積起來後才開始會有動作,工作時管路的封閉狀態,就像油壓千斤頂內蓄積的液態油壓,壓力大到無處可去時,就會把唯一可動之排水管光滑管壁上的堵塞物沖走了。水的最大壓力頂多就是 50 psi ,不會沖爆排水管,我操作過 250 ~ 300 psi 的水系統管路疏通工作。

操作記錄:
200241125 廚房排水

2024年10月15日 星期二

論語是中國文化的基本教材

論語是中國文化的基本教材


曾慶潭 Ching-Tang Tseng
ilikeforth@gmail.com
Hamilton, New Zealand
16 October 2024


2500 年前,中國出現了一位所有老師的老師,名叫孔夫子,他教過三千多個優秀的學生,他不寫書、只教課,學生就把他的言論一條一條的記錄下來,累積成了 20 個章節的書本,書名叫作"論語",書裡面記錄的內容,全是教大家做人處事的方法。

台灣以前的高中學生都要專門學習"論語"這本書,學習"論語"的重點在要求學生進入社會之後,要身體力行書中所講的道理,並力行使用從"論語"學到的精神來做人處事。

兩千多年來,很多中國人都說,只要能做到半本"論語"書中的要求,全世界就和平了。


舉兩個簡單的例子來簡單介紹論語這本書:


論語書中記錄的第一章第一條,孔子就告訴大家如何求學。孔子說:求學的時候能夠經常複習、實際練習是一件非常愉快的事情。

本網頁的貼文,偏重於與 Forth 相關的技術與文藝,顯示我終生學習 Forth ,每篇貼文都代表了我恆在複習、練習著 Forth ,我當然很快樂。


在第十六章第四條中,孔子談到了交朋友方面的事情。

孔子說:

好朋友有三種,個性正直的人,做事講信用的人,知識淵博的人,都是好朋友。
翻譯成簡單的英文,三種好朋友就是:
1. He is always insisted to do the correct things.
2. He is trustable.
3. He knows a lot.

壞朋友有三種:不顧尊嚴而做壞事的人,不誠實而愛討好他人的人,喜歡狡辯而顛倒是非的人,都是壞朋友。
翻譯成簡單的英文,三種壞朋友就是:
1. He is always to do bad things without integrity.
2. His speaking is always untrustable.
3. He likes to argue and always upside down the right things.

每個人的一生都有許多朋友,我就是根據孔子的言論來區分朋友的好壞。

論語原為古時候以最精簡扼要的文言文記錄,後人讀論語都讀經過名人註釋與講解的讀本,現代人則讀經過名家收集了註釋而且已經寫成白話文或稱語體文的書本。


隨文附貼一份溫馨感人的漫畫:『第一天學 FORTH (Starting FORTH)』,我永遠都像是第一天在學 Forth ,而且也希望天天能這樣學 Forth 。

漫畫的原創人並不是我,我只是剪貼原漫畫,再加油添醋般加進了『4th』。原創漫畫的書名及作者的資料如下:
『Strong and aggressively drawing stories』by Herluf Bidstrap, made in Praha, Publishing house ROH, Bratislava, 1961
感謝 Herluf Bidstrap 這份深具學習 Forth 意義的創作。

第一天學 FORTH (Starting FORTH) 連環圖畫:

.

2024年10月6日 星期日

一百個例題 (96 ~ 100)

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

2024年10月5日 星期六

一百個例題 (91 ~95)

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

2024年10月4日 星期五

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