Plan 9 from Bell Labs’s /usr/web/sources/contrib/fgb/root/sys/src/cmd/4th/lib/gamma.4th

Copyright © 2021 Plan 9 Foundation.
Distributed under the MIT License.
Download the Plan 9 distribution.


\ 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]




Bell Labs OSI certified Powered by Plan 9

(Return to Plan 9 Home Page)

Copyright © 2021 Plan 9 Foundation. All Rights Reserved.
Comments to [email protected].