Функция, запрограммированная на Фортране 95 для вычисления значений гамма-функции из математики, не выдает правильных значений.
Я пытаюсь реализовать рекурсивную функцию в Fortran 95, которая вычисляет значения гамма-функции, используя приближение Ланцоша (да, я знаю, что для этого есть встроенная функция в стандарте 2003 года и более поздних версиях). Я очень внимательно следил за стандартной формулой, поэтому не уверен, что не так. Правильные значения для гамма-функции имеют решающее значение для некоторых других численных вычислений, которые я делаю, включая численное вычисление полиномов Якоби с помощью соотношения рекурсии.
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
В то время как этот код должен выдавать правильные значения для гамма-функции по всей реальной строке, кажется, что он выдает только постоянное значение, близкое к 122, независимо от ввода. Я подозреваю, что есть какая-то странная проблема с арифметикой с плавающей запятой, которую я не вижу.
pi
? Было бы полезно добавитьimplicit none
. - person Alexander Vogt   schedule 23.12.2018dp
определены в другом месте модуля следующим образом: целое число, параметр :: dp = selected_real_kind (15, 307) real (dp), параметр :: pi = 3,1415926535897 Я отредактирую вопрос чтобы это было ясно. - person pr0gramR   schedule 23.12.2018pi
имеет только одинарную точность. Я предлагаю вам использовать что-то вродеpi = 4 * atan(1._dp)
для его определения. - person evets   schedule 23.12.2018