2026年10月2日 星期五

高級數學函數

高級數學函數


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


本文記錄高級數學函數在ABC FORTH數學計算系統中的發展方法、設計格式、以及一些相關問題的探討,是一篇實用教材。ABC FORTH數學計算系統讓這類問題容易研究發展,也容易實際應用。

相對於高級數學函數者,則被稱作『一般數學函數』。在『ABC FORTH數學計算系統使用說明』第17頁中,提及一項重要觀念,它是設計出這個系統所依據的基本原理。數學式子內,除了數字、括弧、四則運算外、其他符號全被視同為一般數學函數來處理,包括等號『=』亦然,sine與cosine則是大家耳熟能詳的最典型一般數學函數。既然根據這樣的原理,才容易設計出這個系統,繼續延伸應用同一原理,就能大量擴充系統的功能。換句話說,此後所有發展出來的高級數學函數,只要設計格式與系統規格要求一致,就都可以輕而易舉的納入ABC FORTH數學計算系統,而成為系統的一個函數。因此,系統的發展是無限的,完全符合FORTH的基本哲理。

函數是用來描述變數間關係的數學表示式(Function is a mathematical expression describing the relation between variables.)。數學計算系統必須提供直接使用函數的功能,才能方便的設計程式,函數的數量卻是無限的。因此,除了一般數學函數以外的所有函數,都可以稱之為高級數學函數,ABC FORTH數學計算系統則根基本Forth系統建立了一般數學函數的功能。

現行個人電腦中,一般數學函數的計算,都是由具有浮點功能的數學處理器(Coprocessor)硬體元件協助完成的,硬體不能造的太大,所以功能就有所限制,通常均以提供一般數學函數為主,ABC FORTH數學計算系統的原始功能也僅限於此。近代與將來的電腦,則會根據使用需求,設計出許多具有高級數學函數功能的硬體元件,例如:為影像處理或高速信號處理的需要,而增加之特別函數功能加速處理器(Accelerator boards),增添的功能也只能設計出多少就加多少,永遠都還是有無限多的高級數學函數,不能以硬體實現其功能。ABC FORTH數學計算系統具備了健全的系統規劃,可以直接了當的添加任何新增硬體功能,因此,我們不討論如何增加硬體性的高級數學函數問題,只討論軟體性的程式設計技術,這些技術永遠需要,永遠有效。

適合以電腦處理的高級數學函數,是工程或研究方面需要使用、不易人工手解、需要大量數學計算的函數,它們需要電腦協助,電腦才予以協助,並不是所有高級數學函數,都需要設計程式實現功能,有些函數可以直接手解得到答案,不須藉助於電腦。本文僅舉一個不太可能會被考慮使用硬體實現其功能,卻又非得使用電腦計算不可的高級數學函數,當作程式設計範例,它被稱為指數積分函數(Exponential integral function),型式如下:

y=∫(exp(-x)/x)dx 積分的範圍從任意x值到無限大

與此高級數學函數相類似,且在工程應用與科學研究方面,普及性需求相當高的高級數學函數有許多種,而且在電腦中被處理的方式也完全相近,例如:

(1) gamma function
(2) Bessel function
(3) Cosine and sine integrals
(4) Fresnel integrals
(5) Legendre polynomial
(6) confluent hyper-geometric function
(7) hyper-geometric functions
(8) Hermite polynomial
(9) Laguerre polynomial
(10)Chebyshev polynomial

以指數積分函數當作基礎,透過觸類旁通的方式,利用ABC FORTH數學計算系統的便利操作性能,可以迅速設計出這些高級數學函數的功能程式。因此,ABC FORTH系統除了可以協助使用者更為深入的了解這些函數,或核驗這些程式的品質及準確性外,還可以進一步提供有興趣的人,以ABC FORTH系統設計出更新、更好、更可靠的實用程式。

Fig總會網頁的科學庫存程式(Forth Scientific Library Algorithm)中,第一個可以自由下載的純FORTH程式,就是指數積分高級數學函數,它具有高級數學函數的全面代表性,我們把原程式轉置此處,再仔細討論與比較。

 
\ expint     Real Exponential Integral         ACM Algorithm #20

\ Forth Scientific Library Algorithm #1

\ Evaluates the Real Exponential Integral,
\     E1(x) = - Ei(-x) =   int_x^\infty exp^{-u}/u du      for x > 0
\ using a rational approximation

\ This code conforms with ANS requiring:
\      1. The Floating-Point word set
\      2. The immediate word '%' which takes the next token
\         and converts it to a floating-point literal
\ 

\ Collected Algorithms from ACM, Volume 1 Algorithms 1-220,
\ 1980; Association for Computing Machinery Inc., New York,
\ ISBN 0-89791-017-6

\ (c) Copyright 1994 Everett F. Carter.  Permission is granted by the
\ author to use this software for any application provided the
\ copyright notice is preserved.

CR .( EXPINT     V1.1                  21 September 1994   EFC )

: expint ( --, f: x -- expint[x] )
        FDUP
        % 1.0 F< IF
                    FDUP % 0.00107857 F* % 0.00976004 F-
                    FOVER F*
                    % 0.05519968 F+
                    FOVER F*
                    % 0.24991055 F-
                    FOVER F*
                    % 0.99999193 F+
                    FOVER F*
                    % 0.57721566 F-
                    FSWAP FLN F-
                ELSE
                    FDUP % 8.5733287401 F+
                    FOVER F*
                    % 18.059016973 F+
                    FOVER F*
                    % 8.6347608925 F+
                    FOVER F*
                    % 0.2677737343 F+

                    FOVER
                    FDUP % 9.5733223454 F+
                    FOVER F*
                    % 25.6329561486 F+
                    FOVER F*
                    % 21.0996530827 F+
                    FOVER F*
                    % 3.9584969228 F+

                    FSWAP FDROP
                    F/
                    FOVER F/
                    FSWAP % -1.0 F* FEXP
                    F*
                THEN
;


\ test code generates a small table of selected E1 values.
\ most comparison values are from Abramowitz & Stegun,
\ Handbook of Mathematical Functions, Table 5.1

: expint_test ( -- )

        CR
        ."   x    E1(x) exact      ExpInt[x] " CR

      ."  0.5   0.5597736      "
      % 0.5 expint  F. CR

      ."  1.0   0.2193839      "
      % 1.0 expint   F. CR

      ."  2.0   0.0489005      "
      % 2.0 expint    F. CR

      ."  5.0   0.001148296    "
      % 5.0 expint    F. CR

      ." 10.0   0.4156969e-5   "
      % 10.0 expint   F. CR

;


這個高級數學函數程式在Win32Forth系統中可以執行,但在程式起頭處要自行設計一個%指令:

: % BL WORD COUNT >FLOAT NOT ABORT" not a fp#"
STATE @ IF POSTPONE FLITERAL THEN ; IMMEDIATE

程式的作者是Everett F. Carter,年代是1994年,算是近代新作品,現在還很實用,大家都應該感謝原作者不吝公開個人技術的作風,不可批評程式中的任何不是之處,否則您大可不必使用這個程式。因此,本文以下的論述,凡涉及測試數據檢驗、程式改善比較、資料考古查證…等,與這個程式記錄內容有所違背之處,均不代表對原作者的不敬,反而更應該感謝這位作者,賜給我們可以仔細探討這個問題的機會。

我抱著追根究底的精神遍尋資料,研究這個高級數學函數。在圖書館找資料,追蹤指數積分函數問題,結果發現設計這個程式的最早依據,是來自於更早的數學技術參考文獻,而且程式中原原本本的使用所有的原始常數係數。

1975年出版的計算器科學分析『Scientific analysis on the pocket calculator』,作者:Jon M. Smith,書中第四章明確的列出了設計上列程式時,所須要的兩個式子,記錄如下:



   
以E(x)=∫(exp(-x)/x)dx表示指數積分高級數學函數
當x介於0與1之間或x < 1時 E(x)+ln(x) = a(0)+x(a(1)+x(a(2)+x(a(3)+x(a(4)+a(5)x)))) + eps(x) a(0)= -0.57721566 a(3)= 0.05519968 a(1)= 0.99999193 a(4)= -0.00976004 a(2)= -0.24991055 a(5)= 0.00107857 eps(x) < 2/10^7 當x大於或等於1亦即x >= 1時 x*exp(x)*E(x)=( (a(4)+x(a(3)+x(a(2)+x(a(1)+x)))) /(b(4)+x(b(3)+x(b(2)+x(b(1)+x)))) ) + eps(x) a(1)= 8.5733287401 b(1)= 9.573322454 a(2)= 18.0590169730 b(2)= 25.6329561486 a(3)= 8.6347608925 b(3)= 21.0996530827 a(4)= 0.2677737343 b(4)= 3.9584969228 eps(x) < 2/10^8


找到上列數學式子,您就立刻明白指數積分高級函數FORTH程式是如何設計的了。進一步仔細查核上列兩個有理多項近似式的來源,書中參考資料詳實的記載為美國國家標準局出版的數學函數使用手冊:the Handbook of Mathematical Functions, U.S. Department of Commence, National Bureau of Standards, Applied Mathematics Series 55, 1900.

我對此書出版的年代『1900』存疑,可能是『1960』的手植之誤,但沒有能力追查到這本使用手冊予以證實。我們可以合理的懷疑,19世紀以前,可能還沒有能力算出那麼多個小數點後面十位數的精確數字,1960年附近才有可能,因為,計算機技術,1960年代以後才開始崛起。

Everett F. Carter心地善良,免費提供大家好用的程式,被引用的紐約計算機協會(ACM:Association for Computing Machinery Inc., New York),只是一個換了包裝後的單位,書籍雖是1980年出版,實際技術出處則仍來自於更早的數學研究結果。

如果要直接改寫傳統的FORTH程式成ABC FORTH系統要求的格式,是一件非常痛苦的事情,所以我才會追根究底的查找技術的原始出處,找到之後仔細閱讀,然後,只花了半小時,就用ABC FORTH數學計算系統設計完成了下列完全一樣的程式:


 

REAL X REAL Y REAL RES REAL AUX : RESET-VAR {{ X = 0 }} {{ Y = 0 }} {{ RES = 0 }} {{ AUX = 0 }} ; : (EXPINT2) BASIC 10 IF { X > 1 } THEN 40 20 LET { RES = ( -0.57721566 + X * ( 0.99999193 + X * ( -0.24991055 + X * ( 0.05519968 + X * ( -0.00976004 + X * 0.00107857 ) ) ) ) ) - LN ( X ) } 30 GOTO 60 40 LET { AUX = ( 0.2677737343 + X * ( 8.6347608925 + X * ( 18.0590169730 + X * ( 8.5733287401 + X ) ) ) ) / ( 3.9584969228 + X * ( 21.0996530827 + X * ( 25.6329561486 + X * ( 9.573322454 + X ) ) ) ) } 50 LET { RES = AUX / ( X * EXP ( X ) ) } 60 END ; : EXPINT2 ( --, f: x – EXPINT2[x] ) RESET-VAR ADDRESS-OF X F! (EXPINT2) RES ;


您看完對比如此鮮明的兩個程式,應該可以感受出ABC FORTH數學計算系統的特別之處。

Fig總會科學庫存程式提供的資源絕對是有價值的,參考那些程式的設計規格,讓我知道自己設計的程式最終格式應該如何?所以,上列EXPINT2指令所設計的內容雖然很簡單,卻代表了將來可以被引用的一種標準格式,指令必須被設計成像標準sine函數那樣的使用格式,計算前,堆疊必須先放好數字,計算後,所得結果,也必須放在堆疊上。FORTH程式永遠不能忘記堆疊是傳遞參數的標準場所,ABC FORTH系統雖然設計出了BASIC式數學計算程式的使用便利性,但卻仍然不能忽略堆疊。

本文論述的指數積分高級數學函數,如果去除積分,就會成為一個相當整齊的數學函數,例如:寫成y=exp(-x)/x。這樣的函數可以算是一個非線性一階微分方程式的解,此微分方程式為:x*(dy/dx)+y –y*ln(x*y)=0。

所有y=exp(cx)/x型式的函數都是這個微分方程式的解,其中c為任意常數。關於這一段敘述,以及詳細的探討,並附有插圖的文獻,您可以從1971年出版的『Elementary Differential Equations』,作者:Donald L. Kreider, Robert G, Kuller, Donald R, Ostberg,第8章p.328中找到。看了y的曲線圖,您立刻可以明白,它的積分沒有辦法手解,因為函數不連續,當x=0與很大的值時,函數曲線出現了正與負的無限大情況,所以是一個非解析性的函數,不能手解。

查積分函數表,例如:1966年出版的『Tables of integrals and other mathematical data』,作者:Herbert Bristol Dwight,p.134中列出的積分結果為:

∫(exp(ax)/x)dx=Ln|x|+a*x/1!+(a^2*x^2)/(2*2!)+(a^3*x^3)/(3*3!)+…+(a^n*x^n)/(n*n!)+……

它真的是一個沒有辦法用手算得到答案的問題,所以才求助於電腦。古時候的數學家研究這些問題時確實辛苦,上述的說明要查許多本書才寫得出來。

ABC FORTH數學計算系統易於處理這些問題,首先,就利用『函數曲線繪圖程式』,繪製這個指數函數的曲線圖,正或負的無限大對程式繪圖功能完全沒有影響,設計的程式FX6,與執行操作的方式為:

: FX6 {{ Y = EXP ( NEGATE X ) / X }} ;
' FX6 IS F(X)

您可以自行調整顯示所想要的函數曲線圖,此處直接附上我的操作結果,我所選擇的繪圖設定範圍:x為-6到6,y為-66到66,66大順,所以都用6,函數曲線如圖所示。

從顯示的圖形直接可以看出,在x為0時,此函數確實是不連續而為非解析性函數,其積分也就當然無法手解獲得。就幾何意義而言,所謂函數的積分,係表示函數曲線與x軸之間所圍住的面積,圖20.1顯示在x為正值時,積分的結果可得正的結果,但在x為負值時,積分的結果就會是負值,這也就是為什麼fig總會提供的FORTH程式中說明只適用於x > 0 的原因,事實上,x < 0的部份也可積分,只不過積分所得應該是負值罷了。

y = exp(-x)/x 函數曲線圖

如果要擴充這個高級數學函數成x全域使用範圍,上述的兩個程式都辦不到,必須另覓資源。我在1970年版的IBM System/360 Scientific Subroutine Package中找到這樣的應用程式,Programmer’s Manual的書碼是GH20-0205-4,次程式的名稱為:SUBROUTINE EXPI(X,RES,AUX)。改寫成ABC FORTH系統可以執行的程式,結果如下:


 

: (EXPI)   BASIC
10 IF { X >= 1 } THEN 30
20 GOTO 70
30 LET { Y = 1 / X }
40 LET { AUX = 1 - Y * ( ( ( Y + 3.377358 ) * Y + 2.052156 ) * Y
   + 0.2709479 ) / ( ( ( ( Y * 1.072553 + 5.716543 ) * Y + 6.945239 ) * Y
   + 2.593888 ) * Y + 0.2709456 ) }
50 LET { RES = AUX * Y * EXP ( NEGATE X ) }
60 GOTO 9999
70 IF { X > -3 } THEN 90
80 GOTO 150
90 LET { AUX = ( ( ( ( ( ( ( 7.122452E-7 * X - 1.766345E-6 ) * X
   + 2.928433E-5 ) * X - 2.335379E-4 ) * X + 1.664156E-3 ) * X
   - 1.041576E-2 ) * X + 5.555682E-2 ) * X - 2.500001E-1 ) * X
   + 0.9999999 }
100 LET { RES = -1E75 }
110 IF { X <> 0 } THEN 130
120 GOTO 140
130 LET { RES = X * AUX - LN ( ABS ( X ) ) - 0.5772157 }
140 GOTO 9999
150 IF { X > -9 } THEN 170
160 GOTO 190
170 LET { AUX = 1 - ( ( ( ( 5.176245E-2 * X + 3.061037 ) * X + 32.43665 )
    * X  + 2.244234E2 ) * X + 2.486657E2 ) / ( ( ( ( X + 3.995161 ) * X
    + 3.893944E1 ) * X + 2.263818E1 ) * X + 1.807837E2 ) }
180 GOTO 9999
190 LET { Y = 9 / X }
200 LET { AUX = 1 - Y * ( ( ( Y + 7.659824E-1 ) * Y - 7.271019E-1 ) * Y
    - 1.080653 ) / ( ( ( ( Y * 2.518750 + 1.122927E1 ) * Y + 5.921405 )
    * Y - 8.668702 ) * Y - 9.724216 ) }
210 LET { RES = AUX * EXP ( NEGATE X ) / X }
9999 END
;

: EXPI ( -- , f: x - - expint[x] )     \ 注意:不可以用expint(x)當註解。
  RESET-VAR
  ADDRESS-OF X F!
  (EXPI)
  RES  ;


以上所有程式都可以執行,彼此可以交互驗證所得結果正確與否?或精確度差多少?在查究許多相關書籍的同時,我還獲得了幾個更進一步,細部分段計算的數學近似有理多項式,但不具有可以普及使用的意義,所以就不用來設計程式。

根據積分手冊顯示,指數積分函數還有它的一個大家族:
如果最幼小的函數寫成:
E1(x)=∫(exp(-x)/x)dx
那麼他的所有哥哥、姐姐們就可以寫成:
En(x)=∫(exp(-x)/x^n)dx
而兄弟姐妹之間的關係是:
En+1(x)=(exp(-x)-x*En(x))/n 您能不能使用ABC FORTH數學計算系統算一下,
E10(x)=∫(exp(-x)/x^10)dx
當x = 0.0001的時候,這個高級數學函數值應該是多少?
這並不是一個毫無意義的練習題,是相關基本性能延伸性的應用。
ABC FORTH數學計算系統也可以用來發展找尋有理多項近似式,以求解類似高級數學函數的解答,而且更為快速且方便。無法直接積分的非解析性函數,一旦可由有理多項近似式取代表示時,多項式的積分就可以直接用手解而得到答案。例如:多項式X^3-4*X^2+X-5的積分,直接就可以寫成
(X^4/4)-(4*X^3/3)+(X^2/2)-5X。
積分值變成可以直接計算出來,這就是為什麼前人喜歡使用有理多項近似式來取代非可解析性函數的原因。

常用的高級數學函數,通常都是從一些特別、但很有用的微分方程式解出來的答案,它們分別在各個工程實務或科學研究的領域具有重要意義,而且被使用時一定要計算出函數之值,單用手算卻又非常困難,甚至於有的時候還得知道整個高級數學函數的變化趨勢,或全面性的分佈狀況,才足以完整解釋或說明它所代表的科學現象。本文以單一簡例詳細說明了研究此類問題的現行狀況,也是其他所有高級數學函數程式的設計典範,這是一篇健全且實用的教材,舉一反三就能設計出其他所有的高級數學函數,我們可以根據ABC FORTH系統,建立起像美國國家標準局一樣,但純屬於自己的全套系統,不必抄用上述那些固定係數,還有可能做得比他們更好。

ABC FORTH數學計算系統很容易繪製函數曲線圖,也很容易將許多個函數曲線同繪在一個圖內,您就很容易進行圖形比較,也就容易找到品質較佳的有理多項近似式。非解析性函數雖不能手解求得積分,未積分前的函數值卻是可以直接計算而得的,因此,直接算出函數值後,供作求解假設多項式問題的所需數據,有許多方法,可以用來獲得多項式的係數。所得結果的優劣,取決於全面適用範圍的誤差有多少?透過看圖方式來驗證,可以產生直接性的感覺,當然要比尋找計算誤差的理論式子要快,也當然比一點一點計算的結果要好。有興趣的讀者,可以使用這套系統進行這種發展,得到好的結果公諸於世,就是新的貢獻。