dppt Function

   real(wp) pure function dppt(p, idf) result(dpptr)
   !! Compute the percent point function value for the student's T distribution with `idf`
   !! degrees of freedom.
      ! Adapted from DATAPAC subroutine TPPF, with modifications to facilitate conversion to
      ! real(wp) automatically.
      ! Routines Called  DPPNML
      ! Date Written   901207   (YYMMDD)
      ! Revision Date  920304   (YYMMDD)
      !***Author  Filliben, James J.,
      !       Statistical Engineering Division
      !       National Bureau of Standards
      !       Washington, D. C. 20234
      !       (Original Version--October   1975.)
      !       (Updated         --November  1975.)
      !***Description
      !       --For IDF = 1 AND IDF = 2, the percent point function
      !       for the T distribution exists in simple closed form
      !       and so the computed percent points are exact.
      !       --For IDF between 3 and 6, inclusively, the approximation
      !       is augmented by 3 iterations of Newton's method to
      !       improve the accuracy, especially for P near 0 or 1.
      !***References  National Bureau of Standards Applied Mathmatics
      !       Series 55, 1964, Page 949, Formula 26.7.5.
      !       Johnson and Kotz, Continuous Univariate Distributions,
      !       Volume 2, 1970, Page 102, Formula 11.
      !       Federighi, "Extended Tables of the Percentage Points
      !       of Student"S T Distribution, Journal of the American
      !       Statistical Association, 1969, Pages 683-688.
      !       Hastings and Peacock, Statistical Distributions, A
      !       Handbook for Students and Practitioners, 1975,
      !       Pages 120-123.

      use odrpack_kinds, only: pi, zero, half, one, two, three, eight, fiftn

      real(wp), intent(in) :: p
         !! The probability at which the percent point is to be evaluated. `p` must lie between
         !! 0.0 and 1.0, exclusive.
      integer, intent(in) :: idf
         !! The (positive integer) degrees of freedom.

      ! Local scalars
      real(wp), parameter :: b21 = 4.0E0_wp, &
                             b31 = 96.0E0_wp, &
                             b32 = 5.0E0_wp, &
                             b33 = 16.0E0_wp, &
                             b34 = 3.0E0_wp, &
                             b41 = 384.0E0_wp, &
                             b42 = 3.0E0_wp, &
                             b43 = 19.0E0_wp, &
                             b44 = 17.0E0_wp, &
                             b45 = -15.0E0_wp, &
                             b51 = 9216.0E0_wp, &
                             b52 = 79.0E0_wp, &
                             b53 = 776.0E0_wp, &
                             b54 = 1482.0E0_wp, &
                             b55 = -1920.0E0_wp, &
                             b56 = -945.0E0_wp
      real(wp) :: arg, c, con, d1, d3, d5, d7, d9, df, ppfn, s, term1, term2, term3, &
                  term4, term5, z
      integer :: ipass, maxit

      ! Variable definitions (alphabetically)
      !  ARG:    A value used in the approximation.
      !  B21:    A parameter used in the approximation.
      !  B31:    A parameter used in the approximation.
      !  B32:    A parameter used in the approximation.
      !  B33:    A parameter used in the approximation.
      !  B34:    A parameter used in the approximation.
      !  B41:    A parameter used in the approximation.
      !  B42:    A parameter used in the approximation.
      !  B43:    A parameter used in the approximation.
      !  B44:    A parameter used in the approximation.
      !  B45:    A parameter used in the approximation.
      !  B51:    A parameter used in the approximation.
      !  B52:    A parameter used in the approximation.
      !  B53:    A parameter used in the approximation.
      !  B54:    A parameter used in the approximation.
      !  B55:    A parameter used in the approximation.
      !  B56:    A parameter used in the approximation.
      !  C:      A value used in the approximation.
      !  CON:    A value used in the approximation.
      !  DF:     The degrees of freedom.
      !  D1:     A value used in the approximation.
      !  D3:     A value used in the approximation.
      !  D5:     A value used in the approximation.
      !  D7:     A value used in the approximation.
      !  D9:     A value used in the approximation.
      !  IDF:    The (positive integer) degrees of freedom.
      !  IPASS:  A value used in the approximation.
      !  MAXIT:  The maximum number of iterations allowed for the approx.
      !  P:      The probability at which the percent point is to be
      !          evaluated.  P must lie between 0.0DO and 1.0E0_wp, exclusive.
      !  PPFN:   The normal percent point value.
      !  S:      A value used in the approximation.
      !  TERM1:  A value used in the approximation.
      !  TERM2:  A value used in the approximation.
      !  TERM3:  A value used in the approximation.
      !  TERM4:  A value used in the approximation.
      !  TERM5:  A value used in the approximation.
      !  Z:      A value used in the approximation.

      df = idf
      maxit = 5

      if (idf <= 0) then
         !Treat the IDF < 1 case
         dpptr = zero

      elseif (idf == 1) then
         !Treat the IDF = 1 (Cauchy) case
         arg = pi*p
         dpptr = -cos(arg)/sin(arg)

      elseif (idf == 2) then
         !  Treat the IDF = 2 case
         term1 = sqrt(two)/two
         term2 = two*p - one
         term3 = sqrt(p*(one - p))
         dpptr = term1*term2/term3

      elseif (idf >= 3) then
         ! Treat the IDF greater than or equal to 3 case
         ppfn = dppnml(p)
         d1 = ppfn
         d3 = ppfn**3
         d5 = ppfn**5
         d7 = ppfn**7
         d9 = ppfn**9
         term1 = d1
         term2 = (one/b21)*(d3 + d1)/df
         term3 = (one/b31)*(b32*d5 + b33*d3 + b34*d1)/(df**2)
         term4 = (one/b41)*(b42*d7 + b43*d5 + b44*d3 + b45*d1)/(df**3)
         term5 = (one/b51)*(b52*d9 + b53*d7 + b54*d5 + b55*d3 + b56*d1)/(df**4)
         dpptr = term1 + term2 + term3 + term4 + term5

         if (idf == 3) then
            ! Augment the results for the IDF = 3 case
            con = pi*(p - half)
            arg = dpptr/sqrt(df)
            z = atan(arg)
            do ipass = 1, maxit
               s = sin(z)
               c = cos(z)
               z = z - (z + s*c - con)/(two*c**2)
            end do
            dpptr = sqrt(df)*s/c

         elseif (idf == 4) then
            ! Augment the results for the IDF = 4 case
            con = two*(p - half)
            arg = dpptr/sqrt(df)
            z = atan(arg)
            do ipass = 1, maxit
               s = sin(z)
               c = cos(z)
               z = z - ((one + half*c**2)*s - con)/((one + half)*c**3)
            end do
            dpptr = sqrt(df)*s/c

         elseif (idf == 5) then
            ! Augment the results for the IDF = 5 case
            con = pi*(p - half)
            arg = dpptr/sqrt(df)
            z = atan(arg)
            do ipass = 1, maxit
               s = sin(z)
               c = cos(z)
               z = z - (z + (c + (two/three)*c**3)*s - con)/ &
                   ((eight/three)*c**4)
            end do
            dpptr = sqrt(df)*s/c

         elseif (idf == 6) then
            !  Augment the results for the IDF = 6 case
            con = two*(p - half)
            arg = dpptr/sqrt(df)
            z = atan(arg)
            do ipass = 1, maxit
               s = sin(z)
               c = cos(z)
               z = z - ((one + half*c**2 + (three/eight)*c**4)*s - con)/ &
                   ((fiftn/eight)*c**5)
            end do
            dpptr = sqrt(df)*s/c
         end if
      end if

   end function dppt