f_dst_betai_core Function

elemental function f_dst_betai_core(x, a, b) result(betai)

! ==== Description
!! Computes the regularised incomplete beta function. beta_inc and beta_cf
!! algorithms are based on several public domain Fortran and C code,
!! Lentz's algorithm (1976), and modified to use 2008+ intrinsics.

! ==== Declarations
  real(wp), intent(in) :: x      !! upper limit of integral
  real(wp), intent(in) :: a, b   !! shape parameters for beta dist.
  real(wp)             :: betai  !! regularised incomplete beta function
  real(wp)             :: lnbeta !! log of beta function
  real(wp)             :: bt     !! pre-multiplier
  real(wp)             :: cf     !! continued fraction

! ==== Instructions

  ! handle edge cases
  if (x .le. 0.0_wp) then
     betai = 0.0_wp
     return
  elseif (x .ge. 1.0_wp) then
     betai = 1.0_wp
     return
  elseif (a .le. 0.0_wp .or. b .le. 0.0_wp) then
     betai = 0.0_wp  ! alternative: signal an error
     return
  endif

  ! compute log(beta(a, b)) for numerical stability
  ! NOTE: no intrinsic log_beta function; could improve
  lnbeta = log_gamma(a) + log_gamma(b) - log_gamma(a + b)

  ! avoid problems from large a/b
  bt = exp(a * log(x) + b * log(1.0_wp - x) - lnbeta)

  ! switch x with 1-x for num. stability if x > (a+1)/(a+b+2)
  if (x .lt. (a + 1.0_wp) / (a + b + 2.0_wp)) then
    cf = betai_cf_core(x, a, b)
    betai = bt * cf / a
  else
    cf = betai_cf_core(1.0_wp - x, b, a)
    betai = 1.0_wp - bt * cf / b
  endif

  contains

  ! ------------------------------------------------------------------ !
  elemental function betai_cf_core(x, a, b) result(cf)

     ! ==== Description
     !! computes the continued fraction expansion of incomplete beta function.
     !! Based on Lentz's algorithm (1976)

     ! ==== Declarations
     real(wp), intent(in)   :: x                  !! upper limit of integral
     real(wp), intent(in)   :: a, b               !! shape parameters for beta dist.
     real(wp)               :: cf                 !! continued fraction
     real(wp)               :: c, d
     real(wp)               :: aa, del, qab, qam, qap
     real(wp)   , parameter :: eps   = 1.0e-14_wp !! Convergence threshold (how close to 1 the fractional delta must be to stop iterating)
     real(wp)   , parameter :: fpmin = 1.0e-30_wp !! small number to prevent division by zero
     integer(i4), parameter :: i_max = c_bisect_i        !! max. number of iterations
     integer(i4)            :: i

     ! ==== Instructions

     ! starting conditions
     qab = a + b
     qap = a + 1.0_wp
     qam = a - 1.0_wp

     c = 1.0_wp
     d = 1.0_wp - qab * x / qap
     if (abs(d) .lt. fpmin) d = fpmin
     d = 1.0_wp / d
     cf = d

     ! iterative approx.
     do i = 1, i_max
        ! even term
        aa = i * (b - i) * x / ((qam + 2.0_wp * i) * (a + 2.0_wp * i))
        d = 1.0_wp + aa * d
        if (abs(d) .lt. fpmin) d = fpmin
        c = 1.0_wp + aa / c
        if (abs(c) .lt. fpmin) c = fpmin
        d = 1.0_wp / d
        cf = cf * d * c

        ! odd term
        aa = -(a + i) * (qab + i) * x / ((a + 2.0_wp * i) * (qap + 2.0_wp * i))
        d = 1.0_wp + aa * d
        if (abs(d) .lt. fpmin) d = fpmin
        c = 1.0_wp + aa / c
        if (abs(c) .lt. fpmin) c = fpmin
        d = 1.0_wp / d
        del = d * c
        cf = cf * del

        ! check if fractional delta against convergence threshold
        if (abs(del - 1.0_wp) .lt. eps) exit
     enddo

  end function betai_cf_core

end function f_dst_betai_core