前言:
前面幾篇是來利用最小平方法來配適實驗數據, 求一個多階的多項式. 今天的這一篇, 要倒過來, 假定已經已知一個多階多項式(例如: f(x)=2x^3-4x^2+5x-7), 要如何有效率的求值呢?
就來舉例吧, 假定這個多項式是 f(x)=2x^3-4x^2+5x-7, 要求 x=123 處的值 f(x=123)=?
傳統的算法:
(1)
x=123, 那先算 2x^3 =?
123^2 = 123*123 = 15129
123^3 = 123*123*123 = 15129*123 = 1860867
所以 2x^3 = 2*1860867 = 3721734
(2)
x=123, 接下來算 -4x^2 =?
123^2 = 123*123 = 15129
所以 -4x^2 = -4*15129 = -60516
(3)
x=123, 接下來算 5x =?
123^1 = 123
所以 5x = 5*123 = 615
(4) 全部加起來
f(x=123) = 2*123^3 -4*123^2+5*123-7 = 3721734 -60516 +615 -7 = 3661826
這裡你可以注意到, 當計算 123^3, 123^2, 123^1 的多次方的數值時, 素乎不用重複計算耶. 只要前面的結果拿出來套, 這樣就可以囉!
123^3 = 123^2*123 = 15129*123
19世紀初的英國數學家 威廉·喬治·霍納 (William George Horner), 也發現了這個現象, 因此證明並提出了一個快速計算多項式值的方法. 後來被稱為 Horner's Method 或 Horner's Rule.
再進一步地追朔, 其實13世紀，中國南宋數學家秦九韶在《數書九章》中已經開始使用這樣的方法了, 所以這個方法也被稱之為秦九韶算法. Qin Jiushao's Method (Rule)
比較聰明的算法 (Horner's Method):
f(x)=2x^3-4x^2+5x-7, f(x=123)=?
先把係數列成如下表 (藍色)
| x^3 x^2 x^1 x^0
x=123 | +2 -4 +5 -7
| +246 +29766 +3661833
____________________________________
+2 +242 +29771 +3661826 <-- Ans
計算規則如下, (紅色)
(1) 第一欄 +2 放下來, 寫到最下面第四列 +2
(2) 第一欄最下面乘上 x=123, 2*123 = 246 填寫到第二欄中間第三列
(3) 第二欄第二列係數 -4 跟第三列的數值 +246 相加, 結果寫到最下面第四列 -4 + 246 = +242
(4) 重複剛剛的程序, 最下面第四列的數值乘上 x=123, 放入下一欄中間第三列, 下一欄第二列係數跟第三列的數值相加, 結果寫到最下面第四列.
(5) 重複同樣的程序直到最後, 最後一欄的最下面一列的數值, 即為答案. 所以 f(x=123) = 3661826
為什麼可以這樣計算? 這是什麼魔法?? 我們用代數來推導示範一下!
一樣, 同樣的例子!
f(x)=2x^3-4x^2+5x-7, 將x 不斷提出去, 化簡後變成,
所以, 注意看括弧內的東西, 都是 x乘上( ) + 係數, 所以可以變成一種疊代法的形式囉!
(1) 先讓 b0=2 帶入 b1 = b0*x -4 得到 b1
(2) 同樣的運算步驟, b2 = b1*x +5 得到 b2
(3) 再同樣的程序, b3 = b2*x -7 得到 b3
就這樣, 不斷的將剛剛得到新值乘上 x 加上新的係數項變成下一個新值, 最後就得到解答! 這就是我們剛剛操作的規則, 這個方法, 就叫 Horner's method!
以更一般廣義的數學語言描述如下,
FORTH 程式
因為多項式函數, 不管是在最小平方法配適實驗數據, 或是用來建造一些特殊函數都是非常好用的. 所以我們就來把這個 Horner's method 用 FORTH 語言給實作一下吧!
多項式函數在程式裡面, 可能會需要大量不同的多項式函數, 它們的運作方式都是相同的, 只是裡面的冪次方跟係數不同. 所以這回, 我們來用非常有 FORTH 特色的, FORTH 特有的 create/does> 來實作吧! FORTH 本身不是物件導向的語言, 但是這個 create/does> 卻是有點物件的特性, 它可以寫出一個軟體模板的指令, 使用者可以利用這個指令在程式裡面, 複製出大量具有相同行為跟資料結構但裡面的資料不太相同的軟體元件出來.
create/does> 也是 FORTH 語言裡面重要的一環, 所以這個剛好這是個絕佳的使用範例囉!
語法設計
Horner's Method 本來就是要來計算多項式的值的一種方法. 所以, 我們來用 create/does> 來創造出一個可以創造出任意階多項式的軟體套版跟其語法吧!
PolyDef 用來創造定義一個多項式的新詞
Coef, 用來建構儲存這個多項式裡的係數資料至新詞
PolyEnd 用來結束這個多項式新詞的定義.
使用範例:
定義一個 f(x) = 2 x^3 -4 x^2 +5 x -7 的多項式元件
PolyDef f1(x) 定義一個叫做 f1(x)的新詞
2e0 Coef, x^3的係數
-4e0 Coef, x^2的係數
5e0 Coef, x^1的係數
-7e0 Coef, x^0的係數
PolyEnd 結束定義
這樣就可以定義出一個叫做 f1(x)的多項式函數的詞,
只要在浮點堆疊放入 x 的值, 然後呼叫這個詞, 它會自動算出這個多項式的結果.
例如,
123e0 f1(x) f. 會得到 f1(x=123) = 3661826 的數值.
程式解說:
對初學者, 先複習一下 create/does>
FORTH 的 create does> 語法如下,
: 創造軟體模板的指令
create 從輸入緩衝區取出一個新字串, 在字典裡創造一個新詞, 這個詞執行的時候會把詞的位址丟出來
程式碼, 定義在翻譯時要做什麼事.
does> 編譯一個指令, 讓這個新詞執行時可以跳到這裡之後的程式碼.
程式碼, 定義執行時要做什麼事.
;
我們的語法設計, PolyDef 是這個可以用來創造多項式函數軟體模板的詞, 執行這個詞會創造並建構一個多項式函數出來. 所以 create/does> 的結構放在這裡!
先定義一個特殊旗號 CheckFlag, 來做數據最後收尾時候確認用的
-414159 constant CheckFlag
定義開始時 PolyDef 會先用 create 造新詞, 然後將這個 CheckFlag 旗號推上堆疊, 做個記號. 最後收尾, 所有定義結束時, PolyEnd 會檢查一下這個特殊的旗號 CheckFlag 還在堆疊上嗎? 如果沒抓到它, 這代表堆疊上數字的數目有問題, 有不正常的殘留或短缺, 所以 abort 吧!
PolyDef 會起頭整個多項式新詞的定義, 除了在FORTH 字典裡面創造這個多項式新詞的頭外, 還會留個 CheckFlag 檢查旗號, 跟多項式新詞的起始位址 (here) 在堆疊上. 整個定義的過程, 這兩個資料會一直留在堆疊上. Coef, 除了將係數編入字典外, 還會利用這個堆疊上的這個新詞的起始位址位址來更新多項式最高冪次方的數字. 所以邊編譯係數入字典, 邊數已經編入的係數數目當作最高冪次方的數目. PolyEnd 結束定義時, 做最後的收尾動作, 會利用堆疊上的這兩個資料最最後的總檢查, 確認一下檢查旗號還在嗎, 最高羃次的數字還符合邏輯嗎? 如果都沒有就 abort, 避免粗心的程式設計師犯了一些奇怪的錯誤!
整個 Horner's Method 實作在PolyDef 裡面 does> 之後的程式碼, 所有PolyDef 所創造出來的多項式指令, 都會利用這段程式碼來做 Horner's Method 運算, 根據它們的係數資料, 計算出最後多項式的值.
這裡的程式碼就看程式設計師囉, 大部分的FORTH程式設計師喜歡直接用堆疊操作啦, 不喜歡額外再定義變數, 因為這樣寫出來的程式碼比較簡潔! 當然, 因為堆疊的數字是以堆疊上的順序來被參照的, 所以需要適當的註解來保持一定的可讀性!
: PolyDef
create 翻譯時態,用PolyDef後的字串造個新詞,這個新詞執行的時候會放上它的位址
CheckFlag 放個特殊旗號,數據最後收尾時候確認用的
[ CheckFlag s>d d>f ] fliteral 浮點堆疊也放一個吧,最後收尾時確認用
here 把這個新詞的位址推上堆疊,這個位址存放多項式的階數
0 , ( num)把here這個位址存個零,這是一開始多項式的階數(零)
does> ( addr F: x — result)
執行時,被執行的多項式會先留個位址在堆疊後跳到這裡來執行
後面的程式碼就是用來執行 Homer's Method
dup cell+ 位址往後移一位,這是係數1 最高冪次係數的位置
swap @ 0 迴圈,掃過所有係數 0e0 一開始的 b=0,然後跳入迴圈
do 堆疊狀態: addr' (係數位址) F: x b
fover f* ( addr' F: x b*x) 就 b*x
dup f@ ( addr’ F: x b*x Coef) 取出多項式係數Coef
f+ ( addr’ F: x b*x+Coef) 就 b*x + Coef
float+ 位址往後移一位,指向下一個係數的位址
loop
drop fswap fdrop 結束Homer's Method疊代,丟棄不用的參數,只留最後結果
;
: Coef, 翻譯時態, 用來將係數編入多項式詞裡面. 同時將最高冪次數目加一
f, 將係數編入多項式詞裡面
1 over +! 最高冪次數目加一
;
: PolyEnd 翻譯時態, 結束多項式定義的收尾指令
@ 1 < 取出這個多項式的冪次,小於1的話代表什麼值都沒有,有問題啊!
swap CheckFlag <> 檢查一下,那個特殊檢查旗號還在嗎,不見了就是有問題的!
[ CheckFlag s>d d>f ] fliteral f= invert 也檢查一下浮點堆疊吧!
or or 三個條件有任一符合的話就 abort 吧!
abort" Errror... Something went wrong..."
;
不廢話, 最後執行的結果!
三個範例:
第一個, 就是我們的舉例 f(x) = 2x^3 -4x^2 +5x -7
PolyDef f1(x)
2e0 Coef,
-4e0 Coef,
5e0 Coef,
-7e0 Coef,
PolyEnd
以 x=-2 來確認, 結果是 f(x=-2) = -49
以 x=123 來確認, 結果是 f(x=123) = 3661826
兩個答案都正確無誤!
第二個, f(x) = 2x^4 -3x^2 +3x -4
PolyDef f2(x)
2e0 Coef,
0e0 Coef,
-3e0 Coef,
3e0 Coef,
-4e0 Coef,
PolyEnd
以 x=-2 來確認, 結果是 f(x=-2) = 10
這個結果也是正確無誤的!
最後一個, 回到我們最小平方法多項式配適的那篇 BLOG, 我們配適出來撞擊能量 Energy 跟溫度T的關係. E(T) = -7.12085781E-7 T^4 -7.04053872E-5 T^3 +0.0103986946 T^2 +1.46923973 T +49.2061305
所以來定義它囉, 這樣隨時就可以很方便計算!
PolyDef Energy(T) ( F:Temperature -- F:Energy)
-7.12085781E-7 Coef,
-7.04053872E-5 Coef,
+0.0103986946E0 Coef,
+1.46923973E0 Coef,
+49.2061305E0 Coef,
PolyEnd
要計算 T=25 DegreeC 下, 撞擊能量多少呢?? 只要在終端機鍵入下面指令,
25E0 Energy(T) f. 馬上回答 91.058065...
為什麼要這麼麻煩?
最後來解答, 為什麼要使用 Horner's method?? 不是直接計算就好囉!!
原來, 是計算量的問題! 來計算一下, 一般的計算法跟 Horner's method 兩者計算量的差異吧!
首先, 電腦裡面 CPU 做加減法運算跟做乘除法法運算, 所耗的時間不太一樣. 乘除法比較耗 CPU 的時間, 假設這是 1 比 5, 也就是假設加減法運算需要1個單位時間, 乘除法法運算需要5個單位時間.
(+), (-): 1個單位時間, (*), (/): 5個單位時間. 假定一個n階冪次方的多項式
1. 直接運算法
乘法 (*) 運算的次數: (n+1) + n + (n-1) + ... + 2 + 1 = (n+1+1)*(n+1)/2 = (n+2)(n+1)/2
加法 (+) 運算的次數: 1 + 1 + .... + 1 = n
所以, 總運算時間 = 乘法 (*) 運算的次數 * 5 + 加法 (+) 運算的次數 * 1 = 5(n+2)(n+1)/2 + n = (5n^2 +17n+5)/2 = 2.5n^2 + 8.5n +2.5 ~ O(n^2)
2. Horner's method
乘法 (*) 運算的次數: 1 + 1 + .... + 1 = n
加法 (+) 運算的次數: 1 + 1 + .... + 1 + 1 = n+1
所以, 總運算時間 = 乘法 (*) 運算的次數 * 5 + 加法 (+) 運算的次數 * 1 = 5n + (n+1) = 6n+1 ~ O(n)
直接運算法所花的時間 order 是 O(n^2), 所以是以平方的形式隨著多項式的階數增加而瞬間暴增.
Horner's method 所花的時間是 O(n), 是隨著多項式的階數增加而線性和緩增加.
兩者所花的時間真的天差地遠呀! 所以當然要用 Horner's method 這種聰明的方法來計算多項式的值囉!
原始程式列表:
\
\ horner's method for polynomial
\ Frank Lin 2018.06.01
\
-414159 constant CheckFlag
: PolyDef
create CheckFlag [ CheckFlag s>d d>f ] fliteral
here 0 , ( num)
does> ( addr F: x — result)
dup cell+
swap @ 0 0e0 ( addr’ F: x b )
do fover f* ( addr' F: x b*x)
dup f@ ( addr’ F: x b*x Coef)
f+ ( addr’ F: x b*x+Coef)
float+ ( step to next coef)
loop
drop fswap fdrop
;
: Coef,
f, 1 over +!
;
: PolyEnd
@ 1 < ( no data)
swap CheckFlag <> ( Checkflag disappear)
[ CheckFlag s>d d>f ] fliteral f= invert ( Checkflag disappear on f stack)
or or
abort" Errror... Something went wrong..."
;
\
\ Example Codes
\
\
\ ex1: f1(x) = 2x^3 -4x^2 +5x -7
\ f(x=-2) = -49
\
true
[if]
PolyDef f1(x) ( F:x -- F:y)
2e0 Coef,
-4e0 Coef,
5e0 Coef,
-7e0 Coef,
PolyEnd
.( f1 = 2x^3 -4x^2 +5x -7, when x=-2, it's ) -2e0 f1(x) f. cr
.( correct answer is -49.)
cr cr
.( f1 = 2x^3 -4x^2 +5x -7, when x=123, it's ) 123e0 f1(x) f. cr
.( correct answer is 3661826.)
cr cr
\
\ ex2: f2(x) = 2x^4-3x^2+3x-4
\ f(x=-2) = 10
\
PolyDef f2(x) ( F:x -- F:y)
2e0 Coef,
0e0 Coef,
-3e0 Coef,
3e0 Coef,
-4e0 Coef,
PolyEnd
.( f2 = 2x^4-3x^2+3x-4, when x=-2, it's ) -2e0 f2(x) f. cr
.( correct answer is 10.)
cr cr
\
\ ex3:
\ E = -7.12085781E-7 T^4 -7.04053872E-5 T^3 +0.0103986946 T^2 +1.46923973 T +49.2061305
\
\
PolyDef Energy(T) ( F:Temperature -- F:Energy)
-7.12085781E-7 Coef,
-7.04053872E-5 Coef,
+0.0103986946E0 Coef,
+1.46923973E0 Coef,
+49.2061305E0 Coef,
PolyEnd
.( Energy when T=25 DegreeC, it's ) 25e0 Energy(T) f. cr
.( correct answer is 91.058065 )
cr cr
[then]
xxx
留言列表
