【发布时间】：2018-12-23 17:44:36
【问题描述】：
在 Fortran 95 中编程以从数学中计算 Gamma 函数值的函数没有产生正确的值。
我正在尝试在 Fortran 95 中实现一个递归函数，该函数使用 Lanczos approximation 计算 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 的常数值。我怀疑有一些我没有看到的奇怪的浮点算术问题。
【问题讨论】：
-
pi在哪里定义？添加implicit none会有所帮助。 -
哦，对了，我很抱歉。常量 pi 和双精度修饰符
dp在模块的其他地方定义如下：整数，参数 :: dp = selected_real_kind(15, 307) real(dp)，参数 :: pi = 3.1415926535897 我将编辑这个问题要弄清楚。 -
你的常量
pi只有单精度。我建议你使用类似pi = 4 * atan(1._dp)的东西来定义它。 -
顺便说一句，我忘了说计算真正的伽马函数比你想象的要困难得多。建议你搜索 Netlib 或者查看 FreeBSD 的数学库中的 C 语言实现。
-
是的，我知道还有其他用于逼近 Gamma 函数的预设选项，我只是好奇为什么这个实现没有产生正确的值，希望能捕捉到一些我没有捕捉到的编程错误如果在我的其他项目中复制可能会导致错误。
标签： recursion fortran gfortran fortran95 gamma-function