Gamma 函数实现未产生正确的值

Question

在 Fortran 95 中编写的用于从数学计算 Gamma 函数值的函数未生成正确的值。

我正在尝试在 Fortran 95 中实现一个递归函数，该函数使用 Lanczos approximation 计算 Gamma 函数的值（是的，我知道在 2003 标准及更高版本中有一个内部函数）。我非常严格地遵循标准公式，所以我不确定哪里出了问题。 Gamma 函数的正确值对于我正在进行的其他一些数值计算至关重要，这些计算涉及通过递归关系对雅可比多项式进行数值计算。

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

虽然这段代码应该在整条实线中为 Gamma 函数生成正确的值，但它似乎只生成一个接近 122 的常数值，而不管输入如何。我怀疑有一些我没有看到的奇怪的浮点运算问题。

Answer 1

您的代码有两个明显的问题

最严重的是，代码在第 42 行越界访问数组，即在循环中
```
do n = 1,9
    t = t + q(n) / (xx + real(n, dp))
end do
```
你有点混淆了你的精度，一些常量是 dp 类型，其他是默认类型

对这些问题进行我认为适当的修复后，您的程序可以正确编译、链接和运行，至少据我所知是这样。见下文：

ian@eris:~/work/Whosebug$ 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/Whosebug$ gfortran -O -std=f2008 -Wall -Wextra -fcheck=all g.f90 
ian@eris:~/work/Whosebug$ ./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/Whosebug$

Gamma 函数实现未产生正确的值

Gamma function implementation not producing correct values

recursion

fortran

gfortran

fortran95

gamma-function