\ gamma The Gamma, loggamma and reciprocal gamma functions
\ Calulates Gamma[x], Log[ Gamma[x] ] and 1/Gamma[x] functions
\ (for real arguments)
\ Forth Scientific Library Algorithm #18
\ This is an ANS Forth program requiring:
\ 1. The Floating-Point word set
\ 2. The word '}' to dereference a one-dimensional array.
\ 3. Uses the word '}Horner' for fast polynomial evaluation.
\ 4. The FCONSTANT PI (3.1415926536...)
\ 5. The words 'S>F' and 'F>S' to convert between floats and singles
\ 6. The word F>
\ : F> FSWAP F< ;
\ 7. The compilation of the test code is controlled by the
\ the conditional compilation words in the Programming-Tools wordset.
\ Baker, L., 1992; C Mathematical Function Handbook, McGraw-Hill,
\ New York, 757 pages, ISBN 0-07-911158-0
\ The reciprocal Gamma function is ACM Algorithm #80
\ 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 this
\ copyright notice is preserved.
\ ( GAMMA V1.2 6 October 1994 EFC )
\ include lib/ansfloat.4th
\ include lib/fsl-util.4th
[UNDEFINED] gamma [IF]
[UNDEFINED] }horner [IF] include lib/horner.4th [THEN]
[UNDEFINED] fsin [IF] include lib/fsinfcos.4th [THEN]
[UNDEFINED] fln [IF] include lib/flnflog.4th [THEN]
[UNDEFINED] fexp [IF] include lib/fexp.4th [THEN]
[UNDEFINED] f> [IF] include lib/fequals.4th [THEN]
9 FLOAT [*] 1 [+] ARRAY b{
4 FLOAT [*] 1 [+] ARRAY ser{
14 FLOAT [*] 1 [+] ARRAY b-inv{
FLOAT b{ FSL-ARRAY
FLOAT ser{ FSL-ARRAY
FLOAT b-inv{ FSL-ARRAY
:THIS b{ DOES> (FSL-ARRAY) ;
:THIS ser{ DOES> (FSL-ARRAY) ;
:THIS b-inv{ DOES> (FSL-ARRAY) ;
FLOAT ARRAY X-TMP \ scratch space to be kind on the fstack
FLOAT ARRAY Z-TMP
: logsr2pi s" 0.918938533" s>float ;
: init-b-ser
1 s>f b{ 0 } F! S" -0.577191652" s>float b{ 1 } F!
S" 0.988205891" s>float b{ 2 } F! S" -0.897056937" s>float b{ 3 } F!
S" 0.918206857" s>float b{ 4 } F! S" -0.756704078" s>float b{ 5 } F!
S" 0.482199394" s>float b{ 6 } F! S" -0.193527818" s>float b{ 7 } F!
S" 0.035868343" s>float b{ 8 } F!
S" 0.08333333333333" s>float ser{ 0 } F!
S" -0.002777777777" s>float ser{ 1 } F!
S" 0.000793650793" s>float ser{ 2 } F!
S" -0.000595238095" s>float ser{ 3 } F!
1 s>f b-inv{ 0 } F! S" -0.422784335092" s>float b-inv{ 1 } F!
S" -0.233093736365" s>float b-inv{ 2 } F! S" 0.191091101162" s>float b-inv{ 3 } F!
S" -0.024552490887" s>float b-inv{ 4 } F! S" -0.017645242118" s>float b-inv{ 5 } F!
S" 0.008023278113" s>float b-inv{ 6 } F! S" -0.000804341335" s>float b-inv{ 7 } F!
S" -0.000360851496" s>float b-inv{ 8 } F! S" 0.000145624324" s>float b-inv{ 9 } F!
S" -0.000017527917" s>float b-inv{ 10 } F! S" -0.000002625721" s>float b-inv{ 11 } F!
S" 0.000001328554" s>float b-inv{ 12 } F! S" -0.000000181220" s>float b-inv{ 13 } F!
;
: non-negative-x ( -- , f: x -- loggamma{x} )
FDUP 1 s>f F> IF
FDUP 2 s>f F> IF
X-TMP F!
1 s>f X-TMP F@ F/
FDUP Z-TMP F! FDUP F*
ser{ 3 }Horner Z-TMP F@ F*
logsr2pi F+ X-TMP F@ F-
X-TMP F@ FLN
X-TMP F@ s" 0.5" s>float F- F*
F+
ELSE
1 s>f F- b{ 8 }Horner FLN
THEN
ELSE
FDUP F0= 0= IF
FDUP X-TMP F!
b{ 8 }Horner
X-TMP F@ F/ FLN
THEN
THEN
;
: ?negative-int ( -- t/f , f: x -- x )
\ check to see if x is a negative integer, or zero
FDUP F0< IF
FDUP FDUP F>S S>F F- F0=
ELSE
FDUP F0=
THEN
;
: rgam ( -- , f: x -- z )
FDUP
b-inv{ 13 }Horner
FOVER 1 s>f F+ F* F*
;
: rgam-large-x ( -- , f: x -- z )
1 s>f \ the AA loop
BEGIN
FSWAP 1 s>f F-
FSWAP FOVER F*
FOVER 1 s>f F> 0=
UNTIL
FOVER 1 s>f F= IF FSWAP FDROP 1 s>f FSWAP F/
ELSE
FSWAP rgam FSWAP F/
THEN
;
: rgam-small-x ( -- , f: x -- z )
FDUP -1 s>f F= IF FDROP 0 s>f
ELSE
FDUP -1 s>f F> IF rgam
ELSE
FDUP \ the CC loop
BEGIN
FSWAP 1 s>f F+
FDUP -1 s>f F<
WHILE
FSWAP
FOVER F*
REPEAT
rgam F*
THEN
THEN
;
\ Log Gamma function
: loggam ( -- ) ( f: x -- loggamma{x} ) \ input arg is returned if routine aborts
\ check to make sure x is not a negative integer or zero
?negative-int ABORT" loggam has 0 or negative integer argument"
FDUP F0< IF
FABS 1 s>f F+ Z-TMP F!
PI Z-TMP F@ F* FSIN FABS PI FSWAP F/ FLN
Z-TMP F@
non-negative-x
F-
ELSE
non-negative-x
THEN
;
\ Gamma function
: gamma ( -- ) ( f: x -- g{x} ) \ input arg is returned if routine aborts
\ check to make sure x is not a negative integer or zero
?negative-int ABORT" gamma has 0 or negative integer argument"
FDUP loggam FEXP
FOVER F0< IF
FOVER F>S NEGATE 2 MOD
2* 1- S>F F*
THEN
FSWAP FDROP
;
: rgamma ( -- ) ( f: x -- 1/g{x} ) \ reciprocal gamma function
FDUP F0= FDUP 1 s>f F= OR 0= \ will return x if x is zero or one
IF
FDUP 1 s>f F< IF rgam-small-x
ELSE
rgam-large-x
THEN
THEN
;
fclear
init-b-ser
[DEFINED] 4TH# [IF]
hide b{
hide ser{
hide b-inv{
hide Z-TMP
hide X-TMP
hide logsr2pi
hide init-b-ser
hide non-negative-x
hide ?negative-int
hide rgam
hide rgam-large-x
hide rgam-small-x
[THEN]
false [IF] \ test code =============================================
: gamma-test ( -- )
fclear
10 set-precision
CR
." gamma( 5.0 ) = " 5 s>f gamma F. ." (should be 24.0) " CR
." gamma( -1.5 ) = " s" -1.5" s>float gamma F. ." (should be 2.36327) " CR
." gamma( -0.5 ) = " s" -0.5" s>float gamma F. ." (should be -3.54491) " CR
." gamma( 0.5 ) = " s" 0.5" s>float gamma F. ." (should be 1.77245) " CR
CR
." rgamma( 5.0 ) = " 5 s>f rgamma F. ." (should be 0.0416667) " CR
." rgamma( -1.5 ) = " s" -1.5" s>float rgamma F. ." (should be 0.423142) " CR
." rgamma( -0.5 ) = " s" -0.5" s>float rgamma F. ." (should be -0.282095) " CR
." rgamma( 0.5 ) = " s" 0.5" s>float rgamma F. ." (should be 0.564190) " CR
;
gamma-test
[THEN]
[THEN]
|
