stdlib_specialfunctions_gamma.fypp Source File


Contents


Source Code

#:set WITH_QP = False
#:set WITH_XDP = False
#:include "common.fypp"
#:set CI_KINDS_TYPES = INT_KINDS_TYPES + CMPLX_KINDS_TYPES
module stdlib_specialfunctions_gamma
    use iso_fortran_env, only : qp => real128
    use stdlib_kinds, only :  sp, dp, int8, int16, int32, int64
    use stdlib_error, only : error_stop

    implicit none
    private

    integer(int8), parameter :: max_fact_int8 = 6_int8
    integer(int16), parameter :: max_fact_int16 = 8_int16
    integer(int32), parameter :: max_fact_int32 = 13_int32
    integer(int64), parameter :: max_fact_int64 = 21_int64

    #:for k1, t1 in REAL_KINDS_TYPES
    ${t1}$, parameter :: tol_${k1}$ = epsilon(1.0_${k1}$)
    #:endfor
    real(qp), parameter :: tol_qp = epsilon(1.0_qp)



    public :: gamma, log_gamma, log_factorial
    public :: lower_incomplete_gamma, log_lower_incomplete_gamma
    public :: upper_incomplete_gamma, log_upper_incomplete_gamma
    public :: regularized_gamma_p, regularized_gamma_q



    interface gamma
    !! Gamma function for integer and complex numbers
    !!
        #:for k1, t1 in CI_KINDS_TYPES
        module procedure gamma_${t1[0]}$${k1}$
        #:endfor
    end interface gamma



    interface log_gamma
    !! Logarithm of gamma function
    !!
        #:for k1, t1 in CI_KINDS_TYPES
        module procedure l_gamma_${t1[0]}$${k1}$
        #:endfor
    end interface log_gamma



    interface log_factorial
    !! Logarithm of factorial n!, integer variable
    !!
        #:for k1, t1 in INT_KINDS_TYPES
        module procedure l_factorial_${t1[0]}$${k1}$
        #:endfor
    end interface log_factorial



    interface lower_incomplete_gamma
    !! Lower incomplete gamma function
    !!
        #:for k1, t1 in INT_KINDS_TYPES
          #:for k2, t2 in REAL_KINDS_TYPES
        module procedure ingamma_low_${t1[0]}$${k1}$${k2}$
          #:endfor
        #:endfor

        #:for k1, t1 in REAL_KINDS_TYPES
        module procedure ingamma_low_${t1[0]}$${k1}$
        #:endfor
    end interface lower_incomplete_gamma



    interface log_lower_incomplete_gamma
    !! Logarithm of lower incomplete gamma function
    !!
        #:for k1, t1 in INT_KINDS_TYPES
          #:for k2, t2 in REAL_KINDS_TYPES
        module procedure l_ingamma_low_${t1[0]}$${k1}$${k2}$
          #:endfor
        #:endfor

        #:for k1, t1 in REAL_KINDS_TYPES
        module procedure l_ingamma_low_${t1[0]}$${k1}$
        #:endfor
    end interface log_lower_incomplete_gamma



    interface upper_incomplete_gamma
    !! Upper incomplete gamma function
    !!
        #:for k1, t1 in INT_KINDS_TYPES
          #:for k2, t2 in REAL_KINDS_TYPES
        module procedure ingamma_up_${t1[0]}$${k1}$${k2}$
          #:endfor
        #:endfor

        #:for k1, t1 in REAL_KINDS_TYPES
        module procedure ingamma_up_${t1[0]}$${k1}$
        #:endfor
    end interface upper_incomplete_gamma



    interface log_upper_incomplete_gamma
    !! Logarithm of upper incomplete gamma function
    !!
        #:for k1, t1 in INT_KINDS_TYPES
          #:for k2, t2 in REAL_KINDS_TYPES
        module procedure l_ingamma_up_${t1[0]}$${k1}$${k2}$
          #:endfor
        #:endfor

        #:for k1, t1 in REAL_KINDS_TYPES
        module procedure l_ingamma_up_${t1[0]}$${k1}$
        #:endfor
    end interface log_upper_incomplete_gamma



    interface regularized_gamma_p
    !! Regularized (normalized) lower incomplete gamma function, P
    !!
        #:for k1, t1 in INT_KINDS_TYPES
          #:for k2, t2 in REAL_KINDS_TYPES
        module procedure regamma_p_${t1[0]}$${k1}$${k2}$
          #:endfor
        #:endfor

        #:for k1, t1 in REAL_KINDS_TYPES
        module procedure regamma_p_${t1[0]}$${k1}$
        #:endfor
    end interface regularized_gamma_p



    interface regularized_gamma_q
    !! Regularized (normalized) upper incomplete gamma function, Q
    !!
        #:for k1, t1 in INT_KINDS_TYPES
          #:for k2, t2 in REAL_KINDS_TYPES
        module procedure regamma_q_${t1[0]}$${k1}$${k2}$
          #:endfor
        #:endfor

        #:for k1, t1 in REAL_KINDS_TYPES
        module procedure regamma_q_${t1[0]}$${k1}$
        #:endfor
    end interface regularized_gamma_q



    interface gpx
    ! Incomplete gamma G function.
    ! Internal use only
    !
        #:for k1, t1 in REAL_KINDS_TYPES
        module procedure gpx_${t1[0]}$${k1}$         !for real p and x
        #:endfor

        #:for k1, t1 in INT_KINDS_TYPES
          #:for k2, t2 in REAL_KINDS_TYPES
        module procedure gpx_${t1[0]}$${k1}$${k2}$   !for integer p and real x
          #:endfor
        #:endfor
    end interface gpx



    interface l_gamma
    ! Logarithm of gamma with integer argument for designated output kind.
    ! Internal use only
    !
        #:for k1, t1 in INT_KINDS_TYPES
          #:for k2, t2 in REAL_KINDS_TYPES
          module procedure l_gamma_${t1[0]}$${k1}$${k2}$
          #:endfor
        #:endfor
    end interface l_gamma





contains

    #:for k1, t1 in INT_KINDS_TYPES
    impure elemental function gamma_${t1[0]}$${k1}$(z) result(res)
        ${t1}$, intent(in) :: z
        ${t1}$ :: res, i
        ${t1}$, parameter :: zero = 0_${k1}$, one = 1_${k1}$

        if(z <= zero) call error_stop("Error(gamma): Gamma function argument"  &
            //" must be positive integer.")

        if(z > max_fact_${k1}$) call error_stop("Error(gamma): Gamma function" &
            //" integer argument is greater than the upper limit from which an"&
            //" integer overflow will be generated. Suggest switch to high "   &
            //" precision or convert to real data type")

        res = one

        do i = one, z - one

            res = res * i

        end do

    end function gamma_${t1[0]}$${k1}$

    #:endfor




    #:for k1, t1 in CMPLX_KINDS_TYPES
      #:if k1 == "sp"
          #:set k2 = "dp"
      #:elif k1 == "dp"
          #:set k2 = "qp"
      #:endif
      #:set t2 = "real({})".format(k2)

    impure elemental function gamma_${t1[0]}$${k1}$(z) result(res)
        ${t1}$, intent(in) :: z
        ${t1}$ :: res
        integer :: i

        real(${k1}$), parameter :: zero_k1 = 0.0_${k1}$
        ${t2}$, parameter :: zero = 0.0_${k2}$, half = 0.5_${k2}$,             &
                             one = 1.0_${k2}$, pi = acos(- one), sqpi = sqrt(pi)
        complex(${k2}$) :: y, x, sum

        #:if k1 == "sp"
          #! for single precision input, using double precision for calculation

        integer, parameter :: n = 10
        ${t2}$, parameter :: r = 10.900511_${k2}$
        ${t2}$, parameter :: d(0 : n) = [2.48574089138753566e-5_${k2}$,        &
                                             1.05142378581721974_${k2}$,       &
                                            -3.45687097222016235_${k2}$,       &
                                             4.51227709466894824_${k2}$,       &
                                            -2.98285225323576656_${k2}$,       &
                                             1.05639711577126713_${k2}$,       &
                                         -1.95428773191645870e-1_${k2}$,       &
                                          1.70970543404441224e-2_${k2}$,       &
                                         -5.71926117404305781e-4_${k2}$,       &
                                          4.63399473359905637e-6_${k2}$,       &
                                         -2.71994908488607704e-9_${k2}$]
        ! parameters from above referenced source.

        #:elif k1 == "dp"
          #! for double precision input, using quadruple precision for calculation

        integer, parameter :: n = 24
        ${t2}$, parameter :: r = 25.617904_${k2}$
        ${t2}$, parameter :: d(0 : n)=                                         &
                         [1.0087261714899910504854136977047144166e-11_${k2}$,  &
                              1.6339627701280724777912729825256860624_${k2}$,  &
                          -1.4205787702221583745972794018472259342e+1_${k2}$,  &
                           5.6689501646428786119793943350900908698e+1_${k2}$,  &
                          -1.3766376824252176069406853670529834070e+2_${k2}$,  &
                           2.2739972766608392140035874845640820558e+2_${k2}$,  &
                          -2.7058382145757164380300118233258834430e+2_${k2}$,  &
                          2.39614374587263042692333711131832094166e+2_${k2}$,  &
                          -1.6090450559507517723393498276315290189e+2_${k2}$,  &
                          8.27378183187161305711485619113605553100e+1_${k2}$,  &
                          -3.2678977082742592701862249152153110206e+1_${k2}$,  &
                             9.89018079175824824537131521501652931756_${k2}$,  &
                             -2.2762136356329318377213053650799013041_${k2}$,  &
                          3.93265017303573867227590563182750070164e-1_${k2}$,  &
                          -5.0051054352146209116457193223422284239e-2_${k2}$,  &
                          4.57142601898244576789629257292603538238e-3_${k2}$,  &
                          -2.8922592124650765614787233510990416584e-4_${k2}$,  &
                          1.20833375377219592849746118012697473202e-5_${k2}$,  &
                          -3.1220812187551248389268359432609135033e-7_${k2}$,  &
                          4.55117045361638520378367871355819524460e-9_${k2}$,  &
                         -3.2757632817493581828033170342853173968e-11_${k2}$,  &
                         9.49784279240135747819870224486376897253e-14_${k2}$,  &
                         -7.9480594917454410117072562195702526836e-17_${k2}$,  &
                         1.04692819439870077791406760109955648941e-20_${k2}$,  &
                         -5.8990280044857540075384586350723191533e-26_${k2}$]
        ! parameters from above referenced source.

        #:endif



        if(abs(z % im) < tol_${k1}$) then

            res = cmplx(gamma(z % re), kind = ${k1}$)
            return

        end if

        if(z % re < zero_k1) then

            x = cmplx(abs(z % re), - z % im, kind = ${k1}$)
            y = x - one

        else

            y = z - one

        end if

        sum = cmplx(d(0), kind = ${k2}$)

        do i = 1, n

            sum = sum + d(i) / (y + i)

        end do

        y = exp((y + half) * log(y + half + r) - y) * sum

        y = y * 2 / sqpi                         !Re(z) > 0 return

        if(z % re < zero_k1 ) then

            y = - pi / (sin(pi * x) * x * y)     !Re(z) < 0 return

        end if

        res = y
    end function gamma_${t1[0]}$${k1}$

    #:endfor




    #:for k1, t1 in INT_KINDS_TYPES
    impure elemental function l_gamma_${t1[0]}$${k1}$(z) result(res)
    !
    ! Logarithm of gamma function for integer input
    !
        ${t1}$, intent(in) :: z
        real :: res
        ${t1}$ :: i
        ${t1}$, parameter :: zero = 0_${k1}$, one = 1_${k1}$, two = 2_${k1}$

        if(z <= zero) call error_stop("Error(log_gamma): Gamma function"       &
            //" argument must be positive integer.")

        select case(z)

        case (one)

            res = 0.0

        case (two :)

            res = 0.0

            do i = one, z - one

               res = res + log(real(i))

            end do

        end select
    end function l_gamma_${t1[0]}$${k1}$

    #:endfor




    #:for k1, t1 in INT_KINDS_TYPES
      #:for k2, t2 in REAL_KINDS_TYPES

    impure elemental function l_gamma_${t1[0]}$${k1}$${k2}$(z, x) result(res)
    !
    ! Logarithm of gamma function for integer input with defined precision output
    !
        ${t1}$, intent(in) :: z
        ${t2}$, intent(in) :: x
        ${t2}$ :: res
        ${t1}$ :: i
        ${t1}$, parameter :: zero = 0_${k1}$, one = 1_${k1}$, two = 2_${k1}$
        ${t2}$, parameter :: zero_k2 = 0.0_${k2}$

        if(z <= zero) call error_stop("Error(log_gamma): Gamma function"       &
            //" argument must be positive integer.")

        select case(z)

        case (one)

            res = zero_k2

        case (two :)

            res = zero_k2

            do i = one, z - one

               res = res + log(real(i, ${k2}$))

            end do

        end select
    end function l_gamma_${t1[0]}$${k1}$${k2}$

      #:endfor
    #:endfor




    #:for k1, t1 in CMPLX_KINDS_TYPES
      #:if k1 == "sp"
          #:set k2 = "dp"
      #:elif k1 == "dp"
          #:set k2 = "qp"
      #:endif
      #:set t2 = "real({})".format(k2)
    impure elemental function l_gamma_${t1[0]}$${k1}$(z) result (res)
    !
    ! log_gamma function for any complex number, excluding negative whole number
    ! "Computation of special functions", Shanjie Zhang & Jianmin Jin, 1996, p.48
    ! "Computing the principal branch of log-gamma", D.E.G. Hare,
    ! J. of Algorithms, 25(2), 1997 p. 221–236
    !
    ! Fortran 90 program by Jim-215-Fisher
    !
        ${t1}$, intent(in) :: z
        ${t1}$ :: res, z1, z2
        real(${k1}$) :: d
        integer :: m, i
        complex(${k2}$) :: zr, zr2, sum, s
        real(${k1}$), parameter :: z_limit = 10_${k1}$, zero_k1 = 0.0_${k1}$
        integer, parameter :: n = 20
        ${t2}$, parameter :: zero = 0.0_${k2}$, one = 1.0_${k2}$,              &
                             pi = acos(-one), ln2pi = log(2 * pi)
        ${t2}$, parameter :: a(n) = [                                          &
                           .8333333333333333333333333333333333333333E-1_${k2}$,&
                          -.2777777777777777777777777777777777777778E-2_${k2}$,&
                           .7936507936507936507936507936507936507937E-3_${k2}$,&
                          -.5952380952380952380952380952380952380952E-3_${k2}$,&
                           .8417508417508417508417508417508417508418E-3_${k2}$,&
                          -.1917526917526917526917526917526917526918E-2_${k2}$,&
                           .6410256410256410256410256410256410256410E-2_${k2}$,&
                          -.2955065359477124183006535947712418300654E-1_${k2}$,&
                           .1796443723688305731649384900158893966944E+0_${k2}$,&
                          -.1392432216905901116427432216905901116427E+1_${k2}$,&
                           .1340286404416839199447895100069013112491E+2_${k2}$,&
                          -.1568482846260020173063651324520889738281E+3_${k2}$,&
                           .2193103333333333333333333333333333333333E+4_${k2}$,&
                          -.3610877125372498935717326521924223073648E+5_${k2}$,&
                           .6914722688513130671083952507756734675533E+6_${k2}$,&
                          -.1523822153940741619228336495888678051866E+8_${k2}$,&
                           .3829007513914141414141414141414141414141E+9_${k2}$,&
                         -.1088226603578439108901514916552510537473E+11_${k2}$,&
                          .3473202837650022522522522522522522522523E+12_${k2}$,&
                         -.1236960214226927445425171034927132488108E+14_${k2}$]
        ! parameters from above reference

        z2 = z

        if(z % re < zero_k1) then

            z2 = cmplx(abs(z % re), - z % im, kind = ${k1}$) + 1

        end if

        d = hypot(z2 % re, z2 % im)
        z1 = z2
        m = 0

        if(d <= z_limit) then                       !for small |z|

            m = ceiling(z_limit - d)
            z1 = z2 + m

        end if

        zr = one / z1
        zr2 = zr * zr

        sum = (((a(20) * zr2 + a(19)) * zr2 + a(18)) * zr2 + a(17)) * zr2
        sum = (((sum + a(16)) * zr2 + a(15)) * zr2 + a(14)) * zr2
        sum = (((sum + a(13)) * zr2 + a(12)) * zr2 + a(11)) * zr2
        sum = (((sum + a(10)) * zr2 + a(9)) * zr2 + a(8)) * zr2
        sum = (((sum + a(7)) * zr2 + a(6)) * zr2 + a(5)) * zr2
        sum = (((sum + a(4)) * zr2 + a(3)) * zr2 + a(2)) * zr2
        sum = (sum + a(1)) * zr + ln2pi / 2 - z1 + (z1 - 0.5_${k2}$) * log(z1)

        if(m /= 0) then

            s = cmplx(zero, zero, kind = ${k2}$)

            do i = 1, m

                s = s + log(cmplx(z1, kind = ${k2}$) - i)

            end do

            sum = sum - s

        end if

        if(z % re < zero_k1) then

            sum = log(pi) - log(sin(pi * z)) - sum
            m = ceiling((2 * z % re - 3) / 4)
            sum % im = sum % im + 2 * pi * m * sign(1.0_${k1}$, z % im)

        end if

        res = cmplx(sum, kind = ${k1}$)
    end function l_gamma_${t1[0]}$${k1}$

    #:endfor




    #:for k1, t1 in INT_KINDS_TYPES
    impure elemental function l_factorial_${t1[0]}$${k1}$(n) result(res)
    !
    ! Log(n!)
    !
        ${t1}$, intent(in) :: n
        real :: res
        ${t1}$, parameter :: zero = 0_${k1}$, one = 1_${k1}$, two = 2_${k1}$
        real, parameter :: zero_k2 = 0.0

        if(n < zero) call error_stop("Error(l_factorial): Logarithm of"        &
            //" factorial function argument must be non-negative")

        select case(n)

        case (zero)

            res = zero_k2

        case (one)

            res = zero_k2

        case (two : )

            res = l_gamma(n + 1, 1.0D0)

        end select
    end function l_factorial_${t1[0]}$${k1}$

    #:endfor



    #:for k1, t1 in REAL_KINDS_TYPES
      #:if k1 == "sp"
          #:set k2 = "dp"
      #:elif k1 == "dp"
          #:set k2 = "qp"
      #:endif
      #:set t2 = "real({})".format(k2)

    impure elemental function gpx_${t1[0]}$${k1}$(p, x) result(res)
    !
    ! Approximation of incomplete gamma G function with real argument p.
    !
    ! Based on Rémy Abergel and Lionel Moisan "Algorithm 1006, Fast and
    ! Accurate Evaluation of a Generalized Incomplete Gamma Function", ACM
    ! Transactions on Mathematical Software, March 2020.
    !
    ! Fortran 90 program by Jim-215-Fisher
    !
        ${t1}$, intent(in) :: p, x
        integer :: n, m

        ${t2}$ :: res, p_lim, a, b, g, c, d, y, ss
        ${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}$

        if(p <= zero_k1) call error_stop("Error(gpx): Incomplete gamma"        &
            //" function must have a positive parameter p")

        if(x < -9.0_${k1}$) then

            p_lim = 5.0_${k1}$ * (sqrt(abs(x)) - 1.0_${k1}$)

        elseif(x >= -9.0_${k1}$ .and. x <= zero_k1) then

            p_lim = zero_k1

        else

            p_lim = x

        endif

        if(x < zero_k1 .and. p < p_lim .and. abs(anint(p) - p) > tol_${k1}$)   &
            call error_stop("Error(gpx): Incomplete gamma function with "      &
            //"negative x must come with a whole number p not too small")

        if(p >= p_lim) then     !use modified Lentz method of continued fraction
                                !for eq. (15) in the above reference.
            a = one
            b = p
            g = a / b
            c = a / dm
            d = one / b
            n = 2

            do

                if(mod(n, 2) == 0) then
                    a = (one - p - n / 2) * x
                else
                    a = (n / 2) * x
                end if

                b = p - one + n
                d =  d * a + b

                if(d == zero) d = dm

                c = b + a / c

                if(c == zero) c = dm

                d = one / d
                y = c * d
                g = g * y
                n = n + 1

                if(abs(y - one) < tol_${k2}$) exit

            end do

        else if(x >= zero_k1) then       !use modified Lentz method of continued
                                         !fraction for eq. (16) in the reference.
            a = one
            b = x + one - p
            g = a / b
            c = a / dm
            d = one / b
            n = 2

            do

                a = (n - 1) * (1 + p - n)
                b = b + 2
                d = d * a + b

                if(d == zero) d = dm

                c = b + a / c

                if(c == zero) c = dm

                d = one / d
                y = c * d
                g = g * y
                n = n + 1

                if(abs(y - one) < tol_${k2}$) exit

            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

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

        res = g
    end function gpx_${t1[0]}$${k1}$

    #:endfor



    #:for k1, t1 in INT_KINDS_TYPES
      #:for k2, t2 in REAL_KINDS_TYPES
    impure elemental function gpx_${t1[0]}$${k1}$${k2}$(p, x) result(res)
    !
    ! Approximation of incomplete gamma G function with integer argument p.
    !
    ! Based on Rémy Abergel and Lionel Moisan "Algorithm 1006, Fast and
    ! Accurate Evaluation of a Generalized Incomplete Gamma Function", ACM
    ! Transactions on Mathematical Software, March 2020.
    !
        ${t1}$, intent(in) :: p
        ${t2}$, intent(in) :: x
        ${t2}$ :: res, p_lim, a, b, g, c, d, y
        integer :: n, m
        ${t2}$, parameter :: zero = 0.0_${k2}$, one = 1.0_${k2}$
        ${t2}$, parameter :: dm = tiny(1.0_${k2}$) * 10 ** 6
        ${t1}$, parameter :: zero_k1 = 0_${k1}$, two = 2_${k1}$

        if(p <= zero_k1) call error_stop("Error(gpx): Incomplete gamma "       &
            //"function must have a positive parameter p")

        if(x < -9.0_${k2}$) then

            p_lim = 5.0_${k2}$ * (sqrt(abs(x)) - 1.0_${k2}$)

        else if(x >= -9.0_${k2}$ .and. x <= zero) then

            p_lim = zero

        else

            p_lim = x

        end if

        if(real(p, ${k2}$) >= p_lim) then

            a = one
            b = p
            g = a / b
            c = a / dm
            d = one / b
            n = 2

            do

                if(mod(n, 2) == 0) then

                    a = (1 - p - n / 2) * x

                else

                    a = (n / 2) * x

                end if

                b = p - 1 + n
                d =  d * a + b

                if(d == zero) d = dm

                c = b + a / c

                if(c == zero) c = dm

                d = one / d
                y = c * d
                g = g * y
                n = n + 1

                if(abs(y - one) < tol_${k2}$) exit

            end do

        else if(x >= zero) then

            a = one
            b = x + 1 - p
            g = a / b
            c = a / dm
            d = one / b
            n = 2

            do

                a = -(n - 1) * (n - 1 - p)
                b = b + 2
                d = d * a + b

                if(d == zero) d = dm

                c = b + a / c

                if(c == zero) c = dm

                d = one / d
                y = c * d
                g = g * y
                n = n + 1

                if(abs(y - one) < tol_${k2}$) exit

            end do

        else

            a = -x
            c = one / a
            d = p - 1
            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(int(n, ${k1}$) > (p - two) / two .or. y < b * tol_${k2}$) exit

            end do

            if(y >= b * tol_${k2}$ .and. mod(p, two) /= zero_k1)               &
                b = b + d * c / a

            g = ((-1) ** p * exp(-a + l_gamma(p, one) - (p - 1) * log(a))      &
                + b ) / a

        end if

        res = g
    end function gpx_${t1[0]}$${k1}$${k2}$

      #:endfor
    #:endfor



    #:for k1, t1 in REAL_KINDS_TYPES
    impure elemental function ingamma_low_${t1[0]}$${k1}$(p, x) result(res)
    !
    ! Approximation of lower incomplete gamma function with real p.
    !
        ${t1}$, intent(in) :: p, x
        ${t1}$ :: res, s1, y
        ${t1}$, parameter :: zero = 0.0_${k1}$, one = 1.0_${k1}$

        if(x == zero) then

            res = zero

        else if(x > p) then

            s1 = log_gamma(p)
            y = one - exp(-x + p * log(x) - s1) * gpx(p, x)
            res = exp(s1 + log(y))

        else if(x <= p .and. x > zero) then

            s1 = -x + p * log(x)
            res = gpx(p, x) * exp(s1)

        else

            call error_stop("Error(Logarithm of upper incomplete gamma "       &
              //"function): negative x must be with integer p")

        end if
    end function ingamma_low_${t1[0]}$${k1}$

    #:endfor



    #:for k1, t1 in INT_KINDS_TYPES
      #:for k2, t2 in REAL_KINDS_TYPES
    impure elemental function ingamma_low_${t1[0]}$${k1}$${k2}$(p, x)          &
        result(res)
    !
    ! Approximation of lower incomplete gamma function with integer p.
    !
        ${t1}$, intent(in) :: p
        ${t2}$, intent(in) :: x
        ${t2}$ :: res, s1, y
        ${t2}$, parameter :: zero = 0.0_${k2}$, one = 1.0_${k2}$

        if(x == zero) then

            res = zero

        else if(x > real(p, ${k2}$)) then

            s1 = l_gamma(p, one)
            y = one - exp(-x + p * log(x) - s1) * gpx(p, x)
            res = exp(s1 + log(y))

        else if(x <= real(p, ${k2}$) .and. x > zero) then

            s1 = -x + p * log(x)
            res = gpx(p, x) * exp(s1)

        else

            s1 = -x + p * log(abs(x))
            res = gpx(p, x) * exp(s1)
            res = (-1) ** p * res

        end if
    end function ingamma_low_${t1[0]}$${k1}$${k2}$

      #:endfor
    #:endfor



    #:for k1, t1 in REAL_KINDS_TYPES
    impure elemental function l_ingamma_low_${t1[0]}$${k1}$(p, x) result(res)

        ${t1}$, intent(in) :: p, x
        ${t1}$ :: res, s1, y
        ${t1}$, parameter :: zero = 0.0_${k1}$, one = 1.0_${k1}$


        if(x == zero) then

            res = zero

        else if(x > p) then

            s1 = log_gamma(p)
            y = one - exp(-x + p * log(x) - s1) * gpx(p, x)
            res = s1 + log(y)

        else if(x <= p .and. x > zero) then

            s1 = -x + p * log(abs(x))
            res = log(abs(gpx(p, x))) + s1

        else

            call error_stop("Error(Logarithm of upper incomplete gamma "       &
              //"function): negative x must be with integer p")

        end if
    end function l_ingamma_low_${t1[0]}$${k1}$

    #:endfor




    #:for k1, t1 in INT_KINDS_TYPES
      #:for k2, t2 in REAL_KINDS_TYPES
    impure elemental function l_ingamma_low_${t1[0]}$${k1}$${k2}$(p, x)        &
        result(res)

        ${t1}$, intent(in) :: p
        ${t2}$, intent(in) :: x
        ${t2}$ :: res, s1, y
        ${t2}$, parameter :: zero = 0.0_${k2}$, one = 1.0_${k2}$

        if(x == zero) then

            res = zero

        else if(x > real(p, ${k2}$)) then

            s1 = l_gamma(p, one)
            y = one - exp(-x + p * log(x) - s1) * gpx(p, x)
            res = s1 + log(y)

        else if(x <= real(p, ${k2}$)) then

            s1 = -x + p * log(abs(x))
            res = log(abs(gpx(p, x))) + s1

        end if
    end function l_ingamma_low_${t1[0]}$${k1}$${k2}$

      #:endfor
    #:endfor



    #:for k1, t1 in REAL_KINDS_TYPES
    impure elemental function ingamma_up_${t1[0]}$${k1}$(p, x) result(res)
    !
    ! Approximation of upper incomplete gamma function with real p.
    !
        ${t1}$, intent(in) :: p, x
        ${t1}$ :: res, s1, y
        ${t1}$, parameter :: zero = 0.0_${k1}$, one = 1.0_${k1}$

        if(x == zero) then

            res = gamma(p)

        else if(x > p) then

            s1 = -x + p * log(x)
            res = gpx(p, x) * exp(s1)

        else if(x <= p .and. x > zero) then

            y = log_gamma(p)
            s1 = -x + p * log(x) - y
            res = (one - gpx(p, x) * exp(s1)) * exp(y)

        else


            call error_stop("Error(Logarithm of upper incomplete gamma "       &
              //"function): negative x must be with integer p")

        end if
    end function ingamma_up_${t1[0]}$${k1}$

    #:endfor



    #:for k1, t1 in INT_KINDS_TYPES
      #:for k2, t2 in REAL_KINDS_TYPES
    impure elemental function ingamma_up_${t1[0]}$${k1}$${k2}$(p, x)           &
        result(res)
    !
    ! Approximation of upper incomplete gamma function with integer p.
    !
        ${t1}$, intent(in) :: p
        ${t2}$, intent(in) :: x
        ${t2}$ :: res, s1, y
        ${t2}$, parameter :: zero = 0.0_${k2}$, one = 1.0_${k2}$

        if(x == zero) then

            res = gamma(real(p, ${k2}$))

        else if(x > real(p, ${k2}$)) then

            s1 = -x + p * log(x)
            res = gpx(p, x) * exp(s1)

        else if(x <= real(p, ${k2}$) .and. x > zero) then

            y = l_gamma(p, one)
            s1 = -x + p * log(x) - y
            res = gpx(p, x) * exp(s1)
            res = (one - res) * exp(y)

        else

            y = l_gamma(p, one)
            s1 = -x + p * log(abs(x)) - y
            res = gpx(p, x) * exp(s1)
            res = (one - (-1) ** p * res) * exp(y)

        end if
    end function ingamma_up_${t1[0]}$${k1}$${k2}$

      #:endfor
    #:endfor



    #:for k1, t1 in REAL_KINDS_TYPES
    impure elemental function l_ingamma_up_${t1[0]}$${k1}$(p, x) result(res)

        ${t1}$, intent(in) :: p, x
        ${t1}$ :: res, s1, y
        ${t1}$, parameter :: zero = 0.0_${k1}$, one = 1.0_${k1}$


        if(x == zero) then

            res = log_gamma(p)

        else if(x > p) then

            s1 = -x + p * log(x)
            res = log(gpx(p, x)) + s1

        else if(x <= p .and. x > zero) then

            y= log_gamma(p)
            s1 = -x + p * log(x) - y
            res = gpx(p, x) * exp(s1)
            res = log(one - res) + y

        else

            call error_stop("Error(Logarithm of upper incomplete gamma "       &
              //"function): negative x must be with integer p")

        end if
    end function l_ingamma_up_${t1[0]}$${k1}$

    #:endfor




    #:for k1, t1 in INT_KINDS_TYPES
      #:for k2, t2 in REAL_KINDS_TYPES
    impure elemental function l_ingamma_up_${t1[0]}$${k1}$${k2}$(p, x)         &
        result(res)

        ${t1}$, intent(in) :: p
        ${t2}$, intent(in) :: x
        ${t2}$ :: res, s1, y
        ${t2}$, parameter :: zero = 0.0_${k2}$, one = 1.0_${k2}$

        if(x == zero) then

            res = l_gamma(p, one)

        else if(x > real(p, ${k2}$)) then

            s1 = -x + p * log(x)
            res = log(gpx(p, x)) + s1

        else if(x <= real(p, ${k2}$) .and. x > zero) then

            y = l_gamma(p, one)
            s1 = -x + p * log(x) - y
            res = gpx(p, x) * exp(s1)
            res = log(one - res) + y

        else

            y = l_gamma(p, one)
            s1 = -x + p * log(abs(x)) + log(gpx(p, x))
            res = (-1) ** p * exp(s1)
            res = log(abs(exp(y) - res))

        end if
    end function l_ingamma_up_${t1[0]}$${k1}$${k2}$

      #:endfor
    #:endfor




    #:for k1, t1 in REAL_KINDS_TYPES
    impure elemental function regamma_p_${t1[0]}$${k1}$(p, x) result(res)
    !
    ! Approximation of regularized incomplete gamma function P(p,x) for real p
    !
        ${t1}$, intent(in) :: p, x
        ${t1}$ :: res, s1
        ${t1}$, parameter :: zero = 0.0_${k1}$, one = 1.0_${k1}$

        if(x < zero) call error_stop("Error(regamma_p): Regularized gamma_p"   &
            //" function is not defined at x < 0")


        if(x == zero) then

            res = zero

        else if(x > p) then

            s1 = -x + p * log(x) - log_gamma(p)
            res = one - exp(s1 + log(gpx(p,x)))

        else if(x <= p) then

            s1 = -x + p * log(abs(x)) - log_gamma(p)
            res = exp(log(gpx(p, x)) + s1)

        end if
    end function regamma_p_${t1[0]}$${k1}$

    #:endfor



    #:for k1, t1 in INT_KINDS_TYPES
      #:for k2, t2 in REAL_KINDS_TYPES
    impure elemental function regamma_p_${t1[0]}$${k1}$${k2}$(p, x) result(res)
    !
    ! Approximation of regularized incomplete gamma function P(p,x) for integer p
    !
        ${t1}$, intent(in) :: p
        ${t2}$, intent(in) :: x
        ${t2}$ :: res, s1
        ${t2}$, parameter :: zero = 0.0_${k2}$, one = 1.0_${k2}$

        if(x < zero) call error_stop("Error(regamma_p): Regularized gamma_p"   &
            //" function is not defined at x < 0")


        if(x == zero) then

            res = zero

        else if(x > real(p, ${k2}$)) then

            s1 = -x + p * log(x) - l_gamma(p, one)
            res = one - exp(s1 + log(gpx(p,x)))

        else if(x <= real(p, ${k2}$)) then

            s1 = -x + p * log(abs(x)) - l_gamma(p, one)
            res = exp(log(gpx(p, x)) + s1)

        end if
    end function regamma_p_${t1[0]}$${k1}$${k2}$

      #:endfor
    #:endfor



    #:for k1, t1 in REAL_KINDS_TYPES
    impure elemental function regamma_q_${t1[0]}$${k1}$(p, x) result(res)
    !
    ! Approximation of regularized incomplete gamma function Q(p,x) for real p
    !
        ${t1}$, intent(in) :: p, x
        ${t1}$ :: res, s1
        ${t1}$, parameter :: zero = 0.0_${k1}$, one = 1.0_${k1}$

        if(x < zero) call error_stop("Error(regamma_p): Regularized gamma_q"   &
            //" function is not defined at x < 0")


        if(x == zero) then

            res = one

        else if(x > p) then

            s1 = -x + p * log(x) - log_gamma(p)
            res = exp(s1 + log(gpx(p,x)))

        else if(x <= p) then

            s1 = -x + p * log(abs(x)) - log_gamma(p)
            res = one - exp(log(gpx(p, x)) + s1)

        end if
    end function regamma_q_${t1[0]}$${k1}$

    #:endfor



    #:for k1, t1 in INT_KINDS_TYPES
      #:for k2, t2 in REAL_KINDS_TYPES
    impure elemental function regamma_q_${t1[0]}$${k1}$${k2}$(p, x) result(res)
    !
    ! Approximation of regularized incomplet gamma function Q(p,x) for integer p
    !
        ${t1}$, intent(in) :: p
        ${t2}$, intent(in) :: x
        ${t2}$ :: res, s1
        ${t2}$, parameter :: zero = 0.0_${k2}$, one = 1.0_${k2}$

        if(x < zero) call error_stop("Error(regamma_q): Regularized gamma_q"   &
            //" function is not defined at x < 0")


        if(x == zero) then

            res = one

        else if(x > real(p, ${k2}$)) then

            s1 = -x + p * log(x) - l_gamma(p, one)
            res = exp(log(gpx(p,x)) + s1)

        elseif(x <= real(p, ${k2}$)) then

            s1 = -x + p * log(abs(x)) - l_gamma(p, one)
            res = one - exp(s1 + log(gpx(p,x)))

        end if
     end function regamma_q_${t1[0]}$${k1}$${k2}$

      #:endfor
    #:endfor

end module stdlib_specialfunctions_gamma