\ 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 immediate word '%' which takes the next token \ and converts it to a floating-point literal, \ : % BL WORD COUNT >FLOAT 0= ABORT" NAN " \ STATE @ IF POSTPONE FLITERAL THEN ; IMMEDIATE \ 3. Uses words 'Private:', 'Public:' and 'Reset_Search_Order' \ to control the visibility of internal code. \ 4. Uses the words 'REAL*4' and ARRAY to create floating point arrays. \ 5. The word '}' to dereference a one-dimensional array. \ 6. Uses the word '}Horner' for fast polynomial evaluation. \ 7. The FCONSTANT PI (3.1415926536...) \ 8. The words 'S>F' and 'F>S' to convert between floats and singles \ 9 The word F> \ : F> FSWAP F< ; \ 10. The compilation of the test code is controlled by the VALUE TEST-CODE? \ and 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. cr .( GAMMA V1.2 6 October 1994 EFC ) Private: 9 FLOAT ARRAY b{ 4 FLOAT ARRAY ser{ 14 FLOAT ARRAY b-inv{ FVARIABLE X-TMP \ scratch space to be kind on the fstack FVARIABLE Z-TMP % 0.918938533 FCONSTANT logsr2pi : init-b-ser % 1.0 b{ 0 } F! % -0.577191652 b{ 1 } F! % 0.988205891 b{ 2 } F! % -0.897056937 b{ 3 } F! % 0.918206857 b{ 4 } F! % -0.756704078 b{ 5 } F! % 0.482199394 b{ 6 } F! % -0.193527818 b{ 7 } F! % 0.035868343 b{ 8 } F! % 0.08333333333333 ser{ 0 } F! % -0.002777777777 ser{ 1 } F! % 0.000793650793 ser{ 2 } F! % -0.000595238095 ser{ 3 } F! % 1.0 b-inv{ 0 } F! % -0.422784335092 b-inv{ 1 } F! % -0.233093736365 b-inv{ 2 } F! % 0.191091101162 b-inv{ 3 } F! % -0.024552490887 b-inv{ 4 } F! % -0.017645242118 b-inv{ 5 } F! % 0.008023278113 b-inv{ 6 } F! % -0.000804341335 b-inv{ 7 } F! % -0.000360851496 b-inv{ 8 } F! % 0.000145624324 b-inv{ 9 } F! % -0.000017527917 b-inv{ 10 } F! % -0.000002625721 b-inv{ 11 } F! % 0.000001328554 b-inv{ 12 } F! % -0.000000181220 b-inv{ 13 } F! ; init-b-ser : non-negative-x ( -- , f: x -- loggamma{x} ) FDUP % 1.0 F> IF FDUP % 2.0 F> IF X-TMP F! % 1.0 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@ % 0.5 F- F* F+ ELSE % 1.0 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-integer ( -- 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.0 F+ F* F* ; : rgam-large-x ( -- , f: x -- z ) % 1.0 \ the AA loop BEGIN FSWAP % 1.0 F- FSWAP FOVER F* FOVER % 1.0 F> 0= UNTIL FOVER % 1.0 F= IF FSWAP FDROP % 1.0 FSWAP F/ ELSE FSWAP rgam FSWAP F/ THEN ; : rgam-small-x ( -- , f: x -- z ) FDUP % -1.0 F= IF FDROP % 0.0 ELSE FDUP % -1.0 F> IF rgam ELSE FDUP \ the CC loop BEGIN FSWAP % 1.0 F+ FDUP % -1.0 F< WHILE FSWAP FOVER F* REPEAT rgam F* THEN THEN ; Public: \ 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-integer ABORT" loggam has 0 or negative integer argument " FDUP F0< IF FABS % 1.0 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-integer 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.0 F= OR 0= \ will return x if x is zero or one IF FDUP % 1.0 F< IF rgam-small-x ELSE rgam-large-x THEN THEN ; Reset_Search_Order TEST-CODE? [IF] \ test code ============================================= : gamma-test ( -- ) CR ." gamma( 5.0 ) = " % 5.0 gamma F. ." (should be 24.0) " CR ." gamma( -1.5 ) = " % -1.5 gamma F. ." (should be 2.36327) " CR ." gamma( -0.5 ) = " % -0.5 gamma F. ." (should be -3.54491) " CR ." gamma( 0.5 ) = " % 0.5 gamma F. ." (should be 1.77245) " CR CR ." rgamma( 5.0 ) = " % 5.0 rgamma F. ." (should be 0.0416667) " CR ." rgamma( -1.5 ) = " % -1.5 rgamma F. ." (should be 0.423142) " CR ." rgamma( -0.5 ) = " % -0.5 rgamma F. ." (should be -0.282095) " CR ." rgamma( 0.5 ) = " % 0.5 rgamma F. ." (should be 0.564190) " CR ; [THEN]