Skip to content

fix: gamma functions #943

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 3 commits into from
Mar 7, 2025
Merged

fix: gamma functions #943

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
069501e
fix: special gamma functions
wassup05 Mar 6, 2025
c0ce8dc
remove unreachable code
wassup05 Mar 6, 2025
499052f
add else clause to cover all branches
wassup05 Mar 6, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
33 changes: 8 additions & 25 deletions src/stdlib_specialfunctions_gamma.fypp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
#:set CI_KINDS_TYPES = INT_KINDS_TYPES + C_KINDS_TYPES
module stdlib_specialfunctions_gamma
use iso_fortran_env, only : qp => real128
use ieee_arithmetic, only: ieee_value, ieee_quiet_nan
use stdlib_kinds, only : sp, dp, int8, int16, int32, int64
use stdlib_error, only : error_stop

Expand Down Expand Up @@ -575,9 +576,9 @@ contains
! Fortran 90 program by Jim-215-Fisher
!
${t1}$, intent(in) :: p, x
integer :: n, m
integer :: n

${t2}$ :: res, p_lim, a, b, g, c, d, y, ss
${t2}$ :: res, p_lim, a, b, g, c, d, y
${t2}$, parameter :: zero = 0.0_${k2}$, one = 1.0_${k2}$
${t2}$, parameter :: dm = tiny(1.0_${k2}$) * 10 ** 6
${t1}$, parameter :: zero_k1 = 0.0_${k1}$
Expand All @@ -603,6 +604,9 @@ contains
call error_stop("Error(gpx): Incomplete gamma function with " &
//"negative x must come with a whole number p not too small")

if(x < zero_k1) call error_stop("Error(gpx): Incomplete gamma" &
// " function with negative x must have an integer parameter p")

if(p >= p_lim) then !use modified Lentz method of continued fraction
!for eq. (15) in the above reference.
a = one
Expand Down Expand Up @@ -668,30 +672,9 @@ contains

end do

else !Algorithm 2 in the reference

m = nint(ss)
a = - x
c = one / a
d = p - one
b = c * (a - d)
n = 1

do

c = d * (d - one) / (a * a)
d = d - 2
y = c * (a - d)
b = b + y
n = n + 1

if(n > int((p - 2) / 2) .or. y < b * tol_${k2}$) exit

end do

if(y >= b * tol_${k2}$ .and. mod(m , 2) /= 0) b = b + d * c / a
else
g = ieee_value(1._${k1}$, ieee_quiet_nan)

g = ((-1) ** m * exp(-a + log_gamma(p) - (p - 1) * log(a)) + b) / a
end if

res = g
Expand Down
Loading