2024年9月27日 星期五

一百個例題 (71 ~ 75)


Ching-Tang Tseng
Hamilton, New Zealand
28 September 2024

第(71)個範例,收集了四套以高階定義方式設計出來的雙整數乘法指令 D* 。

百例中的後半段範例,大部份都是與我在發展系統的過程中相關事務的研究,如果想學系統設計,這些事務都得搞清楚。能有四套不同的作法實現程式,這個範例就不可能是一次或一天內落實下來的。有的時候,為了怕忘記程式執行的原理,我還得在收集之後,又回頭去查貢獻者當時所提供的說明文獻。

您若仔細研究這四套的各個設計,可以總結出 Forth 設計程式時的基本精神,就是在搞堆疊上數據的翻來覆去以理順秩序。這就是 Forth 跟所有傳統程式語言完全不一樣的地方。範例介紹進度已經超過 70% 了,更不能回頭重談 Forth 系統中最基本的堆疊操作指令。還有,為了暫時安放數據,這些程式是不是也顯現了開始投機使用可以臨時暫放數據的地方?不放堆疊,以減輕堆疊操作的負荷?這些精神,我在前面的範例中都說明過,所有設計系統發展到後來的人,都會走上同樣的研究路途。

我得提醒大家,許多 Forth 系統,在提供雙整數的四則運算指令時,只有 D+ ,就也有 D- ,因為減法只是將減數取負值後再與被減數 D+ ,就得到 D- 了。但是,大部份 Forth 系統不提供 D* 指令,更不提供 D/ 的指令,為什麼?我長期閱讀國際論壇上的討論,真的發現,只有一個理由可以解釋這種現象。就是 D* 還很容易設計, D/ 就是比較不容易設計,所以才都不提供 D/ 指令。

討論範例時主要使用的 Win32Forth 系統,沒有 D* 也沒有 D/ 指令。

能設計出 D/ 指令,就成了一個寶貝資源,我是會設計,我也將成果留進了庫存程式。方法只是一再的執行移位與相減, 128 位元作 128 次就得到結果了,算法沒有特出處,所以程式也只留供自己庫存。至於平常例行的使用,很少用到 D/ ,若非用不可,我設計的系統已經具有涵蓋範圍更為廣泛的大整數除法,使用起來根本沒有位數限制,所以,我自己設計的 D/ 指令就沒有必要再添入系統了。

Forth 系統中有一個為了提高計算精確度而設計的 */ 指令,裡面涉及的除法,只是用雙整數除以單整數得到單整數,不是雙整數除以雙整數得到雙整數。

四則運算是一切計算的基礎,我也只憑此基礎設計出所有的函數,最後形成了所有的數學體系。雙整數的運用,在低位元的時代,也就是 32 位元以下的時代,因為最多只有 10 位數精確度,可能還具有常用價值。但在進入 64 位元時代之後,一個單整數就具有將近 20 位數的精確度,再搞雙整數去得到 40 位數的精確度,就有點不太實際了。
:

 
\ Dmultiply.f

\ In FSL_UTIL the following is defined:
: d*1  ( d1 d2 -- dprod )   \ double multiply
  dup 3 pick xor >r
  dabs 2swap dabs
  udm*
  r> 0= IF exit ENDIF
  dxor 2swap dxor
  1. UD+c >r                \ ud1 ud2 carry1 -- ud1+ud2 carry2
  2swap r> UD+c drop ;

: d*2   ( d1 d2 -- dprod )
  >r >r
  over r@ um*
  rot r> * +
  rot r> * + ;

: d*3 ( d1 d2 -- d1*d2 )
  >r swap >r 2dup um* 2swap r> * swap r> * + + ;

: d*4 	( multiplicand . multiplier . -- product . )
  3 PICK * >R  TUCK * >R  UM*  R> +  R> + ;

第(72)個範例,是我自創之系統中一個典型的檔案操作源程式。

這個主題,是一個可以完全獨立來討論的範例,您無法在全世界網上的任何角落找到這樣的資源。對能搞檔案操作程式的設計者而言,一方面是知者不說,另方面是做法不同,所以網上沒有相關訊息。我公開展示範例,等於是將自己的技術傾囊相授,能否完全吸收?就只能看您自己了。

一般 Forth 系統,能給您的檔案操作指令,就是 create-file, open-file, read-file, write-file, close-file, copy-file, delete-file ..... 等等等之類的現成指令。可是,您若不照我這樣的設計來使用,您可以自己試試看,可能連一個檔案也開啟不了,還可能把某個檔案搞爛了,最後想殺也殺不掉。

我在發展這個範例程式期間,毀掉過 CentOS7 好幾次,最後才放棄了,此後便堅決改用 Ubuntu 。我家人有從事於軟體設計工作者,不接受我批評 CentOS7 不好的意見,說我這樣玩系統,當然會把系統搞壞。他也不肯幫我另裝 Ubuntu ,他認為我老是在搞病毒與駭客的玩意兒,很危險。只肯跟我寫 email 溝通,最好不要安裝我創作的任何程式,以免令他的公司受害,他就會丟掉飯碗。我也不置可否,反正不能跟他談 Forth 。

這就是為什麼程式設計人,搞到後來,免不了就得自行花很多時間,細讀許多技術文獻,搞清各種所謂的『協定』(protocol)之主因。我也是這麼過來的,否則,您不可能憑猜測,就知道可讀可寫檔案,屬性的位元花樣狀態為何?每個作業系統不同版本的規格也都不一樣,我連書也沒買過,學也沒上過,人也沒問過,就在網上長期搜索,找出、讀通、試成、使用這些資源。歸根結底,就是我自己知道該怎麼做,才做得成檔案操作之事。我有個生產鋁合金錠的大亨朋友,當年他就是這麼說我的,公司的錢不怕讓我隨便用,因為我勇於創新。

大家都已經不能去 CentOS7 跑這個範例程式了,可是,將其重新載入 Lina64 ,照樣能用。換句話說,我已經把它寫成通用格式了,是長久發展後的結果。這個範例在 Win32Forth 系統中不能使用,也不是我在 Ubuntu 環境中的落定程式,這是在 CentOS7 中發展系統時產生的可貴記錄,只適合當虛擬教材使用。

被操作檔案可以放置的地方,有很多可以安排的方式,可以放在 Forth 系統之內,也可以放在系統之外。能放在系統的浮動記憶體內,也可以放在系統的固定記憶體內。範例程式內,都留有這些安排方式的程式技巧。這又牽扯到必須對系統內部結構必須了解的知識要求,要講整個系統的記憶圖譜,有點麻煩,此處不深談。

請注意!我在歸納學問之後,發現,只要設計出這三個指令: GetFile , NewFile , SaveFile 。就夠用來寫任何具有檔案操作需求的程式了。您想執行 delet , copy , move ..... 等等額外要求時,別在程式內搞,請自己用滑鼠去點開資料夾,單用點選操作來達到目的,才是合理的做法,因為『眼見為真』,看著做,就不會出差錯。想要設計出所有的功能指令,我也辦得到,讀讀書,寫程式,就能完成,但沒必要。

搞檔案,是個非常龐大的話題,您難得見到,我只講使用這三個指令就夠了的範例。現在,您只需了解這些東西,按照我在範例程式中旁白所給出的使用方法,搞自己的東西,就能全用程式來解決問題。改天,我也會走向上網搞東西的發展,專搞遠方的對象。請特別注意!它們也必還是只能以檔案作為一個最基本的單元,才能成為被遠方來客處理的東西。

此範例雖在 Win32Forth 系統中不能跑,您若有心,使用 see 指令去看我在該系統中的設計,其內容,仍然是這個範例的設計,我的檔案操作設計,基本精神只有一個,永遠也不想再改。至於不同作業系統會產生不同處理現象的問題,不是重要問題。在微軟的作業系統環境中,可以產生虛擬檔案,用完就不見了。在 Linux 公益作業系統環境中,只要產生,就絕對不是虛擬檔案。這種問題,並不重要,重要的是,要知道檔案到底存放在那理?

我用這些設計做過許多事情,例如:

我的系統能使用語音輸出技術,自動讀出文字檔案的內容。
我的系統能動用好幾個檔案,叫系統合奏音樂。
我的系統能透過網路,自行發動送出百萬封信件的攻擊訊息。
我的系統能跟任何共存於系統中的所有其他軟體進行檔案溝通,並行運作。
我的系統能跑別人的可執行檔,也能再跑許多的另一個自己。
我的系統能開啟別人禁止開啟的檔案,細看究竟。

有說不完的功能,端賴於我自己想不想用而已。

百例中有許多範例的主題,不討論數學計算,換句話說,範例程式做了許多事情,不只是數學計算,這一個範例就是與數學計算無關的程式設計。
:

 
\ FILEOP.F

\ This file operation program works in Lina64 under CentOS7
\ Author: Ching Tang Tseng
\ Date  : 20160310, Hamilton, New Zealand

WANT ALLOCATE

\ -rw0rw0rw0 = 110110110b = 438d : all are R/W enable
438 CONSTANT R/W 
VARIABLE FileID
VARIABLE Fptr
VARIABLE Frem
VARIABLE Flen
VARIABLE Fsize   10240 Fsize !

\ (1)floating Fadr
\ : Fadr PAD 4096 +  ;

\ (2)allocate Fadr ???
\ Fsize ALLOCATE DROP CONSTANT Fadr

\ (3)fixed Fadr: 100 KB below EM
\ EM HERE - . --> get 33425420 --> 33 MB free spaces
: Fadr EM 102400 - ;

: SetUpFptrFrem   ( -- )
  Fadr Fptr ! 
  Flen @ Frem ! ;
: (FILE.)   ( -- addr len )
   Fadr Flen @ ;
: FileType  ( -- )
  CR (FILE.) TYPE ;
: FileDump ( -- )
  CR (FILE.) DUMP ;

\ Beware! only a R/W attributed file can be manipulated 
\ S" Filename.f" GetFile = "Filename.f" GetFile  
: GetFile ( addr len -- )
 Fadr Fsize @ 0 FILL
 R/W OPEN-FILE
 IF CR ABORT" OPEN-FILE error?" THEN FileID !
 CR ." File ID is : " FileID @ .
 Fadr Fsize @ FileID @ READ-FILE
 IF CR ABORT" READ-FILE error?" THEN
 DUP Flen ! 
 CR . ." Bytes has been read!"
 FileID @ CLOSE-FILE
 IF CR ABORT" CLOSE-FILE error!" THEN
 SetUpFptrFrem ;

\ use PAD area create all blanks
: NewFile ( addr len n -- )
  Flen !    PAD Flen @ 32 FILL
  R/W CREATE-FILE
  IF CR ABORT" CREATE-FILE error!" THEN FileID !
  PAD Flen @ FileID @ WRITE-FILE
  IF CR ABORT" WRITE-FILE error!" THEN
  CR Flen @ . ." Bytes has been written!"
  FileID @ CLOSE-FILE
  IF CR ABORT" CLOSE-FILE error!" THEN
  SetUpFptrFrem ;

\ Beware! Flen must be set, before you WriteFile
: SaveFile ( addr len -- )
  R/W CREATE-FILE
  IF CR ABORT" CREATE-FILE error!" THEN FileID !
  CR ." FileID is: " FileID @ .
  Fadr Flen @ FileID @ WRITE-FILE
  IF CR ABORT" WRITE-FILE error!" THEN
  CR Flen @ . ." Bytes has been written!"
  FileID @ CLOSE-FILE
  IF CR ABORT" CLOSE-FILE error!" THEN ;

\ for simple testing usage:
\ S" This is a simple test." SendText>F
\ : SendText>F ( adr n -- )
\   DUP Flen ! Fadr SWAP MOVE ;

\ Frem, Fptr are to be used for other testing.

第(73)個範例,仍是用來考驗 Forth 印出資料能力的另種典範。

在很前面的範例中,早就介紹過印出數字或資料的範例了,為什麼這裡又留了一個這樣的範例?當然是有理由才留的。

在電腦系統中,數字是數字,字元是字元,兩者的意義是完全不同的。只印出數字時,程式很容易寫,只印出文字時,程式也很容易寫。想要連起來印,中間不得有空格時,程式就不好寫了。所以才留這個簡化到最精簡程度的範例程式。

問題的來源是一個希臘生物研究中心的網頁,在該國經濟崩盤後,網頁也消失了,所以貼示網址已無意義,只剩這個例題還能參考。

生物學方面的研究,最需要藉著電腦作為工具進行研究的地方,可能就是處理大量的基因圖譜資訊。人體有十萬個以上的基因序列,而基因的最基本元素,則被歸納出只有四種花樣,簡稱 AGTC ,然後,所有的基因序列,就只根據這四個基本單元重複排列而成。大家對這種東西都不熟,只聽說過基因『比對』的術語,可能比較熟悉。這樣講,大家就比較熟了。

附註:生物術語
A :Adenine 腺嘌呤
G :Guanine 鳥糞嘌呤
T :Thymine 胸腺嘧啶
C :Cytosine 胞嘧啶

在大海中撈針,也像在網路上『比對』出指定資訊,也像在十幾萬個基因序列中『比對』出指定花樣。電腦可用的地方,除了計算很快,就是『比對』速度也很快,所以,近代生物科技的研究,電腦可以幫得上很大的忙。這個範例中沒有去搞『比對』之事,只是專搞『字連數目』的印出方法。電腦的比對運作,使用者是看不見的,您除了印出來看,以了解電腦程式做得對不對,別無它法。所以,想發展這方面的程式時,首先該搞的事情,就是把『字連數』的資料印出來看。我在自行發展浮點數計算系統時,也是從這樣的起點開始做起,先準備好印出浮點數的指令,才寫系統程式,測試程式設計的對不對時,只能把浮點數直接印出來看,才能知道自己設計的程式對不對。

範例分別用純粹的 Forth 與用我自創的 BASIC 式寫法設計程式,比對輸出結果後,功能完全一致。輸出時的關鍵,在印字時,使用 emit ,印數目時,使用 dot(『.』)。

在傳統 BASIC 式語言中,沒有 emit 可用,每印一筆資料後,都會跳列,就無法完成這類工作要求。

我自創的 BASIC 中設計了一個天下無敵的 RUN 指令,這在最原始的設計人 Charles H. Moore 提供的原始程式中,也沒有這樣的設計。他的 RUN 指令,只能跑上列單一個剛剛設計好的 BASIC 式程式。這一特點,是令我能夠青出於藍的最主要部份。我相信發明人與另一位改寫成 F83 規格的程式設計者,早就看過我網上發表過的訊息,這兩位前輩深明其中的道理,但我不曾跟他們聯絡過,就讓彼此的技術都埋在心裡。

關於基因『比對』的技術,此後可能會成為一個全世界想要研究的主題。不久前,我在 Euler’s 問題網頁見到一個題目,就是從一堆成千上萬個亂七八糟毫無規則的數陣中,用程式,比對出一個例如是十位數的數字花樣來。這樣的問題,就是『比對』應用的範例,我自己寫成過程式。也自己再重新使用 Fibonancci 級數產生夠亂的大數花樣,繼續考驗程式,設計得非常成功。驗證結果時,也是把比對完成後的位址留下,再根據此址,把數字當文字印出,人工觀察程式執行後看得到結果對不對?

問題若擴充到網頁上的應用,也是這樣。我從駭客技術中,讀到破解進入他人有對談網頁可以操作的方法。就是根據網址,下載該網頁的網頁用語 HTML 程式源碼,把它存進一個檔案,然後用程式去比對出要求輸入密碼時的固定指令在那裡?程式快速找出位址,予以定址後,就能單改該處的程式結構成為不問密碼,立刻回存修改後的源碼,存回同一網址,此時,就能讓想要入侵者,沒輸入密碼,就能進入網頁,遂行操作。快速比對修改源碼結構的動作,可以寫成程式,完全自動執行,侵入速度就能很快,而且程式可以到處快速測試,能侵入時就侵入,侵入完後就恢復沒入侵前的源碼狀態,不留痕跡。有防火牆不能侵入時,就快速離去。
:

 
\ PRINTstudy.f

: main1 ( -- )
  cr
  73 65
  do
     13 1
     do
          j emit i .
     loop
     cr
  loop ;

2 integers i j

: main2 ( -- )
  basic
10 run cr
20 for j = 65 to 72
30 for i = 1 to 12
40 run j emit i .
50 next i
60 run cr
70 next j
80 end ;

\s

main1 
A1 A2 A3 A4 A5 A6 A7 A8 A9 A10 A11 A12 
B1 B2 B3 B4 B5 B6 B7 B8 B9 B10 B11 B12 
C1 C2 C3 C4 C5 C6 C7 C8 C9 C10 C11 C12 
D1 D2 D3 D4 D5 D6 D7 D8 D9 D10 D11 D12 
E1 E2 E3 E4 E5 E6 E7 E8 E9 E10 E11 E12 
F1 F2 F3 F4 F5 F6 F7 F8 F9 F10 F11 F12 
G1 G2 G3 G4 G5 G6 G7 G8 G9 G10 G11 G12 
H1 H2 H3 H4 H5 H6 H7 H8 H9 H10 H11 H12 
 ok
main2 
A1 A2 A3 A4 A5 A6 A7 A8 A9 A10 A11 A12 
B1 B2 B3 B4 B5 B6 B7 B8 B9 B10 B11 B12 
C1 C2 C3 C4 C5 C6 C7 C8 C9 C10 C11 C12 
D1 D2 D3 D4 D5 D6 D7 D8 D9 D10 D11 D12 
E1 E2 E3 E4 E5 E6 E7 E8 E9 E10 E11 E12 
F1 F2 F3 F4 F5 F6 F7 F8 F9 F10 F11 F12 
G1 G2 G3 G4 G5 G6 G7 G8 G9 G10 G11 G12 
H1 H2 H3 H4 H5 H6 H7 H8 H9 H10 H11 H12 
 ok

第(74)個範例是專講數值積分的最基礎範例,雖基礎,卻已足供實際應用。

我大學時教我數值分析的王九龍老師曾經明確的告訴我們:『數值積分可以積出很精確的結果,數值微分的結果則比較不準。』

那時,我們沒有太多的機會跑大型電腦,對老師授業的重點沒有太強烈的感受,五十幾年後的今天,不一樣了,我自己能創作系統,更對這樣的教誨充滿了無限的感激,因為,我能仔細體會得出教師名言的實際意義。我先講實例說明,再討論程式。(貼出此文的日期是教師節 9 月 28 日)

量測中子時的積分技術。

新的原子爐,自己的偵測系統無法量得讀數非常微小的中子源微少的中子通量訊號,少到什麼程度?大約是每秒只能有個位數的中子量。但新爐子就只能用這個量來啟動,邊測量邊添加燃料入爐,一直要測量到爐子能夠自持臨界反應的零功率後,我再測出可以添加多出大約 6% 的燃料量,以便控制系統能涵蓋掉這多出來的 6% 反應能力,才終止實驗,停下爐子,拆除我的實驗設備。

由於爐子在臨界之後,便可以擁有很強的輻射,能持續的激發已放在爐中的中子源,此後的重新起爐,就不必再藉助於另外加裝非常靈敏的中子量測裝置了。我記得很清楚,當時,我啟動測量裝置時,每秒的中子讀數是 2.67 個,系統自己安裝之永久性三套中子偵測裝置,必須要等到我的測量讀數為每秒超過百萬個中子以後,指針才有指示。臨界過一次之後,三套儀器就都能自行工作了。

每秒鐘 2.67 個中子,是怎麼測量出來的?只憑內含武器級鈾(鈾濃縮度大於90%)的中子偵測器,每次累積測量半小時,等到顯示裝置長出一個 Gauss 分佈狀的小饅頭專有信號後,我把曲線下方的面積積分出來,總數除以 30x60=1800 秒鐘,才得到這樣的數目。實驗過程中,我要不斷地搞信號積分,把每秒鐘的中子通量讀數算出來。

這種所謂的個位數讀數,必須調整我的偵測器與中子源的距離來完成。太近了,讀數太大,雖比較靈敏,等到臨界時,我的偵測器讀數就會超飽和,滿標而失效了,就不能決定燃料的可多添加量。太遠了,讀數太小,因為不靈敏,低讀數的測量就得進行更多次,每次加添過一次燃料後,我的測量至少就都得等上半小時,這樣搞爐子,從早上開始搞起,等到天黑了,也搞不出來。所以,我得自己調整出起始讀數就是個位數,這樣的實驗,我做了兩星期。

這就告訴了大家,積分確實可以積得很準,準到能夠搞原子爐啟動,與搞原子彈爆炸。

微分呢?不太準,不準到什麼程度?新聞上有個例子可以解釋。

波音 737 MAX 老摔飛機,最後找出了問題所在,就是機身上的攻角測量儀器讀數有問題,才亂搶控制權,讓飛機栽下來的。攻角怎麼量?量的就是機體跟水平方向的夾角,靠空速管內流體的變量,比對變率後,提供數據,這儀器肯定靠微分(變率)產生信號,結果,微分天生就做不準,電腦再好,就是算不準,就出問題了。結論是,微分的不準,可能不準到會摔飛機。

在原子爐的控制實務中也有個微分很難取得準確值的例子。

控制原子爐時,強調要把中子的倍增週期控制在 10 秒鐘左右,而週期的換算公式是 1/T = dln(N)/dt ,其中: T 是週期,N 是當時量得的中子通量,中子通量之對數標值的變率就是週期的倒數。中子通量是數位式的計量值,如果採用把每秒幾個中子的讀數曲線進行這樣的微分,肯定是來不及反應出準確的週期值來,這樣操控原子爐就很危險。實際的控制系統採用二極體的電壓對電流的對數式特性曲線來對應到中子通量,其變率是類比電腦的反應速度,便可以即時指示出中子倍增週期為幾秒。這種控制信號不能採用對中子通量之數位訊號進行數值分析的技術來獲得微分變率,它會不準也不夠快。

回頭講我們的範例程式,積分很簡單,就是用微小的寬度 dx 乘上該點的函數值 f(x) ,然後從頭到尾,把所有的微小面積加總起來,就得到結果了。很不幸,函數曲線斜率變化較大時,亂算面積,所得結果,誤差就會很大,於是,我就把光算長方形,取低位左邊長、中央長、高位右邊長的結果全算一遍來比較,另外也搞個梯形面積算法也參加進來,全都不準。最後,用到了所謂的辛普森算法,才終於有了比較好的結果,於是,範例就留在這裡了。這些算法在 Rosetta code 網頁上有詳細說明,也有很多程式可以比較,有興趣者,請自行前往。範例(74)是純用 Forth 寫成的程式,附帶方法說明。範例(74-1)則是全用 BASIC 語法寫成的程式,附帶執行結果。

積分準確的程度,甚至於可以用來作為系統的函數值產生源,問題在要細分成萬點以上的面積計算才有效果時,算得太慢了,所以才沒有人這樣使用。例如: 1/x 的積分就是對數函數 ln(x) 之值,還有好幾個函數都能這樣設計。

積分原理,也不僅只有這幾個而已,我還知道例如:辛普森 1/3 法則、 3/8 法則、韋爾德(Weddle)法則...等等,都可以用,數值分析的課本中有更多的方法可以參考,我只能點到為止的提及它們的名稱。

本網頁 20180401 的文章,是全用英文貼出介紹數值積分(Numerical Integration)的文章,其內容就是取材自這個範例。
:

 
\ (74)SimpsonIntegration.f

fvariable step
 
defer method ( fn F: x -- fn[x] )
 
: left                    execute ;
: right  step f@       f+ execute ;
: mid    step f@ 2e f/ f+ execute ;
: trap
  dup fdup  left
      fswap right f+  2e f/ ;
: simpson
  dup fdup  left
  dup fover mid 4e f* f+
      fswap right f+  6e f/ ;
 
: set-step ( n F: a b -- n F: a )
  fover f- dup 0 d>f f/ step f! ;
 
: integrate ( xt n F: a b -- F: sigma )
  set-step
  0e
  0 do
    dup fover method f+
    fswap step f@ f+ fswap
  loop
  drop fnip
  step f@ f* ;
 \ testing similar to the D example
: test
  ' is method  ' 4 -1e 2e integrate f. ;
 
: fn1  fsincos f+ ;
: fn2  fdup f* 4e f* 1e f+ 2e fswap f/ ;
 
7 set-precision
test left    fn2  \ 2.456897
test right   fn2  \ 2.245132
test mid     fn2  \ 2.496091
test trap    fn2  \ 2.351014
test simpson fn2  \ 2.447732

\S
Numerical integration
Write functions to calculate the definite integral of a function f(x)     
using all five of the following methods:

  rectangular
  left
  right
  midpoint
  trapezium
  Simpson's

Your functions should take in the upper and lower bounds (a and b),   
and the number of approximations to make in that range (n).

Assume that your example already has a function that gives values for f(x).

Simpson's method is defined by the following pseudo-code:

h := (b - a) / n
sum1 := f(a + h/2)
sum2 := 0

loop on i from 1 to (n - 1)
sum1 := sum1 + f(a + h * i + h/2)
sum2 := sum2 + f(a + h * i)

answer := (h / 6) * (f(a) + f(b) + 4*sum1 + 2*sum2)
Demonstrate your function by showing the results for:

f(x) = x3,  where x is [0,1], with 100 approximations. The exact result is 1/4, or 0.25.
f(x) = 1/x, where x is [1,100], with 1,000 approximations. The exact result is the natural log of 100, or about 4.605170
f(x) = x,   where x is [0,5000], with 5,000,000 approximations. The exact result is 12,500,000.
f(x) = x,   where x is [0,6000], with 6,000,000 approximations. The exact result is 18,000,000.

\ (74-1)RectangleIntegration.f
\ Author: Ching-Tang Tseng
\ 20161107, Hamilton, Hamilton, NZ

9 reals a b h sum1 sum2 x f(x) f(a) f(b)
2 integers n i

: function ( -- )
basic
10 let { f(x) = 1 / x }
20 end ;

: Setting ( -- )
basic
10 let n = 1000
20 let { a = 1 }
30 let { b = 100 }
40 end  ;
\ All parts above are adjustable 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 }
40 run function
50 let { sum1 = f(x) }

60 let { sum2 = 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 }
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 ) * ( f(a) + f(b) + 4 * sum1 + 2 * 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 }

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 }

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 }

50 for i = 1 to n
60 let { x = x + h / 2 }
70 run function
80 let { sum1 = sum1 + h * f(x) }
90 let { x = x + h / 2 }
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 }

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  }

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 ;

18 sigdigits !
: main ( -- )
cr cr ." Integration of 1/x from 1 to 100 = Ln(100) " cr
L-Rectangle
R-Rectangle
M-Rectangle
Trapezoidal
Simpson
cr cr ." Ln(100) = " 100.0e0 fln f. cr ;

main

\S

Integration of 1/x from 1 to 100 = Ln(100)

  L-Rectangle integration =        4.65499105751467615  
  R-Rectangle integration =        4.55698105751467615  
  M-Rectangle integration =        4.60476254867837518  
  Trapezoidal integration =        4.60598605751467615  
      Simpson integration =        4.60517038495714217

Ln(100) = 4.60517018598809137
 ok

第(75)個範例是以高斯(Gauss)法求解多元一次聯立方程式的實用程式。

我在 1980 年代就已經很熟這種專門用來求解下列矩陣方程式的技術:

[A] * [X] = [Y], 已知 [A] 及 [Y] 求 [X] 。

方法也很簡單,我都能背下來了,根據線性方程式的定義,基本運算原理就是:『任何一列乘以或除以任意常數,加或減到另一列去,整個矩陣方程式的本質不變』。

於是,像前面的範例展示過的矩陣分解法ㄧ樣,把 [A] 搞成上三角矩陣後,就能從右下角最後一個元素的直接對等關係,得到該 X(n) 的解答,然後,反推上去,次一個 X(n-1) 也有了答案,依次類推,最後,所有的 X(i) 的答案就都得到了。其中, i = 1, ... n 。

話雖如此,想寫成程式,並不簡單。一方面是運算時,後面跟著的 Y(i) 值也得跟著變。另方面是為了維持浮點數運算的最佳精確度,矩陣在操作前,最好先挑出數值為最大的元素來,將其轉置到 [A] 矩陣的左上角當樞紐元素,然後才開始演算,這樣,所得的結果才會最精確。還有,免不了有時會碰上無意義的(trivial)線性相依情況,設計程式時,也得一併解決,成果才能萬用。

講完了原理,改為強調一下這個程式的特色,才是今天介紹範例的重點。

請注意程式的寫法,不是 BASIC 語法,每一列前,不用標號,沒有 BASIC 指令。也不是純 FORTH 語法,你看不到堆疊操作指令,例如: dup, over, swap, drop ... 等等。

程式中所有的算式,都可以寫成純粹中算符(infix)的格式,只須記住,在 『{{ ..... }} 』內,只對浮點數運算,在『 [[ ..... ]] 』內,只對單整數運算,記住這樣的規矩就夠了。如果另外還想進行複數的中算符運算,那就多記一個,在 『{[ ..... ]} 』內,只對複數運算。

這樣的系統性能,在我提供的 ABC Forth 系統中,永遠都有,大家可以直接載入此程式來看結果,就能明白。

事實上,我早在 1980 年代,就已經發展出這種單用中算符寫程式的技術了,那時還沒有搞出 BASIC 程式語法中所需要的指令。最大的差異,就在邏輯語法, Forth 與 BASIC 兩者互相相反。固定迴路寫法,也完全不同, Forth 用 DO ..... LOOP , BASIC 用 FOR ..... NEXT ,上下限用法完全顛倒。

這個程式改寫出來的日期,標示為 2009 年,更原始的純 Forth 格式程式,也是我親手設計的,在『FORTH期刊』第 9 期,民國七十六年(1987)年一月份出版, p.58 中就有我寫的程式,該文作者項舉勁在文末參考文獻中,提到使用我的這個設計,是純用傳統 Forth 程式語言寫成的,程式早就在許多地方試用過,都沒有問題。

我從使用純粹的 Forth 開始,於 1980 年代,發展出中算符技術,就像這個範例中所留下來的記錄一樣。 1990 年代以後,我改走 BASIC 語法技術,演變出了今天的成果,就是現在大家所見與所用的 ABC Forth 。

我自創的系統中,還有許多其他的隱藏特性,沒有用到之前,我就沒講,看到了再說。我自知,就算我介紹完ㄧ百個例題,也不可能討論得完所有我自創系統的特點。更何況,我一直還在發展,有些新東西,還都是配合現行作業系統才發展出來的新特性。這無限多種的特殊技術,是講也講不完的東西,我們只能見一個討論一個,這也是 Forth 程式語言的本質,一直都能繼續發展。
:

 
\ (75)高斯消去法解多元聯立一次方程式2009-08-22

INTEGER #ROW INTEGER ITEMP             \ #ROW為方陣量,本例為3X3矩陣
INTEGER M INTEGER M1 INTEGER M2
REAL FTEMP
10 ARRAY Y    10 ARRAY X               \ 矩陣上限宣告,暫設為10X10
10 10 MATRIX AA

: 輸入數據
  [[ #ROW = 3 ]]
  {{ AA ( 1 1 ) = 3 }} {{ AA ( 1 2 ) = 2 }} {{ AA ( 1 3 ) = 1 }}
  {{ AA ( 2 1 ) = 2 }} {{ AA ( 2 2 ) = 4 }} {{ AA ( 2 3 ) = 5 }}
  {{ AA ( 3 1 ) = 4 }} {{ AA ( 3 2 ) = 1 }} {{ AA ( 3 3 ) = 1 }}
  {{  Y ( 1 ) = 17  }} {{  Y ( 2 )  = 41 }} {{  Y ( 3 ) = 16  }}  ;

: 高斯處理
  #ROW 1
  DO
    [[ M = I + 1 ]]
    #ROW 1+ M
    DO
      #ROW 1+ M
      DO
  {{ AA ( J I ) = AA ( J I ) - ( AA ( K I ) * AA ( J K ) ) / AA ( K K ) }}
      LOOP
  {{ Y ( I ) = Y ( I ) - AA ( I J ) * Y ( J ) / AA ( J J ) }}
    LOOP
        M #ROW =
        IF  {{ FTEMP = AA ( M M ) }}
            FTEMP F0=
            IF ." Equation disparity! " QUIT THEN
        THEN
   LOOP ;

: 得到解答
  [[ M = #ROW + 1 ]]
  #ROW 1+ 1
  DO
    [[ M1 = M - I ]]
    {{ X ( M1 ) = Y ( M1 ) / AA ( M1 M1 ) }}
    [[ ITEMP = ( M - I ) - #ROW ]]
    ITEMP DUP
    0>
    IF ." Trival! " QUIT THEN
      0<
      IF #ROW 1+ 1
        DO
        [[ M2 = M - I ]]
        {{ X ( M1 ) = X ( M1 ) - AA ( M1 M2 ) * X ( M2 ) / AA ( M1 M1 ) }}
        [[ ITEMP = ( M - I - 1 ) - ( M - J ) ]]
        ITEMP
        0=
        IF LEAVE THEN
        LOOP
      THEN
  LOOP ;
: 輸出答案
  #ROW 1+ 1
  DO
   {{ FTEMP = X ( I ) }}
   CR ." X( " I . ." )= " FTEMP F.
  LOOP ;
: GAUSS
  輸入數據
  高斯處理
  得到解答
  輸出答案
; 

沒有留言: