[Boost-commit] svn:boost r4489

From: johnmaddock_at_[hidden]
Date: 2007-06-07 13:42:23

Author: johnmaddock
Date: 2007-06-07 13:42:22 EDT (Thu, 07 Jun 2007)
New Revision: 4489
URL: http://svn.boost.org/trac/boost/changeset/4489

Log:
Performance optimisations: added Cornish-Fisher style expansions to give much better initial approximations to the inverses.  Optimised Student's T quantile so that only one Halley step is taken at double precision.

Added:
   sandbox/math_toolkit/boost/math/special_functions/detail/t_distribution_inv.hpp
Text files modified: 
   sandbox/math_toolkit/boost/math/distributions/binomial.hpp                 |   109 +++++++++++++++++++++++++++++++++++---- 
   sandbox/math_toolkit/boost/math/distributions/students_t.hpp               |    25 ++++++--                                
   sandbox/math_toolkit/boost/math/special_functions/detail/gamma_inva.hpp    |    65 ++++++++++++++++++++++-                 
   sandbox/math_toolkit/boost/math/special_functions/detail/ibeta_inv_ab.hpp  |    80 +++++++++++++++++++++++++++--           
   sandbox/math_toolkit/boost/math/special_functions/detail/ibeta_inverse.hpp |    21 ++++++                                  
   5 files changed, 270 insertions(+), 30 deletions(-)

Modified: sandbox/math_toolkit/boost/math/distributions/binomial.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/distributions/binomial.hpp	(original)
+++ sandbox/math_toolkit/boost/math/distributions/binomial.hpp	2007-06-07 13:42:22 EDT (Thu, 07 Jun 2007)
@@ -160,6 +160,42 @@
               return false;
            return true;
         }
+
+         template <class T>
+         T inverse_binomial_cornish_fisher(T n, T sf, T p, T q)
+         {
+            using namespace std;
+            // mean:
+            T m = n * sf;
+            // standard deviation:
+            T sigma = sqrt(n * sf * (1 - sf));
+            // skewness
+            T sk = (1 - 2 * sf) / sigma;
+            // kurtosis:
+            // T k = (1 - 6 * sf * (1 - sf) ) / (n * sf * (1 - sf));
+            // Get the inverse of a std normal distribution:
+            T x = boost::math::erfc_inv(p > q ? 2 * q : 2 * p) * constants::root_two<T>();
+            // Set the sign:
+            if(p < 0.5)
+               x = -x;
+            T x2 = x * x;
+            // w is correction term due to skewness
+            T w = x + sk * (x2 - 1) / 6;
+            /*
+            // Add on correction due to kurtosis.
+            // Disabled for now, seems to make things worse?
+            //
+            if(n >= 10)
+               w += k * x * (x2 - 3) / 24 + sk * sk * x * (2 * x2 - 5) / -36;
+               */
+            w = m + sigma * w;
+            if(w < tools::min_value<T>())
+               return sqrt(tools::min_value<T>());
+            if(w > n)
+               return n;
+            return w;
+         }
+
      }
 
     // typedef binomial_distribution<double> binomial;
@@ -566,9 +602,10 @@
         // Error checks:
         using namespace std;  // ADL of std names
         RealType result;
+        RealType trials = dist.trials();
         if(false == binomial_detail::check_dist_and_prob(
            BOOST_CURRENT_FUNCTION,
-           dist.trials(),
+           trials,
            dist.success_fraction(),
            p,
            &result))
@@ -589,20 +626,40 @@
            // so n is the most sensible answer here:
            return dist.trials();
         }
-        if (p <= pow(1 - dist.success_fraction(), dist.trials()))
+        if (p <= pow(1 - dist.success_fraction(), trials))
         { // p <= pdf(dist, 0) == cdf(dist, 0)
                                         return 0; // So the only reasonable result is zero.
                                 } // And root finder would fail otherwise.
 
         // Solve for quantile numerically:
         //
+        RealType guess = binomial_detail::inverse_binomial_cornish_fisher(trials, dist.success_fraction(), p, 1-p);
+        RealType factor = 8;
+        if(trials > 100)
+           factor = 1.01f; // guess is pretty accurate
+        else if((trials > 10) && (trials - 1 > guess) && (guess > 3))
+           factor = 1.15f; // less accurate but OK.
+        else if(trials < 10)
+        {
+           // pretty inaccurate guess in this area:
+           if(guess > trials / 64)
+           {
+              guess = trials / 4;
+              factor = 2;
+           }
+           else
+              guess = trials / 1024;
+        }
+        else
+           factor = 2; // trials largish, but in far tails.
+
         detail::binomial_functor<RealType> f(dist, p);
         tools::eps_tolerance<RealType> tol(tools::digits<RealType>());
         boost::uintmax_t max_iter = 1000;
         std::pair<RealType, RealType> r = tools::bracket_and_solve_root(
            f,
-           dist.trials() / 2,
-           static_cast<RealType>(8),
+           guess,
+           factor,
            true,
            tol,
            max_iter);
@@ -621,10 +678,11 @@
         // Error checks:
         RealType q = c.param;
         const binomial_distribution<RealType>& dist = c.dist;
+        RealType trials = dist.trials();
         RealType result;
         if(false == binomial_detail::check_dist_and_prob(
            BOOST_CURRENT_FUNCTION,
-           dist.trials(),
+           trials,
            dist.success_fraction(),
            q,
            &result))
@@ -645,23 +703,50 @@
            // so n is the most sensible answer here:
            return dist.trials();
         }
-
-        if (-q <= powm1(1 - dist.success_fraction(), dist.trials()))
-        { // // q <= cdf(complement(dist, 0)) == pdf(dist, 0)
-          return 0; // So the only reasonable result is zero.
-        } // And root finder would fail otherwise.
+        if(trials == 1)
+        {
+           if(-q <= -dist.success_fraction())
+              return 0;
+        }
+        else
+        {
+           if (-q <= expm1(log1p(-dist.success_fraction()) * trials))
+           { // // q <= cdf(complement(dist, 0)) == pdf(dist, 0)
+             return 0; // So the only reasonable result is zero.
+           } // And root finder would fail otherwise.
+        }
 
         // Need to consider the case where
         //
         // Solve for quantile numerically:
         //
+        RealType guess = binomial_detail::inverse_binomial_cornish_fisher(trials, dist.success_fraction(), 1-q, q);
+        RealType factor = 8;
+        if(trials > 100)
+           factor = 1.01f; // guess is pretty accurate
+        else if((trials > 10) && (trials - 1 > guess) && (guess > 3))
+           factor = 1.15f; // less accurate but OK.
+        else if(trials < 10)
+        {
+           // pretty inaccurate guess in this area:
+           if(guess > trials / 64)
+           {
+              guess = trials / 4;
+              factor = 2;
+           }
+           else
+              guess = trials / 1024;
+        }
+        else
+           factor = 2; // trials largish, but in far tails.
+
         detail::binomial_functor<RealType> f(dist, q, true);
         tools::eps_tolerance<RealType> tol(tools::digits<RealType>());
         boost::uintmax_t max_iter = 1000;
         std::pair<RealType, RealType> r = tools::bracket_and_solve_root(
            f,
-           dist.trials() / 2,
-           static_cast<RealType>(8),
+           guess,
+           factor,
            true,
            tol,
            max_iter);

Modified: sandbox/math_toolkit/boost/math/distributions/students_t.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/distributions/students_t.hpp	(original)
+++ sandbox/math_toolkit/boost/math/distributions/students_t.hpp	2007-06-07 13:42:22 EDT (Thu, 07 Jun 2007)
@@ -85,11 +85,9 @@
    if(false == detail::check_df(
          BOOST_CURRENT_FUNCTION, degrees_of_freedom, &error_result))
       return error_result;
-	 // Might conceivably permit df = +infinity and use normal distribution.
-   // TODO fails for t == 0 and df >=1e16 for ALL fp types.
-   // - probably need to use normal distribution - when available.
-   RealType basem1 = t * t / degrees_of_freedom;
+   // Might conceivably permit df = +infinity and use normal distribution.
    RealType result;
+   RealType basem1 = t * t / degrees_of_freedom;
    if(basem1 < 0.125)
    {
       result = exp(-boost::math::log1p(basem1) * (1+degrees_of_freedom) / 2);
@@ -99,7 +97,6 @@
       result = pow(1 / (1 + basem1), (degrees_of_freedom + 1) / 2);
    }
    result /= sqrt(degrees_of_freedom) * boost::math::beta(degrees_of_freedom / 2, RealType(0.5f));
-
    return result;
 } // pdf
 
@@ -139,12 +136,12 @@
    RealType probability;
    if(degrees_of_freedom > 2 * t2)
    {
-      RealType z = t2 / (degrees_of_freedom + t * t);
+      RealType z = t2 / (degrees_of_freedom + t2);
       probability = ibetac(static_cast<RealType>(0.5), degrees_of_freedom / 2, z) / 2;
    }
    else
    {
-      RealType z = degrees_of_freedom / (degrees_of_freedom + t * t);
+      RealType z = degrees_of_freedom / (degrees_of_freedom + t2);
       probability = ibeta(degrees_of_freedom / 2, static_cast<RealType>(0.5), z) / 2;
    }
    // Check 0 <= probability probability <= 1.
@@ -189,6 +186,11 @@
    if (probability == static_cast<RealType>(0.5))
      return 0;
    //
+   // This next block is disabled in favour of a faster method than
+   // incomplete beta inverse, code retained for future reference:
+   //
+#if 0
+   //
    // Calculate quantile of Student's t using the incomplete beta function inverse:
    //
    probability = (probability > 0.5) ? 1 - probability : probability;
@@ -205,6 +207,15 @@
       t = -t;
 
    return t;
+#endif
+   //
+   // Depending on how many digits RealType has, this may forward
+   // to the incomplete beta inverse as above.  Otherwise uses a
+   // faster method that is accurate to ~15 digits everywhere
+   // and a couple of epsilon at double precision and in the central 
+   // region where most use cases will occur...
+   //
+   return boost::math::detail::fast_students_t_quantile(degrees_of_freedom, probability);
 } // quantile
 
 template <class RealType>

Modified: sandbox/math_toolkit/boost/math/special_functions/detail/gamma_inva.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/detail/gamma_inva.hpp	(original)
+++ sandbox/math_toolkit/boost/math/special_functions/detail/gamma_inva.hpp	2007-06-07 13:42:22 EDT (Thu, 07 Jun 2007)
@@ -32,6 +32,37 @@
 };
 
 template <class T>
+T inverse_poisson_cornish_fisher(T lambda, T p, T q)
+{
+   using namespace std;
+   // mean:
+   T m = lambda;
+   // standard deviation:
+   T sigma = sqrt(lambda);
+   // skewness
+   T sk = 1 / sigma;
+   // kurtosis:
+   // T k = 1/lambda;
+   // Get the inverse of a std normal distribution:
+   T x = boost::math::erfc_inv(p > q ? 2 * q : 2 * p) * constants::root_two<T>();
+   // Set the sign:
+   if(p < 0.5)
+      x = -x;
+   T x2 = x * x;
+   // w is correction term due to skewness
+   T w = x + sk * (x2 - 1) / 6;
+   /*
+   // Add on correction due to kurtosis.
+   // Disabled for now, seems to make things worse?
+   //
+   if(lambda >= 10)
+      w += k * x * (x2 - 3) / 24 + sk * sk * x * (2 * x2 - 5) / -36;
+   */
+   w = m + sigma * w;
+   return w > tools::min_value<T>() ? w : tools::min_value<T>();
+}
+
+template <class T>
 T gamma_inva_imp(const T& z, const T& p, const T& q)
 {
    using namespace std;  // for ADL of std lib math functions
@@ -63,9 +94,31 @@
    // to bracket the root from there:
    //
    T guess;
-   if(z > 1.1)
+   T factor = 8;
+   if(z >= 1)
    {
-      guess = z;
+      //
+      // We can use the relationship between the incomplete 
+      // gamma function and the poisson distribution to
+      // calculate an approximate inverse, for large z
+      // this is actually pretty accurate, but it fails badly
+      // when z is very small.  Also set our step-factor according
+      // to how accurate we think the result is likely to be:
+      //
+      guess = 1 + inverse_poisson_cornish_fisher(z, q, p);
+      if(z > 5)
+      {
+         if(z > 1000)
+            factor = 1.01f;
+         else if(z > 50)
+            factor = 1.1f;
+         else if(guess > 10)
+            factor = 1.25f;
+         else
+            factor = 2;
+         if(guess < 1.1)
+            factor = 8;
+      }
    }
    else if(z > 0.5)
    {
@@ -79,7 +132,13 @@
    // Max iterations permitted:
    //
    boost::uintmax_t max_iter = 200;
-   std::pair<T, T> r = bracket_and_solve_root(f, guess, static_cast<T>(2), false, tol, max_iter);
+   //
+   // Use our generic derivative-free root finding procedure.
+   // We could use Newton steps here, taking the PDF of the
+   // Poisson distribution as our derivative, but that's
+   // even worse performance-wise than the generic method :-(
+   //
+   std::pair<T, T> r = bracket_and_solve_root(f, guess, factor, false, tol, max_iter);
    if(max_iter >= 200)
       tools::logic_error<T>(BOOST_CURRENT_FUNCTION, "Unable to locate the root within a reasonable number of iterations, closest approximation so far was %1%", r.first);
    return (r.first + r.second) / 2;

Modified: sandbox/math_toolkit/boost/math/special_functions/detail/ibeta_inv_ab.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/detail/ibeta_inv_ab.hpp	(original)
+++ sandbox/math_toolkit/boost/math/special_functions/detail/ibeta_inv_ab.hpp	2007-06-07 13:42:22 EDT (Thu, 07 Jun 2007)
@@ -34,6 +34,39 @@
 };
 
 template <class T>
+T inverse_negative_binomial_cornish_fisher(T n, T sf, T sfc, T p, T q)
+{
+   using namespace std;
+   // mean:
+   T m = n * (sfc) / sf;
+   T t = sqrt(n * (sfc));
+   // standard deviation:
+   T sigma = t / sf;
+   // skewness
+   T sk = (1 + sfc) / t;
+   // kurtosis:
+   T k = (6 - sf * (5+sfc)) / (n * (sfc));
+   // Get the inverse of a std normal distribution:
+   T x = boost::math::erfc_inv(p > q ? 2 * q : 2 * p) * constants::root_two<T>();
+   // Set the sign:
+   if(p < 0.5)
+      x = -x;
+   T x2 = x * x;
+   // w is correction term due to skewness
+   T w = x + sk * (x2 - 1) / 6;
+   //
+   // Add on correction due to kurtosis.
+   //
+   if(n >= 10)
+      w += k * x * (x2 - 3) / 24 + sk * sk * x * (2 * x2 - 5) / -36;
+
+   w = m + sigma * w;
+   if(w < tools::min_value<T>())
+      return tools::min_value<T>();
+   return w;
+}
+
+template <class T>
 T ibeta_inv_ab_imp(const T& b, const T& z, const T& p, const T& q, bool swap_ab)
 {
    using namespace std;  // for ADL of std lib math functions
@@ -65,21 +98,56 @@
    // are quite sensitive so we should need too many steps
    // to bracket the root from there:
    //
-   T guess;
-   if((p < q) != swap_ab)
+   T guess = 0;
+   T factor = 5;
+   //
+   // Convert variables to parameters of a negative binomial distribution:
+   //
+   T n = b;
+   T sf = swap_ab ? z : 1-z;
+   T sfc = swap_ab ? 1-z : z;
+   T u = swap_ab ? p : q;
+   T v = swap_ab ? q : p;
+   if(u <= pow(sf, n))
    {
-      guess = b * 2;
+      //
+      // Result is less than 1, negative binomial approximation
+      // is useless....
+      //
+      if((p < q) != swap_ab)
+      {
+         guess = (std::min)(b * 2, T(1));
+      }
+      else
+      {
+         guess = (std::min)(b / 2, T(1));
+      }
    }
-   else
+   if(n * n * n * u * sf > 0.005)
+      guess = 1 + inverse_negative_binomial_cornish_fisher(n, sf, sfc, u, v);
+
+   if(guess < 10)
    {
-      guess = b / 2;
+      //
+      // Negative binomial approximation not accurate in this area:
+      //
+      if((p < q) != swap_ab)
+      {
+         guess = (std::min)(b * 2, T(10));
+      }
+      else
+      {
+         guess = (std::min)(b / 2, T(10));
+      }
    }
+   else
+      factor = (v < sqrt(tools::epsilon<T>())) ? 2 : (guess < 20 ? 1.2f : 1.1f);
    BOOST_MATH_INSTRUMENT_CODE("guess = " << guess);
    //
    // Max iterations permitted:
    //
    boost::uintmax_t max_iter = 200;
-   std::pair<T, T> r = bracket_and_solve_root(f, guess, static_cast<T>(5), swap_ab ? true : false, tol, max_iter);
+   std::pair<T, T> r = bracket_and_solve_root(f, guess, factor, swap_ab ? true : false, tol, max_iter);
    if(max_iter >= 200)
       tools::logic_error<T>(BOOST_CURRENT_FUNCTION, "Unable to locate the root within a reasonable number of iterations, closest approximation so far was %1%", r.first);
    return (r.first + r.second) / 2;

Modified: sandbox/math_toolkit/boost/math/special_functions/detail/ibeta_inverse.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/detail/ibeta_inverse.hpp	(original)
+++ sandbox/math_toolkit/boost/math/special_functions/detail/ibeta_inverse.hpp	2007-06-07 13:42:22 EDT (Thu, 07 Jun 2007)
@@ -9,6 +9,7 @@
 #include <boost/math/special_functions/beta.hpp>
 #include <boost/math/special_functions/erf.hpp>
 #include <boost/math/tools/roots.hpp>
+#include <boost/math/special_functions/detail/t_distribution_inv.hpp>
 
 namespace boost{ namespace math{ namespace detail{
 
@@ -444,7 +445,7 @@
 };
 
 template <class T, class L>
-T ibeta_inv_imp(T a, T b, T p, T q, const L& /* l */, T* py)
+T ibeta_inv_imp(T a, T b, T p, T q, const L& l, T* py)
 {
    using namespace std;  // For ADL of math functions.
 
@@ -465,6 +466,16 @@
    T lower = 0;
    T upper = 1;
    //
+   // Student's T with b = 0.5 gets handled as a special case, swap
+   // around if the arguments are in the "wrong" order:
+   //
+   if(a == 0.5f)
+   {
+      std::swap(a, b);
+      std::swap(p, q);
+      invert = !invert;
+   }
+   //
    // Handle trivial cases first:
    //
    if(q == 0)
@@ -482,6 +493,12 @@
       if(py) *py = 1 - p;
       return p;
    }
+   else if((b == 0.5f) && (a >= 0.5f))
+   {
+      //
+      // We have a Student's T distribution:
+      x = estimate_ibeta_inv_from_t_dist(a, p, q, &y, l);
+   }
    else if(a + b > 5)
    {
       //
@@ -734,7 +751,7 @@
    //
    if(lower == 0)
    {
-      if(invert)
+      if(invert && (py == 0))
       {
          //
          // We're not interested in answers smaller than machine epsilon:

Added: sandbox/math_toolkit/boost/math/special_functions/detail/t_distribution_inv.hpp
==============================================================================
--- (empty file)
+++ sandbox/math_toolkit/boost/math/special_functions/detail/t_distribution_inv.hpp	2007-06-07 13:42:22 EDT (Thu, 07 Jun 2007)
@@ -0,0 +1,502 @@
+//  (C) Copyright John Maddock 2007.
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_MATH_SF_DETAIL_INV_T_HPP
+#define BOOST_MATH_SF_DETAIL_INV_T_HPP
+
+#include <boost/math/special_functions/cbrt.hpp>
+
+namespace boost{ namespace math{ namespace detail{
+
+//
+// The main method used is due to Hill:
+//
+// G. W. Hill, Algorithm 396, Student’s t-Quantiles, 
+// Communications of the ACM, 13(10): 619-620, Oct., 1970.
+//
+template <class T>
+T inverse_students_t_hill(T ndf, T u)
+{
+   using namespace std;
+   BOOST_ASSERT(u <= 0.5);
+   
+   T a, b, c, d, q, x, y;
+
+   if (ndf > 1e20f)
+      return -boost::math::erfc_inv(2 * u) * constants::root_two<T>();
+
+   a = 1 / (ndf - 0.5f);
+   b = 48 / (a * a);
+   c = ((20700 * a / b - 98) * a - 16) * a + 96.36f;
+   d = ((94.5f / (b + c) - 3) / b + 1) * sqrt(a * constants::pi<T>() / 2) * ndf;
+   y = pow(d * 2 * u, 2 / ndf);
+
+   if (y > (0.05f + a)) 
+   {
+      //
+      // Asymptotic inverse expansion about normal:
+      //
+      x = -boost::math::erfc_inv(2 * u) * constants::root_two<T>();
+      y = x * x;
+
+      if (ndf < 5)
+         c += 0.3f * (ndf - 4.5f) * (x + 0.6f);
+      c += (((0.05f * d * x - 5) * x - 7) * x - 2) * x + b;
+      y = (((((0.4f * y + 6.3f) * y + 36) * y + 94.5f) / c - y - 3) / b + 1) * x;
+      y = boost::math::expm1(a * y * y);
+   } 
+   else 
+   {
+      y = ((1 / (((ndf + 6) / (ndf * y) - 0.089f * d - 0.822f) 
+              * (ndf + 2) * 3) + 0.5 / (ndf + 4)) * y - 1) 
+              * (ndf + 1) / (ndf + 2) + 1 / y;
+   }
+   q = sqrt(ndf * y);
+
+   return -q;
+}
+//
+// Tail and body series are due to Shaw:
+//
+// www.mth.kcl.ac.uk/˜shaww/web_page/papers/Tdistribution06.pdf
+//
+// Shaw, W.T., 2006, "Sampling Student’s T distribution – use of 
+// the inverse cumulative distribution function." 
+// Journal of Computational Finance, Vol 9 Issue 4, pp 37-73, Summer 2006
+//
+template <class T, class L>
+T inverse_students_t_tail_series(T df, T v, T u, const L& l)
+{
+   using namespace std;
+   // Tail series expansion, see section 6 of Shaw's paper.
+   // w is calculated using Eq 60:
+   T w = detail::tgamma_delta_ratio_imp(df / 2, constants::half<T>(), l)
+      * sqrt(df * constants::pi<T>()) * v;
+   // define some variables:
+   T np2 = df + 2;
+   T np4 = df + 4;
+   T np6 = df + 6;
+   //
+   // Calculate the coefficients d(k), these depend only on the
+   // number of degrees of freedom df, so at least in theory
+   // we could tabulate these for fixed df, see p15 of Shaw:
+   //
+   T d[7] = { 1, };
+   d[1] = -(df + 1) / (2 * np2);
+   np2 *= (df + 2);
+   d[2] = -df * (df + 1) * (df + 3) / (8 * np2 * np4);
+   np2 *= df + 2;
+   d[3] = -df * (df + 1) * (df + 5) * (((3 * df) + 7) * df -2) / (48 * np2 * np4 * np6);
+   np2 *= (df + 2);
+   np4 *= (df + 4);
+   d[4] = -df * (df + 1) * (df + 7) * 
+      ( (((((15 * df) + 154) * df + 465) * df + 286) * df - 336) * df + 64 ) 
+      / (384 * np2 * np4 * np6 * (df + 8));
+   np2 *= (df + 2);
+   d[5] = -df * (df + 1) * (df + 3) * (df + 9)
+            * (((((((35 * df + 452) * df + 1573) * df + 600) * df - 2020) * df) + 928) * df -128)
+            / (1280 * np2 * np4 * np6 * (df + 8) * (df + 10));
+   np2 *= (df + 2);
+   np4 *= (df + 4);
+   np6 *= (df + 6);
+   d[6] = -df * (df + 1) * (df + 11) 
+            * ((((((((((((945 * df) + 31506) * df + 425858) * df + 2980236) * df + 11266745) * df + 20675018) * df + 7747124) * df - 22574632) * df - 8565600) * df + 18108416) * df - 7099392) * df + 884736)
+            / (46080 * np2 * np4 * np6 * (df + 8) * (df + 10) * (df +12));
+   //
+   // Now bring everthing together to provide the result,
+   // this is Eq 62 of Shaw:
+   //
+   T rn = sqrt(df);
+   T div = pow(rn * w, 1 / df);
+   T power = div * div;
+   T result = tools::evaluate_polynomial(d, power);
+   result *= rn;
+   result /= div;
+   return -result;
+}
+
+template <class T, class L>
+T inverse_students_t_body_series(T df, T u, const L& l)
+{
+   using namespace std;
+   //
+   // Body series for small N:
+   //
+   // Start with Eq 56 of Shaw:
+   //
+   T v = detail::tgamma_delta_ratio_imp(df / 2, constants::half<T>(), l)
+      * sqrt(df * constants::pi<T>()) * (u - constants::half<T>());
+   //
+   // Workspace for the polynomial coefficients:
+   //
+   T c[11] = { 0, 1, };
+   //
+   // Figure out what the coefficients are, note these depend
+   // only on the degrees of freedom (Eq 57 of Shaw):
+   //
+   c[2] = T(1) / 6 + T(1) / (6 * df);
+   T in = 1 / df;
+   c[3] = (((T(1) / 120) * in) + (T(1) / 15)) * in + (T(7) / 120);
+   c[4] = ((((T(1) / 5040) * in + (T(1) / 560)) * in + (T(3) / 112)) * in + T(127) / 5040);
+   c[5] = ((((T(1) / 362880) * in + (T(17) / 45360)) * in + (T(67) / 60480)) * in + (T(479) / 45360)) * in + (T(4369) / 362880);
+   c[6] = ((((((T(1) / 39916800) * in + (T(2503) / 39916800)) * in + (T(11867) / 19958400)) * in + (T(1285) / 798336)) * in + (T(153161) / 39916800)) * in + (T(34807) / 5702400));
+   c[7] = (((((((T(1) / 6227020800LL) * in + (T(37) / 2402400)) * in + 
+      (T(339929) / 2075673600LL)) * in + (T(67217) / 97297200)) * in +
+      (T(870341) / 691891200LL)) * in + (T(70691) / 64864800LL)) * in +
+      (T(20036983LL) / 6227020800LL));
+   c[8] = (((((((T(1) / 1307674368000LL) * in + (T(1042243LL) / 261534873600LL)) * in +
+            (T(21470159) / 435891456000LL)) * in + (T(326228899LL) / 1307674368000LL)) * in +
+            (T(843620579) / 1307674368000LL)) * in + (T(332346031LL) / 435891456000LL)) * in +
+            (T(43847599) / 1307674368000LL)) * in + (T(2280356863LL) / 1307674368000LL);
+   c[9] = (((((((((T(1) / 355687428096000LL)) * in + (T(24262727LL) / 22230464256000LL)) * in +
+            (T(123706507LL) / 8083805184000LL)) * in + (T(404003599LL) / 4446092851200LL)) * in +
+            (T(51811946317LL) / 177843714048000LL)) * in + (T(91423417LL) / 177843714048LL)) * in +
+            (T(32285445833LL) / 88921857024000LL)) * in + (T(531839683LL) / 1710035712000LL)) * in +
+            (T(49020204823LL) / 50812489728000LL);
+   c[10] = (((((((((T(1) / 121645100408832000LL) * in +
+            (T(4222378423LL) / 13516122267648000LL)) * in +
+            (T(49573465457LL) / 10137091700736000LL)) * in +
+            (T(176126809LL) / 5304600576000LL)) * in +
+            (T(44978231873LL) / 355687428096000LL)) * in +
+            (T(5816850595639LL) / 20274183401472000LL)) * in +
+            (T(73989712601LL) / 206879422464000LL)) * in +
+            (T(26591354017LL) / 259925428224000LL)) * in +
+            (T(14979648446341LL) / 40548366802944000LL)) * in +
+            (T(65967241200001LL) / 121645100408832000LL);
+   //
+   // The result is then a polynomial in v (see Eq 56 of Shaw):
+   //
+   return tools::evaluate_odd_polynomial(c, v);
+}
+
+template <class T, class L>
+T inverse_students_t(T df, T u, T v, const L& l, bool* pexact = 0)
+{
+   //
+   // df = number of degrees of freedom.
+   // u = probablity.
+   // v = 1 - u.
+   // l = lanczos type to use.
+   //
+   using namespace std;
+   bool invert = false;
+   T result = 0;
+   if(pexact)
+      *pexact = false;
+   if(u > v)
+   {
+      // function is symmetric, invert it:
+      std::swap(u, v);
+      invert = true;
+   }
+   if((floor(df) == df) && (df < 20))
+   {
+      //
+      // we have integer degrees of freedom, try for the special
+      // cases first:
+      //
+      T tolerance = ldexp(1.0f, (2 * tools::digits<T>()) / 3);
+
+      switch(boost::math::tools::real_cast<int>(df))
+      {
+      case 1:
+         {
+            //
+            // df = 1 is the same as the Cauchy distribution, see
+            // Shaw Eq 35:
+            //
+            // FIXME: fails when u is small!!!
+            result = tan(constants::pi<T>() * (u - constants::half<T>()));
+            if(u == 0.5)
+               result = 0;
+            else
+               result = -cos(constants::pi<T>() * u) / sin(constants::pi<T>() * u);
+            if(pexact)
+               *pexact = true;
+            break;
+         }
+      case 2:
+         {
+            //
+            // df = 2 has an exact result, see Shaw Eq 36:
+            //
+            result =(2 * u - 1) / sqrt(2 * u * v);
+            if(pexact)
+               *pexact = true;
+            break;
+         }
+      case 4:
+         {
+            //
+            // df = 4 has an exact result, see Shaw Eq 38 & 39:
+            //
+            T alpha = 4 * u * v;
+            T root_alpha = sqrt(alpha);
+            T r = 4 * cos(acos(root_alpha) / 3) / root_alpha;
+            T x = sqrt(r - 4);
+            result = u - 0.5f < 0 ? -x : x;
+            if(pexact)
+               *pexact = true;
+            break;
+         }
+      case 6:
+         {
+            //
+            // We get numeric overflow in this area:
+            //
+            if(u < 1e-150)
+               return (invert ? -1 : 1) * inverse_students_t_hill(df, u);
+            //
+            // Newton-Raphson iteration of a polynomial case,
+            // choice of seed value is taken from Shaw's online
+            // supplement: 
+            //
+            T a = 4 * (u - u * u);//1 - 4 * (u - 0.5f) * (u - 0.5f);
+            T b = boost::math::cbrt(a);
+            static const T c = 0.85498797333834849467655443627193L;
+            T p = 6 * (1 + c * (1 / b - 1));
+            T p0;
+            do{
+               T p2 = p * p;
+               T p4 = p2 * p2;
+               T p5 = p * p4;
+               p0 = p;
+               // next term is given by Eq 41:
+               p = 2 * (8 * a * p5 - 270 * p2 + 2187) / (5 * (4 * a * p4 - 216 * p - 243));
+            }while(fabs((p - p0) / p) > tolerance);
+            //
+            // Use Eq 45 to extract the result:
+            //
+            p = sqrt(p - df);
+            result = (u - 0.5f) < 0 ? -p : p;
+            break;
+         }
+#if 0
+         //
+         // These are Shaw's "exact" but iterative solutions
+         // for even df, the numerical accuracy of these is
+         // rather less than Hill's method, so these are disabled
+         // for now, which is a shame because they are reasonably
+         // quick to evaluate...
+         //
+      case 8:
+         {
+            //
+            // Newton-Raphson iteration of a polynomial case,
+            // choice of seed value is taken from Shaw's online
+            // supplement: 
+            //
+            static const T c8 = 0.85994765706259820318168359251872L;
+            T a = 4 * (u - u * u); //1 - 4 * (u - 0.5f) * (u - 0.5f);
+            T b = pow(a, T(1) / 4);
+            T p = 8 * (1 + c8 * (1 / b - 1));
+            T p0 = p;
+            do{
+               T p5 = p * p;
+               p5 *= p5 * p;
+               p0 = p;
+               // Next term is given by Eq 42:
+               p = 2 * (3 * p + (640 * (160 + p * (24 + p * (p + 4)))) / (-5120 + p * (-2048 - 960 * p + a * p5))) / 7;
+            }while(fabs((p - p0) / p) > tolerance); 
+            //
+            // Use Eq 45 to extract the result:
+            //
+            p = sqrt(p - df);
+            result = (u - 0.5f) < 0 ? -p : p;
+            break;
+         }
+      case 10:
+         {
+            //
+            // Newton-Raphson iteration of a polynomial case,
+            // choice of seed value is taken from Shaw's online
+            // supplement: 
+            //
+            static const T c10 = 0.86781292867813396759105692122285L;
+            T a = 4 * (u - u * u); //1 - 4 * (u - 0.5f) * (u - 0.5f);
+            T b = pow(a, T(1) / 5);
+            T p = 10 * (1 + c10 * (1 / b - 1));
+            T p0;
+            do{
+               T p6 = p * p;
+               p6 *= p6 * p6;
+               p0 = p;
+               // Next term given by Eq 43:
+               p = (8 * p) / 9 + (218750 * (21875 + 4 * p * (625 + p * (75 + 2 * p * (5 + p))))) / 
+                  (9 * (-68359375 + 8 * p * (-2343750 + p * (-546875 - 175000 * p + 8 * a * p6))));
+            }while(fabs((p - p0) / p) > tolerance); 
+            //
+            // Use Eq 45 to extract the result:
+            //
+            p = sqrt(p - df);
+            result = (u - 0.5f) < 0 ? -p : p;
+            break;
+         }
+#endif
+      default:
+         goto calculate_real;
+      }
+   }
+   else
+   {
+calculate_real:
+      if(df < 3)
+      {
+         //
+         // Use a roughly linear scheme to choose between Shaw's
+         // tail series and body series:
+         //
+         T crossover = 0.2742f - df * 0.0242143f;
+         if(u > crossover)
+         {
+            result = boost::math::detail::inverse_students_t_body_series(df, u, l);
+         }
+         else
+         {
+            result = boost::math::detail::inverse_students_t_tail_series(df, u, v, l);
+         }
+      }
+      else
+      {
+         //
+         // Use Hill's method except in the exteme tails
+         // where we use Shaw's tail series.
+         // The crossover point is roughly exponential in -df:
+         //
+         T crossover = ldexp(1.0f, tools::real_cast<int>(df / -0.654f));
+         if(u > crossover)
+         {
+            result = boost::math::detail::inverse_students_t_hill(df, u);
+         }
+         else
+         {
+            result = boost::math::detail::inverse_students_t_tail_series(df, u, v, l);
+         }
+      }
+   }
+   return invert ? -result : result;
+}
+
+template <class T, class L>
+inline T estimate_ibeta_inv_from_t_dist(T a, T p, T q, T* py, const L& l)
+{
+   T u = (p > q) ? 0.5f - q / 2 : p / 2;
+   T v = 1 - u; // u < 0.5 so no cancellation error
+   T df = a * 2;
+   T t = boost::math::detail::inverse_students_t(df, u, v, l);
+   T x = df / (df + t * t);
+   *py = t * t / (df + t * t);
+   return x;
+}
+
+template <class T, class L>
+inline T fast_students_t_quantile_imp(T df, T p, const L& l, const mpl::false_*)
+{
+   using namespace std;
+   //
+   // Need to use inverse incomplete beta to get
+   // required precision so not so fast:
+   //
+   T probability = (p > 0.5) ? 1 - p : p;
+   T t, x, y;
+   x = ibeta_inv(df / 2, T(0.5), 2 * probability, &y);
+   if(df * y > tools::max_value<T>() * x)
+      t = tools::overflow_error<T>(BOOST_CURRENT_FUNCTION);
+   else
+      t = sqrt(df * y / x);
+   //
+   // Figure out sign based on the size of p:
+   //
+   if(p < 0.5)
+      t = -t;
+   return t;
+}
+
+template <class T, class L>
+T fast_students_t_quantile_imp(T df, T p, const L& l, const mpl::true_*)
+{
+   using namespace std;
+   bool invert = false;
+   if((df < 2) && (floor(df) != df))
+      return boost::math::detail::fast_students_t_quantile_imp(df, p, l, static_cast<mpl::false_*>(0));
+   if(p > 0.5)
+   {
+      p = 1 - p;
+      invert = true;
+   }
+   //
+   // Get an estimate of the result:
+   //
+   bool exact;
+   T t = inverse_students_t(df, p, 1-p, l, &exact);
+   if((t == 0) || exact)
+      return invert ? -t : t; // can't do better!
+   //
+   // Change variables to inverse incomplete beta:
+   //
+   T t2 = t * t;
+   T xb = df / (df + t2);
+   T y = t2 / (df + t2);
+   T a = df / 2;
+   //
+   // t can be so large that x underflows, 
+   // just return our estimate in that case:
+   //
+   if(xb == 0)
+      return t;
+   //
+   // Get incomplete beta and it's derivative:
+   //
+   T f1;
+   T f0 = xb < y ? ibeta_imp(a, constants::half<T>(), xb, l, false, true, &f1)
+      : ibeta_imp(constants::half<T>(), a, y, l, true, true, &f1);
+
+   // Get cdf from incomplete beta result:
+   T p0 = f0 / 2  - p;
+   // Get pdf from derivative:
+   T p1 = f1 * sqrt(y * xb * xb * xb / df);
+   //
+   // Second derivative divided by p1:
+   //
+   // yacas gives:
+   //
+   // In> PrettyForm(Simplify(D(t) (1 + t^2/v) ^ (-(v+1)/2)))
+   // 
+   //  |                        | v + 1     |     |
+   //  |                       -| ----- + 1 |     |
+   //  |                        |   2       |     |
+   // -|             |  2     |                   |
+   //  |             | t      |                   |
+   //  |             | -- + 1 |                   |
+   //  | ( v + 1 ) * | v      |               * t |
+   // ---------------------------------------------
+   //                       v   
+   //
+   // Which after some manipulation is:
+   //
+   // -p1 * t * (df + 1) / (t^2 + df)
+   //
+   T p2 = t * (df + 1) / (t * t + df);
+   // Halley step:
+   t = fabs(t);
+   t += p0 / (p1 + p0 * p2 / 2);
+   return !invert ? -t : t;
+}
+
+template <class T>
+inline T fast_students_t_quantile(T df, T p)
+{
+   typedef typename lanczos::lanczos_traits<T>::value_type value_type;
+   typedef typename lanczos::lanczos_traits<T>::evaluation_type evaluation_type;
+   typedef mpl::bool_<
+      (std::numeric_limits<T>::digits <= 53)
+       &&
+      (std::numeric_limits<T>::is_specialized)> tag_type;
+   return tools::checked_narrowing_cast<T>(fast_students_t_quantile_imp(static_cast<value_type>(df), static_cast<value_type>(p), evaluation_type(), static_cast<tag_type*>(0)), BOOST_CURRENT_FUNCTION);
+}
+
+}}} // namespaces
+
+#endif // BOOST_MATH_SF_DETAIL_INV_T_HPP