f_dst_gamma_ppf_core Function

elemental function f_dst_gamma_ppf_core(p, alpha, beta, loc) result(x)

! ==== Description
!! Percent point function/quantile function for gamma distribution.

! ==== Declarations
  real(wp)   , intent(in) :: p                  !! probability between 0.0 - 1.0
  real(wp)   , intent(in) :: alpha              !! shape  parameter
  real(wp)   , intent(in) :: beta               !! scale parameter
  real(wp)   , intent(in) :: loc                !! location parameter
  integer(i4), parameter  :: i_max = c_bisect_i !! max. number of iterations
  real(wp)   , parameter  :: tol = c_bisect_tol !! p deviation tolerance
  real(wp)                :: a, b               !! section bounds for bisection algorithm
  real(wp)                :: x_mid, p_mid       !! x and p mid points in bisection algorithm
  integer(i4)             :: i                  !! for iteration
  real(wp)                :: x                  !! sample position

! ==== Instructions

! ---- compute PPF

  ! set initial section
  a = loc
  b = loc + alpha * beta * 10.0_wp

  ! iteratively refine with bisection method
  do i = 1, i_max
     x_mid = 0.5_wp * (a + b)
     ! difference between passed p and new mid point p
     p_mid = f_dst_gamma_cdf_core(x_mid, alpha=alpha&
          &, beta=beta, loc=loc, tail="left") - p
     ! check if difference is acceptable, update section if not
     if (abs(p_mid) .lt. tol) then
        ! pass final x value
        x = x_mid
        return
     elseif (p_mid .lt. 0.0_wp) then
        a = x_mid
     else
        b = x_mid
     endif
  enddo

  ! if x not found in iterations, pass sentinel
  if (i .ge. i_max) x = c_sentinel_r

end function f_dst_gamma_ppf_core