1

在 Fortran 95 中编程以从数学中计算 Gamma 函数值的函数不会产生正确的值。

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

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 的常数值。我怀疑有一些我没有看到的奇怪的浮点算术问题。

4

1 回答 1

3

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

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

我认为对这些程序进行适当的修复可以正确编译、链接和运行，至少就我所见而言。见下文：

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$
于 2018-12-24T09:56:27.827 回答