44#:set CI_KINDS_TYPES = INT_KINDS_TYPES + C_KINDS_TYPES
55module stdlib_specialfunctions_gamma
66 use iso_fortran_env, only : qp => real128
7+ use ieee_arithmetic, only: ieee_value, ieee_quiet_nan
78 use stdlib_kinds, only : sp, dp, int8, int16, int32, int64
89 use stdlib_error, only : error_stop
910
@@ -575,9 +576,9 @@ contains
575576 ! Fortran 90 program by Jim-215-Fisher
576577 !
577578 ${t1},ドル intent(in) :: p, x
578- integer :: n, m
579+ integer :: n
579580
580- ${t2}$ :: res, p_lim, a, b, g, c, d, y, ss
581+ ${t2}$ :: res, p_lim, a, b, g, c, d, y
581582 ${t2},ドル parameter :: zero = 0.0_${k2},ドル one = 1.0_${k2}$
582583 ${t2},ドル parameter :: dm = tiny(1.0_${k2}$) * 10 ** 6
583584 ${t1},ドル parameter :: zero_k1 = 0.0_${k1}$
@@ -603,6 +604,9 @@ contains
603604 call error_stop("Error(gpx): Incomplete gamma function with " &
604605 //"negative x must come with a whole number p not too small")
605606
607+ if(x < zero_k1) call error_stop("Error(gpx): Incomplete gamma" &
608+ // " function with negative x must have an integer parameter p")
609+ 606610 if(p >= p_lim) then !use modified Lentz method of continued fraction
607611 !for eq. (15) in the above reference.
608612 a = one
@@ -668,30 +672,9 @@ contains
668672
669673 end do
670674
671- else !Algorithm 2 in the reference
672- 673- m = nint(ss)
674- a = - x
675- c = one / a
676- d = p - one
677- b = c * (a - d)
678- n = 1
679- 680- do
681- 682- c = d * (d - one) / (a * a)
683- d = d - 2
684- y = c * (a - d)
685- b = b + y
686- n = n + 1
687- 688- if(n > int((p - 2) / 2) .or. y < b * tol_${k2}$) exit
689- 690- end do
691- 692- if(y >= b * tol_${k2}$ .and. mod(m , 2) /= 0) b = b + d * c / a
675+ else
676+ g = ieee_value(1._${k1},ドル ieee_quiet_nan)
693677
694- g = ((-1) ** m * exp(-a + log_gamma(p) - (p - 1) * log(a)) + b) / a
695678 end if
696679
697680 res = g
0 commit comments