Реализация гамма-функции не дает правильных значений

Question

Функция, запрограммированная в Fortran 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, независимо от входных данных. Я подозреваю, что есть какая-то странная арифметическая проблема с плавающей точкой, которую я не вижу.

Answer 1

Есть два очевидных проблемы с вашим кодом

1. Наиболее серьезно код обращается к массиву вне границ в строке 42, то есть в цикле

   do n = 1,9
       t = t + q(n) / (xx + real(n, dp))
   end do
   

2. Вы немного перепутали свою точность, причем некоторые константы имеют тип dp, а другие имеют тип по умолчанию

Создание того, что я считаю, является подходящими исправлениями для вашей программы, которая компилирует, связывает и запускает правильно, по крайней мере, насколько я могу видеть. Увидеть ниже:

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