一百個例題 (91 ~95)
Ching-Tang Tseng
Hamilton, New Zealand
6 October 2024
FORTH:慶的世界。 本網頁為中華民國(Republic Of China)FORTH語言協會 永久性網頁
2024年10月5日 星期六
第(91)個範例是專門探討以二分法(bisection)求解方程式實數根的技術。
請注意前面的第(76)個範例，也是利用二分法技術，用來算出浮點數的平方根。
這個方法很好用，也絕對可靠，只要對此技術認識得夠清楚，就能用來解決所有方程式求根的問題。
我個人認為，二分法是當今全世界從事於電腦數學運算應用者都應該熟悉的技術。
這個範例的全套英文解說，曾經貼示於我的個人網頁，是 20180416 所刊出的文章： Bisection Method To Solve Equation 。（該文因網頁改版，原貼內容已遭破壞，我會找時間進行修改。）
範例所解的題目，是經過刻意挑選的，您無法手解出它所有的三個實根。題目中的方程式，是代表某種化學反應趨勢的實測數據描述出來的表示式。當初，我大女兒念生化系時，拿這個方程式回家來問我，如何才能把正的兩個實根給算出來？那兩個根，要用來定出激烈化學反應的界限。於是，我發展出了這套程式。
我用英文貼文，大部分台灣人都看不懂，教大家看英文資料，不是我的責任，我貼文的目的，是讓全世界知道我的系統能解這種問題。很多人看得懂英文，但也只看出了表面意義。唯一看過此文，並在國際論壇回應網貼意見者，只有 ciForth 的作者，他才意識到系統能解此題的一個很重要的特點，他貼文表示此特點必須為後世所師法。
我若不刻意指出這個特點，大家可能也很難注意到，但如果從函數圖來看問題，就能看出問題所在了。函數圖中在 x = 80 的位置，是一個奇點，它的函數值是正的無限大，過去，所有的數學計算系統，在演算過程中，若碰到這種計算，會立刻執行例外處理、發出警告、停止執行。可是，我的系統就故意從這一個奇點開始計算，還進行了兩次，系統碰到了，照樣繼續算，直到產生結果，不需要發出任何警告。我的系統可以接受無限大的數字，而且，這個數字恆就是個無限大，使用者若想印出來看，系統就印出是帶有 + 或 - 符號的 infinit 。
二分法的基本精神，前文已經介紹過，這個範例則專挑學校老師也無法手解的問題來示範二分法的功能，它的發揮，主要是必須先靠繪出函數圖來輔助，您才能確認函數的實根大約會落在那裡？現行電腦的功能，已經能夠非常快速輕鬆的解決函數繪圖問題，求解實根時所需提供給程式使用的上下限二值，根據圖形就很容易取得，剩下來的冗長計算，全部交給電腦。
我在發展到 Lina64 系統之後，有感於這種可用的寶貴程式值得庫存。因此，開始把程式寫得非常整齊，並添加了許多必要的註解，讓自己在日後想要再度使用時，一目了然：
新給的方程式定義，就只需修改 FX 的定義內容，而且，寫法就是日常書寫的格式。
新給的上下限二個值，就只需修改 root1 , root2 , root3 ..... 等等就夠了。
仍然只需執行 main 指令，就能得到所有的答案。
我在 Lina64 中也發展了語音輸出及單音樂譜的執行能力，所以在貼文時，故意展示出這兩項功能。
語音播出是『莫斯科廣播電台的台號音樂與呼叫』，人聲與樂聲均能使用程式輸出。
繪圖的輸出，是藉 Forth 系統控制 Python 程式去操控 GnuＰlot ，是個論秒就能繪出函數圖的設計，那個環境，可以直接操作圖形歸檔。
所有執行程序均在 Lina64 中一氣呵成。
浮點運算部分與 ABC Forth 的功能，是我個人的創作，此系統未曾釋出，也不想釋出，因為大家可以注意看所有的浮點數輸入方式，跟此前全世界的任何系統都不一樣，假數部分與指數部分是被中間插入一個 e 或 E 來完成的，執行速度很快，因為數字都不需要轉換，就能直接算出浮點運算的結果。我不想擾亂世界的規矩，自己使用就無所謂，所以不釋出系統。
系統的所有功能， Lina64 的原作者大約也很難完成，所以他仔細讀過我的貼文後刻意回應。他說，他就是希望他的系統能夠辦到這些事情。他講的一點也沒錯，能實現這些功能的基本要素，是他的系統內已經設計好了跟作業系統所需要的介面工具指令，有好幾個，各個都能正確執行出對的結果。只是這位作者可能沒有像我一樣去普讀許多系統資料，供自己套進 Lina64 使用。這位作者的出身是高能物理的博士，我跟他一樣，都懂核子科學。我未曾與他聯繫過，只在國際論壇上神交，但我很感激他公開所有設計的系統。
:
請注意前面的第(76)個範例，也是利用二分法技術，用來算出浮點數的平方根。
這個方法很好用，也絕對可靠，只要對此技術認識得夠清楚，就能用來解決所有方程式求根的問題。
我個人認為，二分法是當今全世界從事於電腦數學運算應用者都應該熟悉的技術。
這個範例的全套英文解說，曾經貼示於我的個人網頁，是 20180416 所刊出的文章： Bisection Method To Solve Equation 。（該文因網頁改版，原貼內容已遭破壞，我會找時間進行修改。）
範例所解的題目，是經過刻意挑選的，您無法手解出它所有的三個實根。題目中的方程式，是代表某種化學反應趨勢的實測數據描述出來的表示式。當初，我大女兒念生化系時，拿這個方程式回家來問我，如何才能把正的兩個實根給算出來？那兩個根，要用來定出激烈化學反應的界限。於是，我發展出了這套程式。
我用英文貼文，大部分台灣人都看不懂，教大家看英文資料，不是我的責任，我貼文的目的，是讓全世界知道我的系統能解這種問題。很多人看得懂英文，但也只看出了表面意義。唯一看過此文，並在國際論壇回應網貼意見者，只有 ciForth 的作者，他才意識到系統能解此題的一個很重要的特點，他貼文表示此特點必須為後世所師法。
我若不刻意指出這個特點，大家可能也很難注意到，但如果從函數圖來看問題，就能看出問題所在了。函數圖中在 x = 80 的位置，是一個奇點，它的函數值是正的無限大，過去，所有的數學計算系統，在演算過程中，若碰到這種計算，會立刻執行例外處理、發出警告、停止執行。可是，我的系統就故意從這一個奇點開始計算，還進行了兩次，系統碰到了，照樣繼續算，直到產生結果，不需要發出任何警告。我的系統可以接受無限大的數字，而且，這個數字恆就是個無限大，使用者若想印出來看，系統就印出是帶有 + 或 - 符號的 infinit 。
二分法的基本精神，前文已經介紹過，這個範例則專挑學校老師也無法手解的問題來示範二分法的功能，它的發揮，主要是必須先靠繪出函數圖來輔助，您才能確認函數的實根大約會落在那裡？現行電腦的功能，已經能夠非常快速輕鬆的解決函數繪圖問題，求解實根時所需提供給程式使用的上下限二值，根據圖形就很容易取得，剩下來的冗長計算，全部交給電腦。
我在發展到 Lina64 系統之後，有感於這種可用的寶貴程式值得庫存。因此，開始把程式寫得非常整齊，並添加了許多必要的註解，讓自己在日後想要再度使用時，一目了然：
新給的方程式定義，就只需修改 FX 的定義內容，而且，寫法就是日常書寫的格式。
新給的上下限二個值，就只需修改 root1 , root2 , root3 ..... 等等就夠了。
仍然只需執行 main 指令，就能得到所有的答案。
我在 Lina64 中也發展了語音輸出及單音樂譜的執行能力，所以在貼文時，故意展示出這兩項功能。
語音播出是『莫斯科廣播電台的台號音樂與呼叫』，人聲與樂聲均能使用程式輸出。
繪圖的輸出，是藉 Forth 系統控制 Python 程式去操控 GnuＰlot ，是個論秒就能繪出函數圖的設計，那個環境，可以直接操作圖形歸檔。
所有執行程序均在 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)，日後再度使用時，只須在數據輸入區塊內更換數據，就可以執行出整個程式得到結果。這個程式我已用過無數次。
:
這樣的安排，對一般人而言，不太具有什麼意義，對我自己而言，則有深層的意義。
在第(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 函數的係數，作為語音輸出時的指定係數，套裝軟體就能發出各種聲音。
究其原因，就是各種聲音的波譜函數，均能重新表示成傅立葉級數展開式，所以，只要取得級數各項的係數，代入程式，就能令語音合成晶片發出同樣波譜的聲音。
音調的高低，只需設定單個頻率參數就能改變。聲音的大小，只是音量放大的設定。持續的時間，也只是設定共振週期的長度便可辦到。套裝軟體的功能一應俱全，說明資料也很完整，軟體還有許多套，套套都很容易使用，規格也都相當近似。
這麼簡單的語音輸出道理，引起我的興趣，於是就寫了這個程式。為了解釋效果，我也再度利用圖示，顯示理論效果。發現，以此方式完成的工作，簡直就跟原來的波譜完全相似，效果非凡。
換句話說，如果我有心從事於這項工作，絕對能把任何歌星的聲音重現，而且有把握讓聲音以清唱的方式，從語音輸出晶片發送出去。
這一生可以研究的主題很多，不缺這種課題，因為我寫出這個範例程式、看過效果之後，覺得學問並不高深，但要做這方面的事情時，必須花費非常多的時間才能完成。
語音的花樣跟數學題目一樣，是無限多個。有用的一整套規模，至少也得發展出一套中文語音輸出所需的數量，才有值得稱道的地方。我在網文上表示過，這應該是國家政策上的計畫，不是我個人該做的事情。想找齊標準男生與女生的音源，就不是憑我個人能力所能及的事情。還要系統性的造出許多必要的波譜，才能給我這種程式把係數分離出來，提供套裝軟體使用。後續的控制，就更有得設計的了，最後才能成為有用的產品。可是，整個工作，需要用到電腦計算，取得係數的地方，只需像我寫的這個範例程式，就夠用了。別人也許也有更特別的方法，最後，也不外乎取得係數。
這麼說來，請看看我的貼圖，體會一下數學計算系統的能力，把一個很尖銳的週期性三角形波譜，全用傅立葉級數展開式重新表示時的差別，也就夠了。
我們在大學裡學過的傅立葉級數，是能夠直接手解解出三角波的理論係數的，這個範例程式沒有直接取用手解解答所得係數來繪圖，是真正使用電腦的數位解法所得之數值係數繪出的圖形，只取五項，我只在單一個週期內繪圖，就已經能夠有這種效果了。
程式的輸入數據，可以是任意的波譜讀數。點數也沒有限制，點取得多了，解析度就高。所用到的理論式子，大二的工程數學中都有。我因長期從事過信號處理的軟體研發工作，熟悉這方面數學與軟體結合起來的應用技術。
:
這個範例程式的詳細解說，包括引用公式，請參考我個人網頁 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 系統的測試結果列示於後：
請從間隔量為一千點，與間隔量提升為一萬點後的積分結果中，體會所得積分之值，與實際本函數之函數值比較，驗證其精確度是否提高了好幾位數？
:
積分技術的細部討論，已經在第(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 便是。
這個範例是個很特別的設計，要讀過很多資料後，才設計得出來。如果覺得沒有用，請自己印個方陣拆掉一些數字，當個電玩遊戲玩一玩就算了。
:
關於魔術方陣，網上可以找到許多參考文獻，魔術方陣是中國人最先發現的，三階魔術方陣又稱『洛書』，已有兩、三千年的歷史。魔術方陣的定義是把由 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 ***************
訂閱：
張貼留言 (Atom)
沒有留言:
張貼留言