f_dst_f_ppf_core Function

elemental function f_dst_f_ppf_core(p, d1, d2, loc, scale) result(x)

! ==== Description
!! Percent point function / quantile function for the F distribution.

! ==== Declarations
  real(wp), intent(in)   :: p                  !! probability (0.0 < p < 1.0)
  real(wp), intent(in)   :: d1                 !! numerator degrees of freedom
  real(wp), intent(in)   :: d2                 !! denominator degrees of freedom
  real(wp), intent(in)   :: loc                !! location parameter
  real(wp), intent(in)   :: scale              !! scale parameter
  integer(i4), parameter :: i_max = c_bisect_i !! max. number of iterations
  real(wp)   , parameter :: tol = c_bisect_tol !! tolerance for convergence
  real(wp)               :: a, b               !! search bounds
  real(wp)               :: x_mid, p_mid       !! midpoint and its CDF
  integer(i4)            :: i                  !! iteration counter
  real(wp)               :: x                  !! result: quantile at p

! ==== Instructions

! ---- compute PPF

  ! set initial section
  a = loc
  b = loc + scale * 100.0_wp  ! heuristically large upper bound

  ! iteratively refine with bisection method
  do i = 1, i_max
     x_mid = 0.5_wp * (a + b)
     p_mid = f_dst_f_cdf_core(x_mid, d1, d2&
          &, loc=loc, scale=scale, tail="left") - p
     if (abs(p_mid) .lt. tol) then
        x = x_mid
        return
     else if (p_mid .lt. 0.0_wp) then
        a = x_mid
     else
        b = x_mid
     end if
  end do

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

end function f_dst_f_ppf_core