Gamma function implementation not producing correct values

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.

There are two obvious issues with your code

1. 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
   

2. 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$