一百個例題 (41 ~ 45)
Ching-Tang Tseng
Hamilton, New Zealand
9 September 2024
FORTH:慶的世界。 本網頁為中華民國(Republic Of China)FORTH語言協會 永久性網頁
2024年9月8日 星期日
第(41)個範例程式是一個根據兩個點的座標值，換算出通解為代表一條直線的方程式。
平面幾何的五大公設中第一條就是：任意一點到另外任意一點可以畫出一直线。這個範例就是在實踐出這條公設的程式。
口語問題想寫成程式，必須要先靠代數，而代數就是解方程式的科學，好的程式語言其語法才能接近代數，也就是方程式怎麼解？程式就該怎麼寫。我拿這種觀念，展示我創作之程式語言的語法，範例的主要目的在此。
我念大四時，台灣出版教科書的規矩，在法律上正要起大變化，經濟也正要起飛，可以配得上版權規矩會令書變得很貴的發展。那時，市面上還流通著許多盜版的書籍。我有個清大動機研究所的朋友，勸我趁早買一套大學畢業後才適合去仔細研讀的英文原版好書：『數學的內容方法及意義』，A. D. Aleksandrov, A. N. Kolmogorov, M. A. Lavrent'ev，“Mathematics: Its Content, Methods and Meaning”，Dover Publications，July 7, 1999，Total pages: 1120。他從台大機械系畢業，先考上了清大研究所，但先當兵，幹排長時，我們在馬祖同寢室，一起帶兵。後來，我接受他的建議，買了三本這一套盜版書。直到今天，此書仍要花錢，才能從網上獲得，上列資料就是今天才從網上節錄下來的售書廣告，確實是一本值得終生保留的好書。我在設計這個範例程式時，就受到這本書的啟發，我讀的是當年台灣清大請書商盜版印售，無出版社的書，此書把所有的數學知識，分章寫成濃縮的敘述，其中都是精華。例如：微分與積分都只被寫出五種基本方法，就憑這各五種方法，保證您可以推演出所有一般大一所教過之微積分的公式。我看到這些內容，受到很大的啟示，也體會出這套書不適合大一理工科學生閱讀的道理，卻非常適合理工科系畢業的學生，用來總結自己大學所學。
平面空間有無數個座標點，兩點連成一線的方程式是確定的公設，高中生的數學，會要求學生學會如何手算算出此直線方程式各項的係數，問題是個死問題，算法也是個死計算，只是，有時會碰到特例，兩點會重合或只是退化成水平或者垂直的直線。當年，上高中時，我們確實是耗費不少時間在手算這種數學的問題上，那時沒有電腦，連計算器也沒有，只能挨磨、挨練。
後來，我自己設計了數學運算系統，有了工具之後，覺得程式語言的主要功能，就是能解決冗長的計算問題，它能在學生學會數學道理之後，不要再浪費太多的時間去搞死算，就能在數學的學習上更快速的進展。更甚的是，程式語言最好也能配合得上繪圖，快速直接地讓結果呈現在電腦上繪出圖來，學習記憶就會保證更為強烈，學習就更容易。這些事情，我設計的系統都能快速的辦到，很多商售的數學套裝軟體，後來也都走上了這樣的發展方式，只是，它們都要賣錢。
範例很簡單，我在讀過很多書後，才有被書籍啟發後的體會，將其寫成程式，程式是整理了書中的道理才寫成的，從口語到代數再到此範例程式，有它的過程，也能永久使用。
不要小看這麼一個簡單的範例，它是以正向思考的方式寫成的程式，所以程式很好寫。我卻拿過它來反向思考，寫過兩點間繪出一條直線的指令，叫 LineTo ( x1 y1 x2 y2 — )，就給兩點的座標，在螢幕上繪出直線，它不容易設計。
有人想了一套方法，並向美國專利局申請獲得過專利，我用同樣想法試出了他的問題：只能在第一象限有效，其餘無效。美國政府卻還真給了這傢伙專利權，很可笑。我不想網上發難給人難堪，技術自己留著便是。我得花十頁以上的份量，才能把正確修正出來的萬用程式交代得清楚，但這種已成為專利的題目，不能作為範例教材。
:
平面幾何的五大公設中第一條就是：任意一點到另外任意一點可以畫出一直线。這個範例就是在實踐出這條公設的程式。
口語問題想寫成程式，必須要先靠代數，而代數就是解方程式的科學，好的程式語言其語法才能接近代數，也就是方程式怎麼解？程式就該怎麼寫。我拿這種觀念，展示我創作之程式語言的語法，範例的主要目的在此。
我念大四時，台灣出版教科書的規矩，在法律上正要起大變化，經濟也正要起飛，可以配得上版權規矩會令書變得很貴的發展。那時，市面上還流通著許多盜版的書籍。我有個清大動機研究所的朋友，勸我趁早買一套大學畢業後才適合去仔細研讀的英文原版好書：『數學的內容方法及意義』，A. D. Aleksandrov, A. N. Kolmogorov, M. A. Lavrent'ev，“Mathematics: Its Content, Methods and Meaning”，Dover Publications，July 7, 1999，Total pages: 1120。他從台大機械系畢業，先考上了清大研究所，但先當兵，幹排長時，我們在馬祖同寢室，一起帶兵。後來，我接受他的建議，買了三本這一套盜版書。直到今天，此書仍要花錢，才能從網上獲得，上列資料就是今天才從網上節錄下來的售書廣告，確實是一本值得終生保留的好書。我在設計這個範例程式時，就受到這本書的啟發，我讀的是當年台灣清大請書商盜版印售，無出版社的書，此書把所有的數學知識，分章寫成濃縮的敘述，其中都是精華。例如：微分與積分都只被寫出五種基本方法，就憑這各五種方法，保證您可以推演出所有一般大一所教過之微積分的公式。我看到這些內容，受到很大的啟示，也體會出這套書不適合大一理工科學生閱讀的道理，卻非常適合理工科系畢業的學生，用來總結自己大學所學。
平面空間有無數個座標點，兩點連成一線的方程式是確定的公設，高中生的數學，會要求學生學會如何手算算出此直線方程式各項的係數，問題是個死問題，算法也是個死計算，只是，有時會碰到特例，兩點會重合或只是退化成水平或者垂直的直線。當年，上高中時，我們確實是耗費不少時間在手算這種數學的問題上，那時沒有電腦，連計算器也沒有，只能挨磨、挨練。
後來，我自己設計了數學運算系統，有了工具之後，覺得程式語言的主要功能，就是能解決冗長的計算問題，它能在學生學會數學道理之後，不要再浪費太多的時間去搞死算，就能在數學的學習上更快速的進展。更甚的是，程式語言最好也能配合得上繪圖，快速直接地讓結果呈現在電腦上繪出圖來，學習記憶就會保證更為強烈，學習就更容易。這些事情，我設計的系統都能快速的辦到，很多商售的數學套裝軟體，後來也都走上了這樣的發展方式，只是，它們都要賣錢。
範例很簡單，我在讀過很多書後，才有被書籍啟發後的體會，將其寫成程式，程式是整理了書中的道理才寫成的，從口語到代數再到此範例程式，有它的過程，也能永久使用。
不要小看這麼一個簡單的範例，它是以正向思考的方式寫成的程式，所以程式很好寫。我卻拿過它來反向思考，寫過兩點間繪出一條直線的指令，叫 LineTo ( x1 y1 x2 y2 — )，就給兩點的座標，在螢幕上繪出直線，它不容易設計。
有人想了一套方法，並向美國專利局申請獲得過專利，我用同樣想法試出了他的問題：只能在第一象限有效，其餘無效。美國政府卻還真給了這傢伙專利權，很可笑。我不想網上發難給人難堪，技術自己留著便是。我得花十頁以上的份量，才能把正確修正出來的萬用程式交代得清楚，但這種已成為專利的題目，不能作為範例教材。
:
\ (41)LineEquation.f
\ 使用了=/=符號，因此只適用於V.658以後版本。
((
Theme : Line equation
Author : Ching-Tang Tseng, Hamilton, New Zealand
Date : 2015-03-20
Reference : None
Description :
Equation of a line determined by two points
y = ax+b
a = (y2-y1)/(x2-x1)
b = y1-((y2-y1)*x1/(x2-x1))
yi = y1+(xi-x1)*(y2-y1)/(x2-x1)
))
8 Reals X1 Y1 X2 Y2 Xi Yi a b
: Reset
{{ X1 = 0 }} {{ Y1 = 0 }}
{{ X2 = 0 }} {{ Y2 = 0 }}
{{ Xi = 0 }} {{ Yi = 0 }}
{{ a = 0 }} {{ b = 0 }}
;
: QAEnter
BASIC
10 RUN PAGE CR
20 PRINT " Equation of a line determined by coordinates of two points. "
30 PRINT " Then, enter a new Xi, correspondent Yi can be computed out. "
40 RUN CR CR
50 PRINT " Enter 1st point coordinates X1 Y1 "
60 InputR X1 , Y1
70 PRINT " Enter 2nd point coordinates X2 Y2 "
80 InputR X2 , Y2
90 END
;
: ExceptionTreatment
BASIC
10 IF { ( X1 = X2 ) AND ( Y1 = Y2 ) } THEN 80
20 IF { X1 = X2 } THEN 40
30 IF { Y1 = Y2 } THEN 60
40 PRINT { " This is a vertical line. Line Equation: x = " ; X1 }
50 GOTO 90
60 PRINT { " This is a horizontal line. Line Equation: y = " ; Y1 }
70 GOTO 90
80 PRINT " Enter data are the same point, program stopped! "
90 END
;
: LineEquation
BASIC
10 LET { a = ( Y2 - Y1 ) / ( X2 - X1 ) }
20 LET { b = Y1 - ( ( Y2 - Y1 ) * X1 / ( X2 - X1 ) ) }
30 IF { b >= 0 } THEN 60
40 PRINT CR { " Line Equation: y = " ; a ; " x - " ; abs ( b ) }
50 GOTO 70
60 PRINT CR { " Line Equation: y = " ; a ; " x + " ; b }
70 END
;
Integer K
: hi
BASIC
10 RUN Reset
20 RUN QAEnter
30 IF { ( X1 =/= X2 ) AND ( Y1 =/= Y2 ) } THEN 60
40 RUN ExceptionTreatment
50 GOTO 150
60 RUN LineEquation
70 PRINT " Enter Xi "
80 InputR Xi
90 LET { Yi = Y1 + ( Xi - X1 ) * ( Y2 - Y1 ) / ( X2 - X1 ) }
100 PRINT { " Xi = " ; Xi ; " Yi = " ; Yi }
110 PRINT CR " Repeat(y/Y)? "
120 LET K = KEY
130 RUN K EMIT
140 IF ( K = 89 ) OR ( K = 121 ) THEN -70
150 END
;
cr cr
.( Usage: hi )
第(42)個範例，是一個標準、完整、最簡單的數據調適(fitting)技術應用範例。
程式仍然是在求一條直線方程式的係數，範例(41)也是。採用直線方程式為例，純粹只是取其簡單、易懂、適合當例子解說而已。同樣的技術可以應用在非常複雜的方程式上，這套調適技術稱為最小平方差的調適技術(the least squares curve fitting)。
前一範例，單憑『兩點』座標值，求出直線方程式的係數。這一範例，可憑『無限量』組座標值，根據統計分析的技術，匹配出最理想之直線方程式的係數。這是兩者間最大的區別。
為什麼要學曲線(包括直線)調適(curve fitting)技術？如果您念過研究所，或在正式期刊上發表過帶有實驗數據的論文，您一定知道，論文中的數據，不可以表列實驗所得的裸值就算了。實驗數據通常都代表著變化趨勢，而變化趨勢都可以寫成數學方程式。因此，在論文中，作者必須好好的把實驗數據匹配出一個最適當的方程式，才能寫進論文，方程式中係數的決定方法，就得使用曲線調適技術。
實驗次數，通常的要求是最少也得能對不同點實驗3次以上。為了取得不同條件位置之值，實驗次數還希望越多越好。那麼，如果代表變化趨勢的方程式是一條直線，多於兩點以上的數據，要怎麼處理？這就得用到統計分析技術了，最基礎與最有名的方法，就是最小平方差(least squares error)的分析方法。
統計學上告訴我們最小平方差的分析方法很可靠，學術與實務上都普遍的能被接受，此處不講它的原理，只用它的公式，範例中列有所用公式，不難懂。
這個範例，我在幾十年前，就常用它來測試系統，資料來源已列在範例中。我從自己設計固定點小數的實數計算系統開始，就使用此例測試系統了，早期的文章都發佈在協會的月會上，最早能發現的網文，則大約是 2008 年我開始使用 Google 開闢個人網頁以後的事情。由此可見，範例的有用程度。
剛開始使用此範例時，我強調的是計算次數非常冗長，純用 Forth 格式寫程式，麻煩又痛苦，若強要我用，我寧可不用。於是，我設計了以 BASIC 格式寫程式的技術，處裡這種冗長計算的程式，輕鬆又愉快，能用，為什麼不用？它就是 ABC Forth 必須發展出來的真正理由。
我的系統初次完整公開的日期是 2006 年，至今已 18 年。
早在 1980 年代，我們就根據 Forth 系統發明人的技術，設計過具有中算符(infix)功能的系統，可以改善計算程式不好寫的問題。但長期使用下來，總覺得只具有中算符的功能，性能總是仍有欠缺，我也不太使用。
後來，國際上在 Forth 界很有名氣的 J. V. Noble 公開捐獻了一套可在 Forth 中以 FORTRAN 格式寫程式的技術，這種方法曾經在美國流行，被廣泛使用過一陣子，最後也無疾而終了。考其原因，主要還是 FORTRAN 不如 BASIC 的那麼能被群眾接受，只有理工科系出身者比較能接受 FORTRAN 。
因次，我便著手單搞以 BASIC 格式寫程式的技術，每次發展出胚胎系統時，都首先使用這個範例測試系統。
如果我創作之以 BASIC 格式設計程式的方式不是個很理想的東西，那麼，我在網上公開貼出過的網文，就根本不會有人注意。但實際情況是， comp.lang.forth 國際論壇網頁上，年年都有人提及我創作的 ABC Forth 並進行討論，也有人跟進，搞成過最原始的 BASIC compiler ，但他們都只能跑整數，沒有我設計出來之 BASIC 能叫用 BASIC 的功能，還有許多其他欠缺，我的創作就恆有存在的價值了。至少，我天天用它來解決問題，幾十年如一日。
這個範例程式發展之後，我也理出了統一的最終應用指令規格：
當只執行一次性的測試時，數據就被放在程式中，使用的指令是 test ，想更換數據的人，可以自行在程式中更換後，重新載入。
當想執行對答式輸入數據的測試時，就使用 hi 指令。
當數據多到應該使用另一個檔案來儲存後才輸入系統，以避免數據出錯時，就使用 main 指令。
這是一個很原始也很好用的範例程式。 20150413 本網頁刊出的網文 Curve fitting demo in ABC FORTH 則是採用差分表推演出來的曲線調適展示，方法與此程式有別。
:
程式仍然是在求一條直線方程式的係數，範例(41)也是。採用直線方程式為例，純粹只是取其簡單、易懂、適合當例子解說而已。同樣的技術可以應用在非常複雜的方程式上，這套調適技術稱為最小平方差的調適技術(the least squares curve fitting)。
前一範例，單憑『兩點』座標值，求出直線方程式的係數。這一範例，可憑『無限量』組座標值，根據統計分析的技術，匹配出最理想之直線方程式的係數。這是兩者間最大的區別。
為什麼要學曲線(包括直線)調適(curve fitting)技術？如果您念過研究所，或在正式期刊上發表過帶有實驗數據的論文，您一定知道，論文中的數據，不可以表列實驗所得的裸值就算了。實驗數據通常都代表著變化趨勢，而變化趨勢都可以寫成數學方程式。因此，在論文中，作者必須好好的把實驗數據匹配出一個最適當的方程式，才能寫進論文，方程式中係數的決定方法，就得使用曲線調適技術。
實驗次數，通常的要求是最少也得能對不同點實驗3次以上。為了取得不同條件位置之值，實驗次數還希望越多越好。那麼，如果代表變化趨勢的方程式是一條直線，多於兩點以上的數據，要怎麼處理？這就得用到統計分析技術了，最基礎與最有名的方法，就是最小平方差(least squares error)的分析方法。
統計學上告訴我們最小平方差的分析方法很可靠，學術與實務上都普遍的能被接受，此處不講它的原理，只用它的公式，範例中列有所用公式，不難懂。
這個範例，我在幾十年前，就常用它來測試系統，資料來源已列在範例中。我從自己設計固定點小數的實數計算系統開始，就使用此例測試系統了，早期的文章都發佈在協會的月會上，最早能發現的網文，則大約是 2008 年我開始使用 Google 開闢個人網頁以後的事情。由此可見，範例的有用程度。
剛開始使用此範例時，我強調的是計算次數非常冗長，純用 Forth 格式寫程式，麻煩又痛苦，若強要我用，我寧可不用。於是，我設計了以 BASIC 格式寫程式的技術，處裡這種冗長計算的程式，輕鬆又愉快，能用，為什麼不用？它就是 ABC Forth 必須發展出來的真正理由。
我的系統初次完整公開的日期是 2006 年，至今已 18 年。
早在 1980 年代，我們就根據 Forth 系統發明人的技術，設計過具有中算符(infix)功能的系統，可以改善計算程式不好寫的問題。但長期使用下來，總覺得只具有中算符的功能，性能總是仍有欠缺，我也不太使用。
後來，國際上在 Forth 界很有名氣的 J. V. Noble 公開捐獻了一套可在 Forth 中以 FORTRAN 格式寫程式的技術，這種方法曾經在美國流行，被廣泛使用過一陣子，最後也無疾而終了。考其原因，主要還是 FORTRAN 不如 BASIC 的那麼能被群眾接受，只有理工科系出身者比較能接受 FORTRAN 。
因次，我便著手單搞以 BASIC 格式寫程式的技術，每次發展出胚胎系統時，都首先使用這個範例測試系統。
如果我創作之以 BASIC 格式設計程式的方式不是個很理想的東西，那麼，我在網上公開貼出過的網文，就根本不會有人注意。但實際情況是， comp.lang.forth 國際論壇網頁上，年年都有人提及我創作的 ABC Forth 並進行討論，也有人跟進，搞成過最原始的 BASIC compiler ，但他們都只能跑整數，沒有我設計出來之 BASIC 能叫用 BASIC 的功能，還有許多其他欠缺，我的創作就恆有存在的價值了。至少，我天天用它來解決問題，幾十年如一日。
這個範例程式發展之後，我也理出了統一的最終應用指令規格：
當只執行一次性的測試時，數據就被放在程式中，使用的指令是 test ，想更換數據的人，可以自行在程式中更換後，重新載入。
當想執行對答式輸入數據的測試時，就使用 hi 指令。
當數據多到應該使用另一個檔案來儲存後才輸入系統，以避免數據出錯時，就使用 main 指令。
這是一個很原始也很好用的範例程式。 20150413 本網頁刊出的網文 Curve fitting demo in ABC FORTH 則是採用差分表推演出來的曲線調適展示，方法與此程式有別。
:
\ (42)LinearLeastSquareCurveFit.f
\ 最小平方直線擬合(調適)，利用公式直接計算求得答案。
((
Theme : First-Order Least Squares Curve Fitting
Author : Ching-Tang Tseng, Hamilton, New Zealand
Date : 2015-01-25
Reference : F.R. Ruckdeschel, BASIC Scientific Subroutines Vol.II
BYTE/McGraw-Hill, 1980, p.16
Description :
This program calculates a linear least squares fit to a given data set.
y = a + b x, given (Xi,Yi), i>1, get a, b. i>2, Standard Deviation D as well.
Hint: using formula to do the tedious computation.
a = ( Sum(Xi)*Sum(Yi))-N*(Sum(Xi*Yi) )/( (Sum(Xi))^2-N*Sum((Xi)^2) )
b = ( Sum(Xi)*Sum(Xi*Yi)-Sum((Xi)^2)*Sum(Yi) )/( (Sum(Xi))^2-N*Sum((Xi)^2) )
D = SQRT ( Sum(Yi-a-b*Xi)/(N-2) )
))
20 value Size \ Unique setting in this program for maximum matrix Dimension
Size ARRAY XX
Size ARRAY YY
10 Reals A B D D1 A1 A2 B0 B1 X Y
2 Integers N I
: ResetArray BASIC
10 FOR I = 1 TO Size
20 LET { XX ( I ) = 0 }
30 LET { YY ( I ) = 0 }
40 NEXT I
50 END
;
: ResetAll
ResetArray
{{ A1 = 0 }} {{ A2 = 0 }}
{{ B0 = 0 }} {{ B1 = 0 }}
;
: FixedData
[[ N = 11 ]]
{{ XX ( 1 ) = 0 }} {{ YY ( 1 ) = 0 }}
{{ XX ( 2 ) = 0.157 }} {{ YY ( 2 ) = 0.1563558123 }}
{{ XX ( 3 ) = 0.314 }} {{ YY ( 3 ) = 0.3088655201 }}
{{ XX ( 4 ) = 0.471 }} {{ YY ( 4 ) = 0.4537776271 }}
{{ XX ( 5 ) = 0.628 }} {{ YY ( 5 ) = 0.5875275257 }}
{{ XX ( 6 ) = 0.785 }} {{ YY ( 6 ) = 0.7068251811 }}
{{ XX ( 7 ) = 0.942 }} {{ YY ( 7 ) = 0.8087360606 }}
{{ XX ( 8 ) = 1.099 }} {{ YY ( 8 ) = 0.8907533184 }}
{{ XX ( 9 ) = 1.256 }} {{ YY ( 9 ) = 0.9508594605 }}
{{ XX ( 10 ) = 1.413 }} {{ YY ( 10 ) = 0.9875759713 }}
{{ XX ( 11 ) = 1.550 }} {{ YY ( 11 ) = 0.9997837642 }}
;
: Coefficient
BASIC
10 FOR I = 1 TO N
20 LET { A1 = A1 + XX ( I ) }
30 LET { A2 = A2 + ( XX ( I ) * XX ( I ) ) }
40 LET { B0 = B0 + YY ( I ) }
50 LET { B1 = B1 + ( YY ( I ) * XX ( I ) ) }
60 NEXT I
70 LET { D = ( A1 * A1 ) - ( I>R ( N ) * A2 ) }
80 LET { A = ( A1 * B1 ) - ( A2 * B0 ) }
90 LET { A = A / D }
100 LET { B = ( A1 * B0 ) - ( I>R ( N ) * B1 ) }
110 LET { B = B / D }
120 END
;
: StandardDeviation
BASIC
10 LET { D = 0 }
20 FOR I = 1 TO N
30 LET { D1 = YY ( I ) - A - ( B * XX ( I ) ) }
40 LET { D = D + ( D1 * D1 ) }
50 NEXT I
60 LET { D = SQRT ( D / ( I>R ( N ) - 2 ) ) }
70 END
;
: Outcome
BASIC
10 run CR
20 Print " (1)Input data as following : "
30 for i = 1 to N
40 print i ; { XX ( i ) , YY ( i ) }
50 next i
60 run CR
70 PRINT " (2)Fitted equation is : "
80 IF { B >= 0 } THEN 110
90 PRINT { " y = " ; A ; B ; " x " }
100 GOTO 120
110 PRINT { " y = " ; A ; " + " ; B ; " x " }
120 RUN CR
130 PRINT { " (3)Standard deviation of fit is : " ; D }
140 RUN CR
150 END
;
: Once
BASIC
10 run Coefficient
20 run StandardDeviation
30 run Outcome
40 end
;
: test
BASIC
10 RUN ResetAll
20 RUN FixedData
30 RUN Once
40 END
;
: QAData ( -- )
BASIC
10 PRINT " How many sets of Xi Yi given data will be input "
20 InputI N
30 RUN CR
40 FOR I = 1 TO N
50 PRINT " Enter X( " ; I ; " ) Y( " ; I ; " ) = "
60 InputR X , Y
70 LET { XX ( I ) = X } :: { YY ( I ) = Y }
80 NEXT I
90 END
;
: hi
BASIC
10 run ResetAll
20 run QAData
30 run Once
40 end
;
: FileData ( -- )
BASIC
10 RUN S" (41-1)Data.f" Get-File
20 RUN GetOneLineData
30 LET N = INT ( BRA ( 1 ) )
40 FOR I = 1 TO N
50 RUN GetOneLineData
60 LET { XX ( I ) = BRA ( 2 ) } :: { YY ( I ) = BRA ( 3 ) }
70 NEXT I
80 END
;
: main ( -- )
BASIC
10 run ResetAll
20 RUN FileData
30 RUN Once
40 END
;
page
.( Usage: ) cr
.( 1. test : for fixed input data ) cr
.( 2. hi : for interactive input data ) cr
.( 3. main : for data input from file ) cr
\S
\ (42-1)Data.f 數據檔案的內容
11 \ This line means following are 11 sets Xi Yi data.
1 0 1 \ This line means set 1 (Xi,Yi) are (0,1)
2 0.157 0.9918106 \ This line means set 2 (Xi,Yi) are (0.157,0.9918106)
3 0.314 0.96756368
4 0.471 0.92820596
5 0.628 0.87526015
6 0.785 0.8107458
7 0.942 0.73707522
8 1.099 0.65693061
9 1.256 0.57313037
10 1.413 0.48849118
11 1.57 0.40569572
第(43)個範例，是第(42)個範例的延伸，原為一階的直線方程式，延伸為二階的曲線方程式，求解調適問題時，由這個範例展示其最小平方差的曲線調適程式寫法。
程式後面，列有執行出來的結果，供作參考，這個範例也是一個完整可以廣泛應用的程式。
原始資料中的程式，有兩處迴路指標的使用，採用了從 0 做到 N-1 的設計方法，我自覺可以另用寫法，故改掉了。有很多系統，例如： Metlab 系統，就不允許宣告出來的陣列或矩陣資料結構，使用以 0 當指標來放置數據，因此，我刻意將原始程式修正為從 1 做到 N 。
ABC Forth 系統，就是在經常遇到這種問題後，才擴充了資料結構的設計。當然，從此之後，每次用到陣列、矩陣、張量等宣告時，都必定會自動產生可用 0 作為指標的記憶體保留量，不用時，確實是有點浪費。但在現今的世界，節省記憶體用量，已不再是系統設計者應該強調的重點，給予使用者最大的方便，才是設計系統一開始就該規劃出來的功能。簡而言之，我設計的系統，允許使用者使用 0 作為指標。 Forth 運作的哲理，給了我們很大的使用彈性，您若高興，就算搞出虛擬的負數指標，又有什麼不可以？只要您想得清楚，這樣搞會有用，那就儘管去搞。
我們在抄改寫前人設計出來的寶貴程式資源時，就是在做著技術傳承的工作，我把系統設計成能適應於各種狀況的應用時，是遭遇過許多實例以及長久的磨練後，總結產生的結果。這些改良性的設計，都已經跟最原始的設計有了很大的不同。
目前世界上的 Forth 系統設計名家，在見過我的貼文後，也極力想實現同樣的夢想，把能用 BASIC 格式寫程式的方法，也納入他們設計的 Forth 系統內， gForth , iForth , ciForth 的作者都做過同樣的努力，在隨系統提供給使用者的資料夾內，都能發現可擴充的類似源程式。
但是，大家都沒有下功夫仔細研究當初 Charles H. Moore 設計之 Tiny BASIC Compiler 的主要精神，所以設計不出像我這樣多功能的產品來，全都只是把原作品改成在他們的 Forth 系統內可以 run 就算了，如此簡單的抄襲，當然不會有用。實際上， Tiny BASIC Compiler 程式裡面蘊藏了龐大、複雜的技術，能夠繼續發揮，必定驚人。
我做過了 Tiny BASIC Compiler 原始程式的徹底研究，融會貫通之後，也才知道新時代的新觀念，該如何加進系統，這個可以使用 0 當指標的功能，也是其中之一。
這些曲線調適程式，在程式內所論及的資源中都有，我很輕鬆的就能改寫出來，需要時才會去做，若要把它們都納入範例，那就太多了。我只保留了這些書籍，唾手可得。
從 y = a0 + a1X 進展到
y = a0 + a1X + a2X^2，也能進展到：
Y = a0 + a1X + a2X^2 + ……. + anX^n
更甚的是，還有推衍出來的應用，例如：下列這些方程式，我都有現成資源可以抄改寫來用。
幾何方程式： y = a * X^b
指數方程式： y = a * exp(bX)
修正指數方程式： y = K - a*b^X
還有許多更複雜的方程式，我就不一一列舉了，反正留者很舊的書，書中都有很容易被我抄改寫的現成程式，這些都是我個人的技術財富，只要我擁有 ABC Forth 系統，書本財富就有價值。
:
程式後面，列有執行出來的結果，供作參考，這個範例也是一個完整可以廣泛應用的程式。
原始資料中的程式，有兩處迴路指標的使用，採用了從 0 做到 N-1 的設計方法，我自覺可以另用寫法，故改掉了。有很多系統，例如： Metlab 系統，就不允許宣告出來的陣列或矩陣資料結構，使用以 0 當指標來放置數據，因此，我刻意將原始程式修正為從 1 做到 N 。
ABC Forth 系統，就是在經常遇到這種問題後，才擴充了資料結構的設計。當然，從此之後，每次用到陣列、矩陣、張量等宣告時，都必定會自動產生可用 0 作為指標的記憶體保留量，不用時，確實是有點浪費。但在現今的世界，節省記憶體用量，已不再是系統設計者應該強調的重點，給予使用者最大的方便，才是設計系統一開始就該規劃出來的功能。簡而言之，我設計的系統，允許使用者使用 0 作為指標。 Forth 運作的哲理，給了我們很大的使用彈性，您若高興，就算搞出虛擬的負數指標，又有什麼不可以？只要您想得清楚，這樣搞會有用，那就儘管去搞。
我們在抄改寫前人設計出來的寶貴程式資源時，就是在做著技術傳承的工作，我把系統設計成能適應於各種狀況的應用時，是遭遇過許多實例以及長久的磨練後，總結產生的結果。這些改良性的設計，都已經跟最原始的設計有了很大的不同。
目前世界上的 Forth 系統設計名家，在見過我的貼文後，也極力想實現同樣的夢想，把能用 BASIC 格式寫程式的方法，也納入他們設計的 Forth 系統內， gForth , iForth , ciForth 的作者都做過同樣的努力，在隨系統提供給使用者的資料夾內，都能發現可擴充的類似源程式。
但是，大家都沒有下功夫仔細研究當初 Charles H. Moore 設計之 Tiny BASIC Compiler 的主要精神，所以設計不出像我這樣多功能的產品來，全都只是把原作品改成在他們的 Forth 系統內可以 run 就算了，如此簡單的抄襲，當然不會有用。實際上， Tiny BASIC Compiler 程式裡面蘊藏了龐大、複雜的技術，能夠繼續發揮，必定驚人。
我做過了 Tiny BASIC Compiler 原始程式的徹底研究，融會貫通之後，也才知道新時代的新觀念，該如何加進系統，這個可以使用 0 當指標的功能，也是其中之一。
這些曲線調適程式，在程式內所論及的資源中都有，我很輕鬆的就能改寫出來，需要時才會去做，若要把它們都納入範例，那就太多了。我只保留了這些書籍，唾手可得。
從 y = a0 + a1X 進展到
y = a0 + a1X + a2X^2，也能進展到：
Y = a0 + a1X + a2X^2 + ……. + anX^n
更甚的是，還有推衍出來的應用，例如：下列這些方程式，我都有現成資源可以抄改寫來用。
幾何方程式： y = a * X^b
指數方程式： y = a * exp(bX)
修正指數方程式： y = K - a*b^X
還有許多更複雜的方程式，我就不一一列舉了，反正留者很舊的書，書中都有很容易被我抄改寫的現成程式，這些都是我個人的技術財富，只要我擁有 ABC Forth 系統，書本財富就有價值。
:
\ (43-1) Data.f 輸入數據檔案的內容
11 \ This line means following are 11 sets Xi Yi data.
1 0 1 \ This line means set 1 (Xi,Yi) are (0,1)
2 0.157 0.9918106 \ This line means set 2 (Xi,Yi) are (0.157, 0.9918106)
3 0.314 0.96756368
4 0.471 0.92820596
5 0.628 0.87526015
6 0.785 0.8107458
7 0.942 0.73707522
8 1.099 0.65693061
9 1.256 0.57313037
10 1.413 0.48849118
11 1.57 0.40569572
\ (43)2ndOrderLeastSquareCurveFit.f
\ 二階最小乘方曲線擬合(調適)
\ 利用公式直接計算求得答案。
((
Theme : First-Order Least Squares Curve Fitting
Author : Ching-Tang Tseng, Hamilton, New Zealand
Date : 2015-01-25
Reference : F.R. Ruckdeschel, BASIC Scientific Subroutines Vol.II
BYTE/McGraw-Hill, 1980, p.24
))
\ Size is the unique setting in this program for maximum matrix Dimension
20 value Size
Size ARRAY XX
Size ARRAY YY
16 Reals A B C D D1 A0 A1 A2 A3 A4 B0 B1 B2 X Y NN
3 Integers M N I
: ResetArray BASIC
10 FOR I = 1 TO Size
20 LET { XX ( I ) = 0 }
30 LET { YY ( I ) = 0 }
40 NEXT I
50 END
;
: ResetAll
ResetArray
{{ A0 = 1 }} {{ A1 = 0 }} {{ A2 = 0 }} {{ A3 = 0 }} {{ A4 = 0 }}
{{ B0 = 0 }} {{ B1 = 0 }} {{ B2 = 0 }}
;
: FixedData
[[ N = 11 ]]
{{ XX ( 1 ) = 0 }} {{ YY ( 1 ) = 0 }}
{{ XX ( 2 ) = 0.157 }} {{ YY ( 2 ) = 0.1563558123 }}
{{ XX ( 3 ) = 0.314 }} {{ YY ( 3 ) = 0.3088655201 }}
{{ XX ( 4 ) = 0.471 }} {{ YY ( 4 ) = 0.4537776271 }}
{{ XX ( 5 ) = 0.628 }} {{ YY ( 5 ) = 0.5875275257 }}
{{ XX ( 6 ) = 0.785 }} {{ YY ( 6 ) = 0.7068251811 }}
{{ XX ( 7 ) = 0.942 }} {{ YY ( 7 ) = 0.8087360606 }}
{{ XX ( 8 ) = 1.099 }} {{ YY ( 8 ) = 0.8907533184 }}
{{ XX ( 9 ) = 1.256 }} {{ YY ( 9 ) = 0.9508594605 }}
{{ XX ( 10 ) = 1.413 }} {{ YY ( 10 ) = 0.9875759713 }}
{{ XX ( 11 ) = 1.550 }} {{ YY ( 11 ) = 0.9997837642 }}
;
: Coefficient
BASIC
10 FOR M = 1 TO N \ modified from original 0 --> N-1
20 LET { A1 = A1 + XX ( M ) }
30 LET { A2 = A2 + XX ( M ) * XX ( M ) }
40 LET { A3 = A3 + XX ( M ) * XX ( M ) * XX ( M ) }
50 LET { A4 = A4 + XX ( M ) * XX ( M ) * XX ( M ) * XX ( M ) }
60 LET { B0 = B0 + YY ( M ) }
70 LET { B1 = B1 + YY ( M ) * XX ( M ) }
80 LET { B2 = B2 + YY ( M ) * XX ( M ) * XX ( M ) }
90 NEXT M
100 LET { NN = I>R ( N ) }
110 LET { A1 = A1 / NN }
120 LET { A2 = A2 / NN }
130 LET { A3 = A3 / NN }
140 LET { A4 = A4 / NN }
150 LET { B0 = B0 / NN }
160 LET { B1 = B1 / NN }
170 LET { B2 = B2 / NN }
180 LET { D = A0 * ( A2 * A4 - A3 * A3 ) - A1 * ( A1 * A4 - A3 * A2 ) + A2 * ( A1 * A3 - A2 * A2 ) }
190 LET { A = B0 * ( A2 * A4 - A3 * A3 ) + B1 * ( A3 * A2 - A1 * A4 ) + B2 * ( A1 * A3 - A2 * A2 ) }
200 LET { A = A / D }
210 LET { B = B0 * ( A3 * A2 - A1 * A4 ) + B1 * ( A0 * A4 - A2 * A2 ) + B2 * ( A2 * A1 - A0 * A3 ) }
220 LET { B = B / D }
230 LET { C = B0 * ( A1 * A3 - A2 * A2 ) + B1 * ( A1 * A2 - A0 * A3 ) + B2 * ( A0 * A2 - A1 * A1 ) }
240 LET { C = C / D }
250 END
;
: StandardDeviation
BASIC
10 LET { D = 0 }
20 FOR M = 1 TO N \ modified from original 0 --> N-1
30 LET { D1 = YY ( M ) - A - B * XX ( M ) - C * XX ( M ) * XX ( M ) }
40 LET { D = D + D1 * D1 }
50 NEXT M
60 LET { D = SQRT ( D / ( NN - 3 ) ) }
70 END
;
: Outcome
BASIC
10 run CR
20 Print " (1)Input data as following : "
30 for i = 1 to N
40 print i ; { XX ( i ) , YY ( i ) }
50 next i
60 run CR
70 PRINT " (2)Fitted equation is : " CR
\ smart print
80 run ." y = " A G.
90 if { B >= 0 } then 120
100 run B G. ." x "
110 goto 130
120 run ." + " B G. ." x "
130 if { C >= 0 } then 160
140 run C G. ." x^2 " CR
150 goto 170
160 run ." + " C G. ." x^2 " CR
\ lazy print
\ 80 PRINT { " y =( " ; A ; " )+( " ; B ; " )x+( " ; C ; " )x^2 " CR }
170 PRINT { " (3)Standard deviation of fit is : " ; D }
180 RUN CR
190 END
;
: Once
BASIC
10 run Coefficient
20 run StandardDeviation
30 run Outcome
40 end
;
: test
BASIC
10 RUN ResetAll
20 RUN FixedData
30 RUN Once
40 END
;
: QAData ( -- )
BASIC
10 PRINT " How many sets of Xi Yi given data will be input "
20 InputI N
30 RUN CR
40 FOR I = 1 TO N
50 PRINT " Enter X( " ; I ; " ) Y( " ; I ; " ) = "
60 InputR X , Y
70 LET { XX ( I ) = X } :: { YY ( I ) = Y }
80 NEXT I
90 END
;
: hi
BASIC
10 run ResetAll
20 run QAData
30 run Once
40 end
;
: FileData ( -- )
BASIC
10 RUN S" (42-1)Data.f" Get-File
20 RUN GetOneLineData
30 LET N = INT ( BRA ( 1 ) )
40 FOR I = 1 TO N
50 RUN GetOneLineData
60 LET { XX ( I ) = BRA ( 2 ) } :: { YY ( I ) = BRA ( 3 ) }
70 NEXT I
80 END
;
: main ( -- )
BASIC
10 run ResetAll
20 RUN FileData
30 RUN Once
40 END
;
page
.( Usage: ) cr
.( 1. test : for fixed input data ) cr
.( 2. hi : for interactive input data ) cr
.( 3. main : for data input from file ) cr
\S
Usage:
1. test : for fixed input data
2. hi : for interactive input data
3. main : for data input from file
ok
test
(1)Input data as following :
1 0.00000E0 0.00000E0
2 .157000 .156356
3 .314000 .308866
4 .471000 .453778
5 .628000 .587528
6 .785000 .706825
7 .942000 .808736
8 1.09900 .890753
9 1.25600 .950859
10 1.41300 .987576
11 1.55000 .999784
(2)Fitted equation is :
y = -.016422 + 1.18056 x -.332942 x^2
(3)Standard deviation of fit is : .011879
ok
main
Fileid is: 1452
387 Bytes are read into the Fadr RAM buffer!
(1)Input data as following :
1 0.00000E0 1.00000
2 .157000 .991811
3 .314000 .967564
4 .471000 .928206
5 .628000 .875260
6 .785000 .810746
7 .942000 .737075
8 1.09900 .656931
9 1.25600 .573130
10 1.41300 .488491
11 1.57000 .405696
(2)Fitted equation is :
y = 1.01339 -.121698 x -.175080 x^2
(3)Standard deviation of fit is : .011052
ok
第(44)個範例為矩陣因式分解的可用程式，合理的矩陣都可以因式分解成等效之上三角矩陣(upper triangular matrix)與下三角矩陣(lower triangular matrix)的乘積。
為了讓電腦的計算，能維持較高的精確度，在處理矩陣問題的領域，通常都會以一種挑出某列或某行中最大元素作為樞紐元素的方式進行前處理操作。樞紐元素會被轉置在例如對角線上的頂點位置，這個元素在計算過程中，通常要被作為分母來使用，太小或等於 0 者，就不適合作為分母。此時，為了讓最大元素被調整到樞紐位置而調整後的矩陣仍能等效於原矩陣，就產生了另一個也是在處理矩陣時經常要設計的程式技術，名稱就叫做轉換出樞紐矩陣(pivot matrix)。
矩陣元素的等效調動，就形同聯立方程式的前後調動，不影響聯立方程式的本質，所以等效。
因矩陣運算數學理論的探討已很完整，解矩陣問題時，通常也都喜歡採用在矩陣體系中直接進行矩陣運算來解決問題，而不是每次逐個處理單一元素的方式解決問題。矩陣被因式分解成三角矩陣後，有利於矩陣體系內的直覺觀察、解析與運算，許多計算也因此而可以簡化成只需處理對角線上的元素就夠了，所以才要搞矩陣的因式分解。前面的範例(17)求行列式之值，也用此法。
範例(44)有三個不同的程式，程式後面附有一個純文字的說明檔，其來有自。純文字檔記錄的是我從德州大學網站下載的 FORTRAN 原始程式，改寫後發現其中的陷阱，輸入數據不正確時，可能導致矩陣為線性相依的無解狀況。原程式中也有三處錯誤。這些問題，要在自己設計出程式、執行出結果後，才能除錯，我已在檔案中標示。我常從美國人的網頁下載各種技術資料，經常可以發現類似的陷阱，我覺得，貼網者是故意的。無能搞清楚問題者，就無能應用技術資料，而這些古聖先賢的創作，都很好用，我的 ABC Forth 系統能解決技術傳承時有人故意設置陷阱的問題，這裡的除錯展示，就是最佳的實例。
另兩個程式，則是從網上搜集而來，網頁中都有簡明的矩陣描述，這裡就不贅述。
(44)是單純的只作因式分解，為了簡明顯示效果，只用一個最簡單的三維矩陣為例，容易看得懂程式在幹甚麼。我曾提起過，我沒有進一步設計出矩陣體系，因為在搞矩陣研究之初，我設計的矩陣宣告格式中，沒有把矩陣的兩個維度儲存起來，所以，就算想印出矩陣時，還得先搬運到指定矩陣才能印出，不能碰到任何矩陣都可以直接印出它的所有內容。這項缺點，我在後來再設計系統時，已經改善，但我仍沒有機會繼續鑽研矩陣體系，總有一天，我會做。
(44-1)多增加了求解出樞紐矩陣的程式，也展示了此第(44)個程式可以直接擴充出(44-1)的效果，這就是 Forth 程式語言模組化程式的設計格式，具有自然可擴充的性質，想添加新功能，就儘管繼續往上增加。
(44-2)另添了特色，輸入數據改採使用『虛擬檔案』的方式送入程式。關於這項技術，將在第(45)個範例中提及，所以相關資訊，容留後敘，此乃我自創的技術，另闢專章來討論比較恰當。
:
為了讓電腦的計算，能維持較高的精確度，在處理矩陣問題的領域，通常都會以一種挑出某列或某行中最大元素作為樞紐元素的方式進行前處理操作。樞紐元素會被轉置在例如對角線上的頂點位置，這個元素在計算過程中，通常要被作為分母來使用，太小或等於 0 者，就不適合作為分母。此時，為了讓最大元素被調整到樞紐位置而調整後的矩陣仍能等效於原矩陣，就產生了另一個也是在處理矩陣時經常要設計的程式技術，名稱就叫做轉換出樞紐矩陣(pivot matrix)。
矩陣元素的等效調動，就形同聯立方程式的前後調動，不影響聯立方程式的本質，所以等效。
因矩陣運算數學理論的探討已很完整，解矩陣問題時，通常也都喜歡採用在矩陣體系中直接進行矩陣運算來解決問題，而不是每次逐個處理單一元素的方式解決問題。矩陣被因式分解成三角矩陣後，有利於矩陣體系內的直覺觀察、解析與運算，許多計算也因此而可以簡化成只需處理對角線上的元素就夠了，所以才要搞矩陣的因式分解。前面的範例(17)求行列式之值，也用此法。
範例(44)有三個不同的程式，程式後面附有一個純文字的說明檔，其來有自。純文字檔記錄的是我從德州大學網站下載的 FORTRAN 原始程式，改寫後發現其中的陷阱，輸入數據不正確時，可能導致矩陣為線性相依的無解狀況。原程式中也有三處錯誤。這些問題，要在自己設計出程式、執行出結果後，才能除錯，我已在檔案中標示。我常從美國人的網頁下載各種技術資料，經常可以發現類似的陷阱，我覺得，貼網者是故意的。無能搞清楚問題者，就無能應用技術資料，而這些古聖先賢的創作，都很好用，我的 ABC Forth 系統能解決技術傳承時有人故意設置陷阱的問題，這裡的除錯展示，就是最佳的實例。
另兩個程式，則是從網上搜集而來，網頁中都有簡明的矩陣描述，這裡就不贅述。
(44)是單純的只作因式分解，為了簡明顯示效果，只用一個最簡單的三維矩陣為例，容易看得懂程式在幹甚麼。我曾提起過，我沒有進一步設計出矩陣體系，因為在搞矩陣研究之初，我設計的矩陣宣告格式中，沒有把矩陣的兩個維度儲存起來，所以，就算想印出矩陣時，還得先搬運到指定矩陣才能印出，不能碰到任何矩陣都可以直接印出它的所有內容。這項缺點，我在後來再設計系統時，已經改善，但我仍沒有機會繼續鑽研矩陣體系，總有一天，我會做。
(44-1)多增加了求解出樞紐矩陣的程式，也展示了此第(44)個程式可以直接擴充出(44-1)的效果，這就是 Forth 程式語言模組化程式的設計格式，具有自然可擴充的性質，想添加新功能，就儘管繼續往上增加。
(44-2)另添了特色，輸入數據改採使用『虛擬檔案』的方式送入程式。關於這項技術，將在第(45)個範例中提及，所以相關資訊，容留後敘，此乃我自創的技術，另闢專章來討論比較恰當。
:
\ (44)LUFactorization.f
\ ABC FORTH demo code: Matrix LU factorization
\ Author: Ching-Tang Tseng, Hamilton, New Zealand
\ Date: 2015-03-27
\ Contact: ilikeforth@gmail.com
\ Reference:http://www.ma.utexas.edu/CNA/NA2/programs/f77code/CHP4/genlu.f
\ Matrix A is input.
\ L(Lower triangular matrix)and U{Upper triangular matrix)are output matrices.
\ W(Working matrix) used for checking.
\ MD is the unique setting in this program for Maximum matrix Dimension.
10 value MD
MD MD matrix a
MD MD matrix l
MD MD matrix u
\ Working matrix
MD MD matrix W
MD (array) t
4 Reals sum sum1 sum2 sum3
5 Integers i j k m n
((
: FixedData
[[ n = 3 ]]
{{ a ( 1 1 ) = 60 }} {{ a ( 1 2 ) = 30 }} {{ a ( 1 3 ) = 20 }}
{{ a ( 2 1 ) = 30 }} {{ a ( 2 2 ) = 20 }} {{ a ( 2 3 ) = 15 }}
{{ a ( 3 1 ) = 20 }} {{ a ( 3 2 ) = 15 }} {{ a ( 3 3 ) = 12 }}
;
: FixedData
[[ n = 3 ]]
{{ a ( 1 1 ) = 1 }} {{ a ( 1 2 ) = 3 }} {{ a ( 1 3 ) = 5 }}
{{ a ( 2 1 ) = 2 }} {{ a ( 2 2 ) = 4 }} {{ a ( 2 3 ) = 7 }}
{{ a ( 3 1 ) = 1 }} {{ a ( 3 2 ) = 1 }} {{ a ( 3 3 ) = 0 }}
;
))
: FixedData
[[ n = 3 ]]
{{ a ( 1 1 ) = 2 }} {{ a ( 1 2 ) = 4 }} {{ a ( 1 3 ) = 7 }}
{{ a ( 2 1 ) = 1 }} {{ a ( 2 2 ) = 3 }} {{ a ( 2 3 ) = 5 }}
{{ a ( 3 1 ) = 1 }} {{ a ( 3 2 ) = 1 }} {{ a ( 3 3 ) = 0 }}
;
: ResetAll
basic
10 for i = 1 to n
20 for j = 1 to n
30 let { a ( i j ) = 0 }
40 let { l ( i j ) = 0 }
50 let { u ( i i ) = 0 }
60 next j
70 next i
80 for i = 1 to n
90 let t ( i ) = 0
100 next i
110 end
;
: Initialize
basic
10 for i = 1 to n
20 let { l ( i i ) = 1 }
30 next i
40 let t ( 1 ) = 1
50 let t ( n ) = 1
60 end
;
: LU1
basic
10 for k = 1 to n
20 let { sum = a ( k k ) }
30 for m = 1 to k - 1
40 let { sum = sum - l ( k m ) * a ( m k ) } \ u --> a
50 next m
60 if t ( k ) = 1 then 90
70 let { l ( k k ) = sum / a ( k k ) } \ u --> a
80 goto 100
90 let { u ( k k ) = sum / a ( k k ) } \ l --> a
100 next k
110 end
;
: LU2
basic
10 for k = 1 to n
20 let { sum1 = a ( k k ) }
30 for m = 1 to k - 1
40 let { sum1 = sum1 - l ( k m ) * u ( m k ) }
50 next m
60 let { u ( k k ) = sum1 / l ( k k ) }
70 for j = k + 1 to n
80 let { sum2 = a ( k j ) }
90 for m = 1 to k - 1
100 let { sum2 = sum2 - l ( k m ) * u ( m j ) }
110 next m
120 let { u ( k j ) = sum2 / l ( k k ) }
130 next j
140 for i = k + 1 to n
150 let { sum3 = a ( i k ) }
160 for m = 1 to k - 1
170 let { sum3 = sum3 - l ( i m ) * u ( m k ) }
180 next m
190 let { l ( i k ) = sum3 / u ( k k ) }
200 next i
210 next k
220 end
;
: PrintMatrixW ( -- )
BASIC
10 RUN CR
20 FOR I = 1 TO N
30 FOR J = 1 TO N
40 LET { Sum = W ( I J ) } \ Based on W matrix
40 RUN Sum (F,) \ Right justified printing
50 NEXT J
60 RUN CR
70 NEXT I
80 END
;
: TotalAmount ( -- )
MD 1+ DUP * 3 CELLS * ;
: MoveAtoW ( -- )
{{ A ( 0 0 ) }} {{ W ( 0 0 ) }} TotalAmount CMOVE ;
: MoveLtoW ( -- )
{{ l ( 0 0 ) }} {{ W ( 0 0 ) }} TotalAmount CMOVE ;
: MoveUtoW ( -- )
{{ u ( 0 0 ) }} {{ W ( 0 0 ) }} TotalAmount CMOVE ;
: l*uToW ( -- )
BASIC
\ Square Matrix Multiplication
\ W(i,j) = l(i,j) * u(i,j)
10 FOR I = 1 TO N
20 FOR J = 1 TO N
30 LET { W ( I J ) = 0 }
40 FOR K = 1 TO N
50 LET { W ( I J ) = W ( I J ) + l ( I K ) * u ( K J ) }
60 NEXT K
70 NEXT J
80 NEXT I
90 END
;
: Checking ( -- )
l*uToW PrintMatrixW
;
: Print-alu
basic
10 run MoveAToW
20 print " a ( i j ) = "
30 run PrintMatrixW
40 run MoveLToW
50 print " l ( i j ) = "
60 run PrintMatrixW
70 run MoveUToW
80 print " u ( i j ) = "
90 run PrintMatrixW
100 end
;
: test
page ResetAll FixedData Initialize
cr ." *** Before *** " cr
Print-alu
LU1 LU2
cr ." *** After *** " cr
Print-alu
cr ." *** Checking *** " cr
cr ." W ( i j ) = l ( i j ) * u ( i j ) = a ( i j ) "
Checking
cr ." *** Finished *** " cr
;
test
\ ***** The end of program *****
\S
(1)Computer output:
*** Before ***
a ( i j ) =
60.0000 30.0000 20.0000
30.0000 20.0000 15.0000
20.0000 15.0000 12.0000
l ( i j ) =
1.00000 0.00000E0 0.00000E0
0.00000E0 1.00000 0.00000E0
0.00000E0 0.00000E0 1.00000
u ( i j ) =
0.00000E0 0.00000E0 0.00000E0
0.00000E0 0.00000E0 0.00000E0
0.00000E0 0.00000E0 0.00000E0
*** After ***
a ( i j ) =
60.0000 30.0000 20.0000
30.0000 20.0000 15.0000
20.0000 15.0000 12.0000
l ( i j ) =
1.00000 0.00000E0 0.00000E0
.500000 1.00000 0.00000E0
.333333 1.00000 1.00000
u ( i j ) =
60.0000 30.0000 20.0000
0.00000E0 5.00000 5.00000
0.00000E0 0.00000E0 .333333
*** Checking ***
W ( i j ) = l ( i j ) * u ( i j ) = a ( i j )
60.0000 30.0000 20.0000
30.0000 20.0000 15.0000
20.0000 15.0000 12.0000
*** Finished ***
ok
(2)References:(with Tags for bug fixes)
Website: http://www.ma.utexas.edu/CNA/NA2/programs/f77code/CHP4/genlu.f
c
c Numerical Analysis:
c Mathematics of Scientific Computing
c Third Edition
c D.R. Kincaid, E.W. Cheney
c Brooks/Cole Publ., 2002
c Copyright (c) 1996
c
c Section 4.2
c
c Example of general LU-factorization
c
c
c file: genlu.f
c
parameter (n=3)
dimension a(n,n),l(n,n),u(n,n)
real l
logical t(n)
data (a(1,j),j=1,n) / 60.0,30.0,20.0/
data (a(2,j),j=1,n) / 30.0,20.0,15.0/
data (a(3,j),j=1,n) / 20.0,15.0,12.0/
data l(1,1),u(2,2),l(3,3) /6.0,2.0,4.0/ --- meaningless data
data (t(i),i=1,n) /.true.,.false.,.true./ modified to /1.0,1.0,1.0/
c
print *
print *,' General LU-factorization example'
print *,' Section 4.2, Kincaid-Cheney'
print *
c
do 3 k=1,n
sum = a(k,k)
do 2 m=1,k-1
sum = sum - l(k,m)*u(m,k) ---Error, u(m,k) modified to a(m,k)
2 continue
if (t(k)) then
u(k,k) = sum/l(k,k) ---Error, l(k,k) modified to a(k,k)
else
l(k,k) = sum/u(k,k) ---Error, u(k,k) modified to a(k,k)
endif
3 continue
c
do 9 k=1,n
sum1 = a(k,k)
do 4 m=1,k-1
sum1 = sum1 - l(k,m)*u(m,k)
4 continue
u(k,k) = sum1/l(k,k)
do 6 j=k+1,n
sum2 = a(k,j)
do 5 m=1,k-1
sum2 = sum2 - l(k,m)*u(m,j)
5 continue
u(k,j) = sum2/l(k,k)
6 continue
do 8 i=k+1,n
sum3 = a(i,k)
do 7 m=1,k-1
sum3 = sum3 - l(i,m)*u(m,k)
7 continue
l(i,k) = sum3/u(k,k)
8 continue
9 continue
c
do 11 i=1,n
do 10 j=1,n
print 12,i,j,l(i,j),i,j,u(i,j)
10 continue
11 continue
c
12 format (1x,'l(',i2,',',i2,') =',e13.6,7x,
+ 'u(',i2,',',i2,') =',e13.6)
stop
end
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\ (44-1)PLUFactorization.f
\ ABC FORTH demo code: Matrix PLU factorization
\ Author: Ching-Tang Tseng, Hamilton, New Zealand
\ Date: 2015-03-30
\ Contact: ilikeforth@gmail.com
\ Reference:
\ (1)http://www.ma.utexas.edu/CNA/NA2/programs/f77code/CHP4/genlu.f
\ (2)Rosettacode website as well.
\ Description:
\ Matrix a is an input raw data matrix.
\ (1)l(Lower triangular matrix) and
\ (2)u{Upper triangular matrix) and
\ (3)p(Pivot matrix) are output matrices.
\ UU, VV, WW(Working matrix) are used for checking.
\ MD is the unique setting in this program for Maximum matrix Dimension.
10 value MD
MD MD matrix a
MD MD matrix l
MD MD matrix u
MD MD matrix p
\ All raw data of input matrix a is stored in matrix q for the checking at last.
MD MD matrix q
\ Working matrices for multiplication using
MD MD matrix UU
MD MD matrix VV
MD MD matrix WW
MD (array) t
4 Reals sum sum1 sum2 sum3
5 Integers i j k m n
((
: FixedData
[[ n = 3 ]]
{{ a ( 1 1 ) = 60 }} {{ a ( 1 2 ) = 30 }} {{ a ( 1 3 ) = 20 }}
{{ a ( 2 1 ) = 30 }} {{ a ( 2 2 ) = 20 }} {{ a ( 2 3 ) = 15 }}
{{ a ( 3 1 ) = 20 }} {{ a ( 3 2 ) = 15 }} {{ a ( 3 3 ) = 12 }}
;
))
: FixedData
[[ n = 4 ]]
{{ a ( 1 1 ) = 11 }} {{ a ( 1 2 ) = 9 }} {{ a ( 1 3 ) = 24 }} {{ a ( 1 4 ) = 2 }}
{{ a ( 2 1 ) = 1 }} {{ a ( 2 2 ) = 5 }} {{ a ( 2 3 ) = 2 }} {{ a ( 2 4 ) = 6 }}
{{ a ( 3 1 ) = 3 }} {{ a ( 3 2 ) = 17 }} {{ a ( 3 3 ) = 18 }} {{ a ( 3 4 ) = 1 }}
{{ a ( 4 1 ) = 2 }} {{ a ( 4 2 ) = 5 }} {{ a ( 4 3 ) = 7 }} {{ a ( 4 4 ) = 1 }}
;
: ResetAll
basic
10 for i = 1 to n
20 for j = 1 to n
30 let { a ( i j ) = 0 }
40 let { l ( i j ) = 0 }
50 let { u ( i j ) = 0 }
60 let { p ( i j ) = 0 }
70 next j
80 next i
90 for i = 1 to n
100 let t ( i ) = 0
110 next i
120 end
;
\ l and p are unit matrices at beginning.
: Initialize
basic
10 for i = 1 to n
20 let { l ( i i ) = 1 }
30 let { p ( i i ) = 1 }
40 next i
50 let t ( 1 ) = 1
60 let t ( n ) = 1
70 end
;
\ Pivot computes out a pivot matrix
integer r
: Pivot
basic
10 for i = 1 to n
20 let { sum = a ( i i ) }
30 let r = i
40 for j = i to n
50 if { a ( j i ) > sum } then 70
60 goto 90
70 let { sum = a ( j i ) }
80 let r = j
90 next j
100 if i =/= r then 120
110 goto 170
120 for j = 1 to n
130 let { sum1 = p ( i j ) }
140 let { p ( i j ) = p ( r j ) }
150 let { p ( r j ) = sum1 }
160 next j
170 next i
180 end
;
: LU1
basic
10 for k = 1 to n
20 let { sum = a ( k k ) }
30 for m = 1 to k - 1
40 let { sum = sum - l ( k m ) * a ( m k ) } \ u(m,k) --> a(m,k)
50 next m
60 if t ( k ) = 1 then 90
70 let { l ( k k ) = sum / a ( k k ) } \ u(k,k) --> a(k,k)
80 goto 100
90 let { u ( k k ) = sum / a ( k k ) } \ l(k,k) --> a(k,k)
100 next k
110 end
;
: LU2
basic
10 for k = 1 to n
20 let { sum1 = a ( k k ) }
30 for m = 1 to k - 1
40 let { sum1 = sum1 - l ( k m ) * u ( m k ) }
50 next m
60 let { u ( k k ) = sum1 / l ( k k ) }
70 for j = k + 1 to n
80 let { sum2 = a ( k j ) }
90 for m = 1 to k - 1
100 let { sum2 = sum2 - l ( k m ) * u ( m j ) }
110 next m
120 let { u ( k j ) = sum2 / l ( k k ) }
130 next j
140 for i = k + 1 to n
150 let { sum3 = a ( i k ) }
160 for m = 1 to k - 1
170 let { sum3 = sum3 - l ( i m ) * u ( m k ) }
180 next m
190 let { l ( i k ) = sum3 / u ( k k ) }
200 next i
210 next k
220 end
;
\ All ready to print out matrices are based on WW matrix
: PrintWW ( -- )
BASIC
10 RUN CR
20 FOR I = 1 TO N
30 FOR J = 1 TO N
40 LET { Sum = WW ( I J ) }
40 RUN Sum (F,) \ Right justified printing
50 NEXT J
60 RUN CR
70 NEXT I
80 END
;
: TotalAmount ( -- )
MD 1+ DUP * 3 CELLS * ;
: MatrixCopy ( addr addr -- )
TotalAmount CMOVE ;
: CopyatoWW ( -- )
{{ a ( 0 0 ) }} {{ WW ( 0 0 ) }} MatrixCopy ;
: CopyltoWW ( -- )
{{ l ( 0 0 ) }} {{ WW ( 0 0 ) }} MatrixCopy ;
: CopyutoWW ( -- )
{{ u ( 0 0 ) }} {{ WW ( 0 0 ) }} MatrixCopy ;
: CopypToWW
{{ p ( 0 0 ) }} {{ WW ( 0 0 ) }} MatrixCopy ;
: CopyqToWW
{{ q ( 0 0 ) }} {{ WW ( 0 0 ) }} MatrixCopy ;
: CopyaToQ
{{ a ( 0 0 ) }} {{ q ( 0 0 ) }} MatrixCopy ;
: CopypToUU
{{ p ( 0 0 ) }} {{ UU ( 0 0 ) }} MatrixCopy ;
: CopyaToVV
{{ A ( 0 0 ) }} {{ VV ( 0 0 ) }} MatrixCopy ;
: CopyWWToA
{{ WW ( 0 0 ) }} {{ A ( 0 0 ) }} MatrixCopy ;
: CopylToUU
{{ l ( 0 0 ) }} {{ UU ( 0 0 ) }} MatrixCopy ;
: CopyuToVV
{{ u ( 0 0 ) }} {{ VV ( 0 0 ) }} MatrixCopy ;
\ Matrix Multiplication
\ WW(i,j) = UU(i,j) * VV(i,j)
\ WW = UU * VV
: Matrix* ( -- )
BASIC
10 FOR I = 1 TO N
20 FOR J = 1 TO N
30 LET { WW ( I J ) = 0 }
40 FOR K = 1 TO N
50 LET { WW ( I J ) = WW ( I J ) + UU ( I K ) * VV ( K J ) }
60 NEXT K
70 NEXT J
80 NEXT I
90 END
;
: Print-alup
basic
10 run CopyAToWW
20 print " a ( i j ) = "
30 run PrintWW
40 run CopyLToWW
50 print " l ( i j ) = "
60 run PrintWW
70 run CopyUToWW
80 print " u ( i j ) = "
90 run PrintWW
\ Both of the following code formats are the same meaning.
\ 100 run CopyptoWW
100 run
address-of p ( 0 0 ) fmatrix
address-of WW ( 0 0 ) fmatrix
MatrixCopy
110 print " p ( i j ) = "
120 run PrintWW
130 end
;
: test
page ResetAll FixedData Initialize
cr ." *** Before *** " cr
Print-alup
Copyatoq
pivot CopyptoUU CopyatoVV matrix* CopyWWtoa
cr ." *** After pivot *** " cr
Print-alup
LU1 LU2
cr ." *** After Factorization *** " cr
Print-alup
cr ." *** Checking *** " cr
cr ." a ( i j ) = "
CopyqtoWW PrintWW
cr ." p ( i j ) * a ( i j ) = "
CopyatoWW PrintWW
cr ." l ( i j ) * u ( i j ) = "
CopyltoUU CopyutoVV Matrix* PrintWW
cr ." *** Finished *** " cr
;
test
\ ***** The end of program *****
\S
References: Rosettacode website
DEF PROCpivot(a(), RETURN p())
LOCAL i%, j%, m%, n%, r% : n% = DIM(a(),2)
DIM p(n%,n%) : FOR i% = 0 TO n% : p(i%,i%) = 1 : NEXT
FOR i% = 0 TO n%
m% = a(i%,i%)
r% = i%
FOR j% = i% TO n%
IF a(j%,i%) > m% m% = a(j%,i%) : r% = j%
NEXT
IF i%<>r% THEN
FOR j% = 0 TO n% : SWAP p(i%,j%),p(r%,j%) : NEXT
ENDIF
NEXT i%
ENDPROC
\\\\\\\\\\\\\\\\\\\\\\\\
\ (44-1)PLUFileIO.f
\ ABC FORTH demo code: Matrix PLU factorization
\ Author: Ching-Tang Tseng, Hamilton, New Zealand
\ Date: 2015-03-30
\ Contact: ilikeforth@gmail.com
\ Reference:
\ (1)http://www.ma.utexas.edu/CNA/NA2/programs/f77code/CHP4/genlu.f
\ (2)Rosettacode website as well.
\ Description:
\ Matrix a is an input raw data matrix.
\ (1)l(Lower triangular matrix) and
\ (2)u{Upper triangular matrix) and
\ (3)p(Pivot matrix) are output matrices.
\ UU, VV, WW(Working matrix) are used for checking.
\ MD is the unique setting in this program for Maximum matrix Dimension.
10 value MD
MD MD matrix a
MD MD matrix l
MD MD matrix u
MD MD matrix p
\ All raw data of input matrix a is stored in matrix q for the checking at last.
MD MD matrix q
\ Working matrices, for multiplication using
MD MD matrix UU
MD MD matrix VV
MD MD matrix WW
MD (array) t
4 Reals sum sum1 sum2 sum3
6 Integers i j k m n j+1
: MakeAPseudoFile
S" 4" firstline
\ Leading number marks the line order.
S" 1 11 9 24 2" nextline
S" 2 1 5 2 6" nextline
S" 3 3 17 18 1" nextline
S" 4 2 5 7 1" nextline
S>FileBuffer ;
: ReadData ( -- )
BASIC
10 RUN MakeAPseudoFile
20 RUN GetOneLineData
30 LET N = INT ( BRA ( 1 ) )
40 FOR i = 1 TO N
50 RUN GetOneLineData
60 FOR j = 1 TO N
70 let j+1 = j + 1
\ Ignore all BRA(1)
80 LET { A ( i j ) = BRA ( j+1 ) }
90 NEXT j
100 NEXT i
110 END
;
((
: FixedData
[[ n = 3 ]]
{{ a ( 1 1 ) = 60 }} {{ a ( 1 2 ) = 30 }} {{ a ( 1 3 ) = 20 }}
{{ a ( 2 1 ) = 30 }} {{ a ( 2 2 ) = 20 }} {{ a ( 2 3 ) = 15 }}
{{ a ( 3 1 ) = 20 }} {{ a ( 3 2 ) = 15 }} {{ a ( 3 3 ) = 12 }}
;
: FixedData
[[ n = 4 ]]
{{ a ( 1 1 ) = 11 }} {{ a ( 1 2 ) = 9 }} {{ a ( 1 3 ) = 24 }} {{ a ( 1 4 ) = 2 }}
{{ a ( 2 1 ) = 1 }} {{ a ( 2 2 ) = 5 }} {{ a ( 2 3 ) = 2 }} {{ a ( 2 4 ) = 6 }}
{{ a ( 3 1 ) = 3 }} {{ a ( 3 2 ) = 17 }} {{ a ( 3 3 ) = 18 }} {{ a ( 3 4 ) = 1 }}
{{ a ( 4 1 ) = 2 }} {{ a ( 4 2 ) = 5 }} {{ a ( 4 3 ) = 7 }} {{ a ( 4 4 ) = 1 }}
;
))
: ResetAll
basic
10 for i = 1 to n
20 for j = 1 to n
30 let { a ( i j ) = 0 }
40 let { l ( i j ) = 0 }
50 let { u ( i j ) = 0 }
60 let { p ( i j ) = 0 }
70 next j
80 next i
90 for i = 1 to n
100 let t ( i ) = 0
110 next i
120 end
;
\ l and p are unit matrices at beginning.
: Initialize
basic
10 for i = 1 to n
20 let { l ( i i ) = 1 }
30 let { p ( i i ) = 1 }
40 next i
50 let t ( 1 ) = 1
60 let t ( n ) = 1
70 end
;
\ Pivot computes out a pivot matrix
integer r
: Pivot
basic
10 for i = 1 to n
20 let { sum = a ( i i ) }
30 let r = i
40 for j = i to n
50 if { a ( j i ) > sum } then 70
60 goto 90
70 let { sum = a ( j i ) }
80 let r = j
90 next j
100 if i =/= r then 120
110 goto 170
120 for j = 1 to n
130 let { sum1 = p ( i j ) }
140 let { p ( i j ) = p ( r j ) }
150 let { p ( r j ) = sum1 }
160 next j
170 next i
180 end
;
: LU1
basic
10 for k = 1 to n
20 let { sum = a ( k k ) }
30 for m = 1 to k - 1
40 let { sum = sum - l ( k m ) * a ( m k ) } \ u(m,k) --> a(m,k)
50 next m
60 if t ( k ) = 1 then 90
70 let { l ( k k ) = sum / a ( k k ) } \ u(k,k) --> a(k,k)
80 goto 100
90 let { u ( k k ) = sum / a ( k k ) } \ l(k,k) --> a(k,k)
100 next k
110 end
;
: LU2
basic
10 for k = 1 to n
20 let { sum1 = a ( k k ) }
30 for m = 1 to k - 1
40 let { sum1 = sum1 - l ( k m ) * u ( m k ) }
50 next m
60 let { u ( k k ) = sum1 / l ( k k ) }
70 for j = k + 1 to n
80 let { sum2 = a ( k j ) }
90 for m = 1 to k - 1
100 let { sum2 = sum2 - l ( k m ) * u ( m j ) }
110 next m
120 let { u ( k j ) = sum2 / l ( k k ) }
130 next j
140 for i = k + 1 to n
150 let { sum3 = a ( i k ) }
160 for m = 1 to k - 1
170 let { sum3 = sum3 - l ( i m ) * u ( m k ) }
180 next m
190 let { l ( i k ) = sum3 / u ( k k ) }
200 next i
210 next k
220 end
;
\ All ready to print out matrices are based on WW matrix
: PrintWW ( -- )
BASIC
10 RUN CR
20 FOR I = 1 TO N
30 FOR J = 1 TO N
40 LET { Sum = WW ( I J ) }
40 RUN Sum (F,) \ Right justified printing
50 NEXT J
60 RUN CR
70 NEXT I
80 END
;
: TotalAmount ( -- )
MD 1+ DUP * 3 CELLS * ;
: MatrixCopy ( addr addr -- )
TotalAmount CMOVE ;
: CopyatoWW ( -- )
{{ a ( 0 0 ) }} {{ WW ( 0 0 ) }} MatrixCopy ;
: CopyltoWW ( -- )
{{ l ( 0 0 ) }} {{ WW ( 0 0 ) }} MatrixCopy ;
: CopyutoWW ( -- )
{{ u ( 0 0 ) }} {{ WW ( 0 0 ) }} MatrixCopy ;
: CopypToWW
{{ p ( 0 0 ) }} {{ WW ( 0 0 ) }} MatrixCopy ;
: CopyqToWW
{{ q ( 0 0 ) }} {{ WW ( 0 0 ) }} MatrixCopy ;
: CopyaToQ
{{ a ( 0 0 ) }} {{ q ( 0 0 ) }} MatrixCopy ;
: CopypToUU
{{ p ( 0 0 ) }} {{ UU ( 0 0 ) }} MatrixCopy ;
: CopyaToVV
{{ A ( 0 0 ) }} {{ VV ( 0 0 ) }} MatrixCopy ;
: CopyWWToA
{{ WW ( 0 0 ) }} {{ A ( 0 0 ) }} MatrixCopy ;
: CopylToUU
{{ l ( 0 0 ) }} {{ UU ( 0 0 ) }} MatrixCopy ;
: CopyuToVV
{{ u ( 0 0 ) }} {{ VV ( 0 0 ) }} MatrixCopy ;
\ Matrix Multiplication
\ WW(i,j) = UU(i,j) * VV(i,j)
\ WW = UU * VV
: Matrix* ( -- )
BASIC
10 FOR I = 1 TO N
20 FOR J = 1 TO N
30 LET { WW ( I J ) = 0 }
40 FOR K = 1 TO N
50 LET { WW ( I J ) = WW ( I J ) + UU ( I K ) * VV ( K J ) }
60 NEXT K
70 NEXT J
80 NEXT I
90 END
;
: Print-alup
basic
10 run CopyAToWW
20 print " a ( i j ) = "
30 run PrintWW
40 run CopyLToWW
50 print " l ( i j ) = "
60 run PrintWW
70 run CopyUToWW
80 print " u ( i j ) = "
90 run PrintWW
\ Both of the following formats are the same meaning.
\ 100 run CopyptoWW
100 run
address-of p ( 0 0 ) fmatrix
address-of WW ( 0 0 ) fmatrix
MatrixCopy
110 print " p ( i j ) = "
120 run PrintWW
130 end
;
: test
\ page ResetAll FixedData Initialize
page ResetAll ReadData Initialize
cr ." *** Before *** " cr
Print-alup
Copyatoq
pivot CopyptoUU CopyatoVV matrix* CopyWWtoa
cr ." *** After pivot *** " cr
Print-alup
LU1 LU2
cr ." *** After Factorization *** " cr
Print-alup
cr ." *** Checking *** " cr
cr ." a ( i j ) = "
CopyqtoWW PrintWW
cr ." p ( i j ) * a ( i j ) = "
CopyatoWW PrintWW
cr ." l ( i j ) * u ( i j ) = "
CopyltoUU CopyutoVV Matrix* PrintWW
cr ." *** Finished *** " cr
;
test
\ ***** The end of program *****
\S
References: Rosettacode website
DEF PROCpivot(a(), RETURN p())
LOCAL i%, j%, m%, n%, r% : n% = DIM(a(),2)
DIM p(n%,n%) : FOR i% = 0 TO n% : p(i%,i%) = 1 : NEXT
FOR i% = 0 TO n%
m% = a(i%,i%)
r% = i%
FOR j% = i% TO n%
IF a(j%,i%) > m% m% = a(j%,i%) : r% = j%
NEXT
IF i% <> r% THEN
FOR j% = 0 TO n% : SWAP p(i%,j%),p(r%,j%) : NEXT
ENDIF
NEXT i%
ENDPROC
第(45)個範例介紹 Forth 可以在 Windows’ 作業系統中開啟一種虛擬檔案來使用的特殊應用程式。
此前，全世界沒有人肯在網上公開的場合，介紹 Forth 操控檔案之事。為什麼？因為事情太複雜了，還高度受制於作業系統，換個環境，辛苦設計的東西，就都會不再能使用。微軟的 OS 又不是最好的東西，資料封閉又難用，誰願意搞？除非有人願意付錢，否則免談。不只 Forth 如此，別種程式語言也差不多，我們可從公開場合獲得的資料相當貧乏，還有其他許多原因，令大家對操作檔案之事，越來越陌生了。
換個方式來思考這個問題，能不能不搞檔案？當然可以不搞。只不過，不搞之後，很多事情就都得人工操作才能完成。由於現今的電腦，很容易操作滑鼠來反白選定、點選複製、轉檔貼上，所以大家才懶得再搞檔案之事，也因此才能不搞。
但是，如果想把事情全面自動化，不搞檔案就不行。我們的 Forth 系統，通常都具有一鍵通(turnkey)的功能，也就是指定執行檔案後，所有的事情，一按鍵就能完成。沒有輸入資料時，事情還很好辦，輸入資料量很大時，不能操作檔案就難辦了。前面有好幾個範例，已經介紹過讓數據從檔案輸入的技術，雖都只是小量的數據，但卻代表著也能搞千萬筆資料的意義，只是我沒必要這樣展示而已。另外，如果資料是從別種套裝軟體貢獻過來的，例如： Excel 的數據， Word 的文字， email 傳輸過來的資料、與 Python 一起跑程式 ..... 等等事情，不搞檔案就是不行，或根本不方便。這些事情我全做過，所以，我自己搞出了一套操作檔案的辦法。
我在不同的階段與不同的作業系統中搞過檔案操作，搞到後來，才發現，只有在 Windows’ OS 中，可以搞得了一種形同虛擬檔案的產品，在其他的作業系統中都不行。但這也不是什麼大不了的問題，不行，那就真正的產生實體檔案又何妨？大不了，用後不想要了，就殺掉，效果也是一樣的。
所謂的虛擬，就是用時感覺它存在，用後它就自然消失了，我是怎麼讓它虛擬化的？沒什麼，因為 Forth 很容易操控出將資料存取於記憶體的作業，我在建立系統之初，就規劃出了能安放檔案的記憶體區塊，那麼，就直接對該區域存取數據，不就得了。只是用時，不執行存檔(save-file)指令而已。
如果您仔細想過，我們在寫程式時，所需要之檔案操作指令，實際上三個就夠了，它們是：
Get-File ( address length — )
New-File ( address length — )
Save-File ( address length — )
為什麼？我留給大家自己思考。至於 Forth 系統中提供了許多與 OS 有關的檔案操作指令，大部分都可以不用，例如：複製(copy)、貼上(post)、或刪除(delete)。但上述三個絕對需要的指令，系統卻不提供，只提供了還有後續問題的開檔(open-file)、讀檔(read-file)、寫出檔案(write-file)或關檔(close-file)，我只好自己設計。您如果不用我設計的這三個指令，並照我的規矩來搞檔案，可能得痛苦很久，都搞不成檔案操作之事。
規劃檔案操作時，我不喜歡後來許多人都把檔案開在 Forth 系統外面的方式搞檔案。所以，我都把落實的 Forth 系統開得很大，檔案使用的記憶體區塊就被指定在 Forth 系統中，所有的指令就都可以對該處記憶體進行自由的存取操作。我規劃出來的幾個記憶體指標，列示如下：
Fptr \ file pointer
Frem \ file remainder
Flen \ file length
另外有3個相關變數：
Fadr \ file starting address
FileID \ file identification
Fsize \ file size
知道這些東西，就能操作檔案了。事實上，我設計的系統也允許把檔案開在 Forth 系統之外，只需把檔案的起始位址指定到 Fadr 變數去，就可以了。我不喜歡那樣搞檔案，我只喜歡把所有的東西都放在 Forth 系統內來搞。
:
此前，全世界沒有人肯在網上公開的場合，介紹 Forth 操控檔案之事。為什麼？因為事情太複雜了，還高度受制於作業系統，換個環境，辛苦設計的東西，就都會不再能使用。微軟的 OS 又不是最好的東西，資料封閉又難用，誰願意搞？除非有人願意付錢，否則免談。不只 Forth 如此，別種程式語言也差不多，我們可從公開場合獲得的資料相當貧乏，還有其他許多原因，令大家對操作檔案之事，越來越陌生了。
換個方式來思考這個問題，能不能不搞檔案？當然可以不搞。只不過，不搞之後，很多事情就都得人工操作才能完成。由於現今的電腦，很容易操作滑鼠來反白選定、點選複製、轉檔貼上，所以大家才懶得再搞檔案之事，也因此才能不搞。
但是，如果想把事情全面自動化，不搞檔案就不行。我們的 Forth 系統，通常都具有一鍵通(turnkey)的功能，也就是指定執行檔案後，所有的事情，一按鍵就能完成。沒有輸入資料時，事情還很好辦，輸入資料量很大時，不能操作檔案就難辦了。前面有好幾個範例，已經介紹過讓數據從檔案輸入的技術，雖都只是小量的數據，但卻代表著也能搞千萬筆資料的意義，只是我沒必要這樣展示而已。另外，如果資料是從別種套裝軟體貢獻過來的，例如： Excel 的數據， Word 的文字， email 傳輸過來的資料、與 Python 一起跑程式 ..... 等等事情，不搞檔案就是不行，或根本不方便。這些事情我全做過，所以，我自己搞出了一套操作檔案的辦法。
我在不同的階段與不同的作業系統中搞過檔案操作，搞到後來，才發現，只有在 Windows’ OS 中，可以搞得了一種形同虛擬檔案的產品，在其他的作業系統中都不行。但這也不是什麼大不了的問題，不行，那就真正的產生實體檔案又何妨？大不了，用後不想要了，就殺掉，效果也是一樣的。
所謂的虛擬，就是用時感覺它存在，用後它就自然消失了，我是怎麼讓它虛擬化的？沒什麼，因為 Forth 很容易操控出將資料存取於記憶體的作業，我在建立系統之初，就規劃出了能安放檔案的記憶體區塊，那麼，就直接對該區域存取數據，不就得了。只是用時，不執行存檔(save-file)指令而已。
如果您仔細想過，我們在寫程式時，所需要之檔案操作指令，實際上三個就夠了，它們是：
Get-File ( address length — )
New-File ( address length — )
Save-File ( address length — )
為什麼？我留給大家自己思考。至於 Forth 系統中提供了許多與 OS 有關的檔案操作指令，大部分都可以不用，例如：複製(copy)、貼上(post)、或刪除(delete)。但上述三個絕對需要的指令，系統卻不提供，只提供了還有後續問題的開檔(open-file)、讀檔(read-file)、寫出檔案(write-file)或關檔(close-file)，我只好自己設計。您如果不用我設計的這三個指令，並照我的規矩來搞檔案，可能得痛苦很久，都搞不成檔案操作之事。
規劃檔案操作時，我不喜歡後來許多人都把檔案開在 Forth 系統外面的方式搞檔案。所以，我都把落實的 Forth 系統開得很大，檔案使用的記憶體區塊就被指定在 Forth 系統中，所有的指令就都可以對該處記憶體進行自由的存取操作。我規劃出來的幾個記憶體指標，列示如下：
Fptr \ file pointer
Frem \ file remainder
Flen \ file length
另外有3個相關變數：
Fadr \ file starting address
FileID \ file identification
Fsize \ file size
知道這些東西，就能操作檔案了。事實上，我設計的系統也允許把檔案開在 Forth 系統之外，只需把檔案的起始位址指定到 Fadr 變數去，就可以了。我不喜歡那樣搞檔案，我只喜歡把所有的東西都放在 Forth 系統內來搞。
:
\ (45)MakeAPseudoFile.f
\ ABC FORTH demo code: Make a pseudo file to simplify data input code.
\ Author: Ching-Tang Tseng, Hamilton, New Zealand
\ Date: 2015-03-30
\ Contact: ilikeforth@gmail.com
: FirstLine ( adr len -- )
fadr Fsize 0 FILL
Fadr to Fptr
DUP +to Fptr
Fadr swap cmove
13 Fptr c! 1 +to Fptr
10 Fptr c! 1 +to Fptr
;
: NextLine ( adr len -- )
DUP >R
Fptr SWAP CMOVE
R> +to Fptr
13 Fptr c! 1 +to Fptr
10 Fptr c! 1 +to Fptr
;
: S>FileBuffer ( -- )
26 Fptr c! 1 +to Fptr
0 Fptr c! 1 +to Fptr
Fptr Fadr - to Flen
Fadr to Fptr Flen to frem
;
\ 此套功能建進 V.659 系統永久使用，增加上列三個指令。
\S
Usage:
fload (45)MakeAPseudoFile.f
: MakeAPseudoFile
S" 4" firstline
\ Leading number marks the line order
S" 1 11 9 24 2" nextline
S" 2 1 5 2 6" nextline
S" 3 3 17 18 1" nextline
S" 4 2 5 7 1" nextline
S>FileBuffer ;
: ReadData ( -- )
BASIC
10 RUN MakeAPseudoFile
20 RUN GetOneLineData
30 LET N = INT ( BRA ( 1 ) )
40 FOR i = 1 TO N
50 RUN GetOneLineData
60 FOR j = 1 TO N
70 let j+1 = j + 1
\ Ignore all BRA(1)
80 LET { A ( i j ) = BRA ( j+1 ) }
90 NEXT j
100 NEXT i
110 END
;
: check
MakeAPseudoFile
cr
Fadr Flen type
cr
;
訂閱：
張貼留言 (Atom)
沒有留言:
張貼留言