Gamma function implementation not producing correct values

Question

Function programmed in Fortran 95 to compute values of the Gamma function from mathematics is not producing the correct values.

I am trying to implement a recursive function in Fortran 95 that computes values of the Gamma function using the Lanczos approximation (yes I know that there is an intrinsic function for this in the 2003 standard and later). I've followed the standard formula very closely so I'm not certain what is wrong. Correct values for the Gamma function are crucial for some other numerical computations I am doing involving the numerical computation of the Jacobi polynomials by means of a recursion relation.

program testGam

implicit none

integer, parameter      :: dp = selected_real_kind(15,307)
real(dp), parameter     :: pi = 3.14159265358979324

real(dp), dimension(10) :: xGam, Gam
integer                 :: n

xGam = (/ -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5 /)
do n = 1,10
    Gam(n) = GammaFun(xGam(n))
end do

do n = 1,10
    write(*,*) xGam(n), Gam(n)
end do


contains    

    recursive function GammaFun(x) result(G)

    real(dp), intent(in) :: x
    real(dp)             :: G
    real(dp), dimension(0:8), parameter :: q = &
           (/ 0.99999999999980993, 676.5203681218851, -1259.1392167224028, &
              771.32342877765313, -176.61502916214059, 12.507343278686905, &
              -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7 /)
    real(dp)             :: t, w, xx
    integer              :: n

    xx = x

    if ( xx < 0.5_dp ) then
        G = pi / ( sin(pi*xx)*GammaFun(1.0_dp - xx) )
    else
        xx = xx - 1.0_dp
        t  = q(0)
        do n = 1,9
            t = t + q(n) / (xx + real(n, dp))
        end do
        w = xx + 7.5_dp
        G = sqrt(2.0_dp*pi)*(w**(xx + 0.5_dp))*exp(-w)*t
    end if

    end function GammaFun

end program testGam

Whereas this code should be producing correct values for the Gamma function over the whole real line, it seems only to produce a constant value close to 122 regardless of the input. I suspect that there is some weird floating point arithmetic issue that I am not seeing.

Where is `pi` defined? It would help to add an `implicit none`. — Alexander Vogt, Dec 23 '18 at 17:49
Oh right, my apologies. The constant pi as well as the double precision modifier `dp` are defined elsewhere in the module as follows: integer, parameter :: dp = selected_real_kind(15, 307) real(dp), parameter :: pi = 3.1415926535897 I'll edit the question to make this clear. — pr0gramR, Dec 23 '18 at 17:55
Your constant `pi` has only single precision. I suggest that you use something like `pi = 4 * atan(1._dp)` to define it. — evets, Dec 23 '18 at 18:17
BTW, I forgot to mention that computing the true gamma function is much harder than you may anticipate. I suggest you search Netlib or look at the C language implementation in FreeBSD's math library. — evets, Dec 23 '18 at 18:29
Yes I know there are other canned options for approximating the Gamma function, I was just curious as to why this implementation was not producing correct values which would hopefully catch some programming error that I'm not catching that could cause errors if reproduced in my other projects. — pr0gramR, Dec 23 '18 at 20:22
You've missed the point about the "canned options" and my comment **computing the true gamma function is much harder than you may anticipate**. If you have a poor algorithm, no coding will fix that. — evets, Dec 24 '18 at 17:31

score 3 · Accepted Answer · answered Dec 24 '18 at 09:56

There are two obvious issues with your code

Most seriously the code accesses an array out of bounds at line 42, i.e. in the loop
```
do n = 1,9
    t = t + q(n) / (xx + real(n, dp))
end do
```
You have mixed up your precision somewhat, with some of the constants being of kind dp, other being of default kind

Making what I believe are the appropriate fixes to these your program compiles, links and runs correctly, at least as far as I can see. See below:

ian@eris:~/work/stackoverflow$ cat g.f90
program testGam

implicit none

integer, parameter      :: dp = selected_real_kind(15,307)
real(dp), parameter     :: pi = 3.14159265358979324_dp

real(dp), dimension(10) :: xGam, Gam
integer                 :: n

xGam = (/ -3.5_dp, -2.5_dp, -1.5_dp, -0.5_dp, 0.5_dp, 1.5_dp, 2.5_dp, 3.5_dp, 4.5_dp, 5.5_dp /)
do n = 1,10
    Gam(n) = GammaFun(xGam(n))
end do

do n = 1,10
    write(*,*) xGam(n), Gam(n), gamma( xGam( n ) ), Abs( Gam( n ) - gamma( xGam( n ) ) )
end do


contains    

    recursive function GammaFun(x) result(G)

    real(dp), intent(in) :: x
    real(dp)             :: G
    real(dp), dimension(0:8), parameter :: q = &
           (/ 0.99999999999980993_dp, 676.5203681218851_dp, -1259.1392167224028_dp, &
              771.32342877765313_dp, -176.61502916214059_dp, 12.507343278686905_dp, &
              -0.13857109526572012_dp, 9.9843695780195716e-6_dp, 1.5056327351493116e-7_dp /)
    real(dp)             :: t, w, xx
    integer              :: n

    xx = x

    if ( xx < 0.5_dp ) then
        G = pi / ( sin(pi*xx)*GammaFun(1.0_dp - xx) )
    else
        xx = xx - 1.0_dp
        t  = q(0)
        do n = 1,8
            t = t + q(n) / (xx + real(n, dp))
        end do
        w = xx + 7.5_dp
        G = sqrt(2.0_dp*pi)*(w**(xx + 0.5_dp))*exp(-w)*t
    end if

    end function GammaFun

end program testGam

ian@eris:~/work/stackoverflow$ gfortran -O -std=f2008 -Wall -Wextra -fcheck=all g.f90 
ian@eris:~/work/stackoverflow$ ./a.out
  -3.5000000000000000       0.27008820585226917       0.27008820585226906        1.1102230246251565E-016
  -2.5000000000000000      -0.94530872048294168      -0.94530872048294179        1.1102230246251565E-016
  -1.5000000000000000        2.3632718012073521        2.3632718012073548        2.6645352591003757E-015
 -0.50000000000000000       -3.5449077018110295       -3.5449077018110318        2.2204460492503131E-015
  0.50000000000000000        1.7724538509055159        1.7724538509055161        2.2204460492503131E-016
   1.5000000000000000       0.88622692545275861       0.88622692545275805        5.5511151231257827E-016
   2.5000000000000000        1.3293403881791384        1.3293403881791370        1.3322676295501878E-015
   3.5000000000000000        3.3233509704478430        3.3233509704478426        4.4408920985006262E-016
   4.5000000000000000        11.631728396567446        11.631728396567450        3.5527136788005009E-015
   5.5000000000000000        52.342777784553583        52.342777784553519        6.3948846218409017E-014
ian@eris:~/work/stackoverflow$

Gamma function implementation not producing correct values

1 Answers1