code:
//  (C) Copyright John Maddock 2008.
//  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_SPECIAL_NEXT_HPP
#define BOOST_MATH_SPECIAL_NEXT_HPP

#ifdef _MSC_VER
#pragma once
#endif

#include <boost/math/special_functions/math_fwd.hpp>
#include <boost/math/policies/error_handling.hpp>
#include <boost/math/special_functions/fpclassify.hpp>
#include <boost/math/special_functions/sign.hpp>
#include <boost/math/special_functions/trunc.hpp>

#include <float.h>

#if !defined(_CRAYC) && !defined(__CUDACC__) && (!defined(__GNUC__) || (__GNUC__ > 3) || ((__GNUC__ == 3) && (__GNUC_MINOR__ > 3)))
#if (defined(_M_IX86_FP) && (_M_IX86_FP >= 2)) || defined(__SSE2__)
#include "xmmintrin.h"
#define BOOST_MATH_CHECK_SSE2
#endif
#endif

namespace boost{ namespace math{

namespace detail{

template <class T>
inline T get_smallest_value(mpl::true_ const&)
{
   //
   // numeric_limits lies about denorms being present - particularly
   // when this can be turned on or off at runtime, as is the case
   // when using the SSE2 registers in DAZ or FTZ mode.
   //
   static const T m = std::numeric_limits<T>::denorm_min();
#ifdef BOOST_MATH_CHECK_SSE2
   return (_mm_getcsr() & (_MM_FLUSH_ZERO_ON | 0x40)) ? tools::min_value<T>() : m;;
#else
   return ((tools::min_value<T>() / 2) == 0) ? tools::min_value<T>() : m;
#endif
}

template <class T>
inline T get_smallest_value(mpl::false_ const&)
{
   return tools::min_value<T>();
}

template <class T>
inline T get_smallest_value()
{
#if defined(BOOST_MSVC) && (BOOST_MSVC <= 1310)
   return get_smallest_value<T>(mpl::bool_<std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::has_denorm == 1)>());
#else
   return get_smallest_value<T>(mpl::bool_<std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::has_denorm == std::denorm_present)>());
#endif
}

//
// Returns the smallest value that won't generate denorms when
// we calculate the value of the least-significant-bit:
//
template <class T>
T get_min_shift_value();

template <class T>
struct min_shift_initializer
{
   struct init
   {
      init()
      {
         do_init();
      }
      static void do_init()
      {
         get_min_shift_value<T>();
      }
      void force_instantiate()const{}
   };
   static const init initializer;
   static void force_instantiate()
   {
      initializer.force_instantiate();
   }
};

template <class T>
const typename min_shift_initializer<T>::init min_shift_initializer<T>::initializer;

template <class T>
inline T get_min_shift_value()
{
   BOOST_MATH_STD_USING
   static const T val = ldexp(tools::min_value<T>(), tools::digits<T>() + 1);
   min_shift_initializer<T>::force_instantiate();

   return val;
}

template <class T, class Policy>
T float_next_imp(const T& val, const Policy& pol)
{
   BOOST_MATH_STD_USING
   int expon;
   static const char* function = "float_next<%1%>(%1%)";

   int fpclass = (boost::math::fpclassify)(val);

   if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE))
   {
      if(val < 0)
         return -tools::max_value<T>();
      return policies::raise_domain_error<T>(
         function,
         "Argument must be finite, but got %1%", val, pol);
   }

   if(val >= tools::max_value<T>())
      return policies::raise_overflow_error<T>(function, 0, pol);

   if(val == 0)
      return detail::get_smallest_value<T>();

   if((fpclass != (int)FP_SUBNORMAL) && (fpclass != (int)FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != -tools::min_value<T>()))
   {
      //
      // Special case: if the value of the least significant bit is a denorm, and the result
      // would not be a denorm, then shift the input, increment, and shift back.
      // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
      //
      return ldexp(float_next(T(ldexp(val, 2 * tools::digits<T>())), pol), -2 * tools::digits<T>());
   }

   if(-0.5f == frexp(val, &expon))
      --expon; // reduce exponent when val is a power of two, and negative.
   T diff = ldexp(T(1), expon - tools::digits<T>());
   if(diff == 0)
      diff = detail::get_smallest_value<T>();
   return val + diff;
}

}

template <class T, class Policy>
inline typename tools::promote_args<T>::type float_next(const T& val, const Policy& pol)
{
   typedef typename tools::promote_args<T>::type result_type;
   return detail::float_next_imp(static_cast<result_type>(val), pol);
}

#if 0 //def BOOST_MSVC
//
// We used to use ::_nextafter here, but doing so fails when using
// the SSE2 registers if the FTZ or DAZ flags are set, so use our own
// - albeit slower - code instead as at least that gives the correct answer.
//
template <class Policy>
inline double float_next(const double& val, const Policy& pol)
{
   static const char* function = "float_next<%1%>(%1%)";

   if(!(boost::math::isfinite)(val) && (val > 0))
      return policies::raise_domain_error<double>(
         function,
         "Argument must be finite, but got %1%", val, pol);

   if(val >= tools::max_value<double>())
      return policies::raise_overflow_error<double>(function, 0, pol);

   return ::_nextafter(val, tools::max_value<double>());
}
#endif

template <class T>
inline typename tools::promote_args<T>::type float_next(const T& val)
{
   return float_next(val, policies::policy<>());
}

namespace detail{

template <class T, class Policy>
T float_prior_imp(const T& val, const Policy& pol)
{
   BOOST_MATH_STD_USING
   int expon;
   static const char* function = "float_prior<%1%>(%1%)";

   int fpclass = (boost::math::fpclassify)(val);

   if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE))
   {
      if(val > 0)
         return tools::max_value<T>();
      return policies::raise_domain_error<T>(
         function,
         "Argument must be finite, but got %1%", val, pol);
   }

   if(val <= -tools::max_value<T>())
      return -policies::raise_overflow_error<T>(function, 0, pol);

   if(val == 0)
      return -detail::get_smallest_value<T>();

   if((fpclass != (int)FP_SUBNORMAL) && (fpclass != (int)FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != tools::min_value<T>()))
   {
      //
      // Special case: if the value of the least significant bit is a denorm, and the result
      // would not be a denorm, then shift the input, increment, and shift back.
      // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
      //
      return ldexp(float_prior(T(ldexp(val, 2 * tools::digits<T>())), pol), -2 * tools::digits<T>());
   }

   T remain = frexp(val, &expon);
   if(remain == 0.5)
      --expon; // when val is a power of two we must reduce the exponent
   T diff = ldexp(T(1), expon - tools::digits<T>());
   if(diff == 0)
      diff = detail::get_smallest_value<T>();
   return val - diff;
}

}

template <class T, class Policy>
inline typename tools::promote_args<T>::type float_prior(const T& val, const Policy& pol)
{
   typedef typename tools::promote_args<T>::type result_type;
   return detail::float_prior_imp(static_cast<result_type>(val), pol);
}

#if 0 //def BOOST_MSVC
//
// We used to use ::_nextafter here, but doing so fails when using
// the SSE2 registers if the FTZ or DAZ flags are set, so use our own
// - albeit slower - code instead as at least that gives the correct answer.
//
template <class Policy>
inline double float_prior(const double& val, const Policy& pol)
{
   static const char* function = "float_prior<%1%>(%1%)";

   if(!(boost::math::isfinite)(val) && (val < 0))
      return policies::raise_domain_error<double>(
         function,
         "Argument must be finite, but got %1%", val, pol);

   if(val <= -tools::max_value<double>())
      return -policies::raise_overflow_error<double>(function, 0, pol);

   return ::_nextafter(val, -tools::max_value<double>());
}
#endif

template <class T>
inline typename tools::promote_args<T>::type float_prior(const T& val)
{
   return float_prior(val, policies::policy<>());
}

template <class T, class U, class Policy>
inline typename tools::promote_args<T, U>::type nextafter(const T& val, const U& direction, const Policy& pol)
{
   typedef typename tools::promote_args<T, U>::type result_type;
   return val < direction ? boost::math::float_next<result_type>(val, pol) : val == direction ? val : boost::math::float_prior<result_type>(val, pol);
}

template <class T, class U>
inline typename tools::promote_args<T, U>::type nextafter(const T& val, const U& direction)
{
   return nextafter(val, direction, policies::policy<>());
}

namespace detail{

template <class T, class Policy>
T float_distance_imp(const T& a, const T& b, const Policy& pol)
{
   BOOST_MATH_STD_USING
   //
   // Error handling:
   //
   static const char* function = "float_distance<%1%>(%1%, %1%)";
   if(!(boost::math::isfinite)(a))
      return policies::raise_domain_error<T>(
         function,
         "Argument a must be finite, but got %1%", a, pol);
   if(!(boost::math::isfinite)(b))
      return policies::raise_domain_error<T>(
         function,
         "Argument b must be finite, but got %1%", b, pol);
   //
   // Special cases:
   //
   if(a > b)
      return -float_distance(b, a, pol);
   if(a == b)
      return 0;
   if(a == 0)
      return 1 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol));
   if(b == 0)
      return 1 + fabs(float_distance(static_cast<T>((a < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), a, pol));
   if(boost::math::sign(a) != boost::math::sign(b))
      return 2 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol))
         + fabs(float_distance(static_cast<T>((a < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), a, pol));
   //
   // By the time we get here, both a and b must have the same sign, we want
   // b > a and both postive for the following logic:
   //
   if(a < 0)
      return float_distance(static_cast<T>(-b), static_cast<T>(-a), pol);

   BOOST_ASSERT(a >= 0);
   BOOST_ASSERT(b >= a);

   int expon;
   //
   // Note that if a is a denorm then the usual formula fails
   // because we actually have fewer than tools::digits<T>()
   // significant bits in the representation:
   //
   frexp(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) ? tools::min_value<T>() : a, &expon);
   T upper = ldexp(T(1), expon);
   T result = 0;
   expon = tools::digits<T>() - expon;
   //
   // If b is greater than upper, then we *must* split the calculation
   // as the size of the ULP changes with each order of magnitude change:
   //
   if(b > upper)
   {
      result = float_distance(upper, b);
   }
   //
   // Use compensated double-double addition to avoid rounding
   // errors in the subtraction:
   //
   T mb, x, y, z;
   if(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) || (b - a < tools::min_value<T>()))
   {
      //
      // Special case - either one end of the range is a denormal, or else the difference is.
      // The regular code will fail if we're using the SSE2 registers on Intel and either
      // the FTZ or DAZ flags are set.
      //
      T a2 = ldexp(a, tools::digits<T>());
      T b2 = ldexp(b, tools::digits<T>());
      mb = -(std::min)(T(ldexp(upper, tools::digits<T>())), b2);
      x = a2 + mb;
      z = x - a2;
      y = (a2 - (x - z)) + (mb - z);

      expon -= tools::digits<T>();
   }
   else
   {
      mb = -(std::min)(upper, b);
      x = a + mb;
      z = x - a;
      y = (a - (x - z)) + (mb - z);
   }
   if(x < 0)
   {
      x = -x;
      y = -y;
   }
   result += ldexp(x, expon) + ldexp(y, expon);
   //
   // Result must be an integer:
   //
   BOOST_ASSERT(result == floor(result));
   return result;
}

}

template <class T, class U, class Policy>
inline typename tools::promote_args<T, U>::type float_distance(const T& a, const U& b, const Policy& pol)
{
   typedef typename tools::promote_args<T, U>::type result_type;
   return detail::float_distance_imp(static_cast<result_type>(a), static_cast<result_type>(b), pol);
}

template <class T, class U>
typename tools::promote_args<T, U>::type float_distance(const T& a, const U& b)
{
   return boost::math::float_distance(a, b, policies::policy<>());
}

namespace detail{

template <class T, class Policy>
T float_advance_imp(T val, int distance, const Policy& pol)
{
   BOOST_MATH_STD_USING
   //
   // Error handling:
   //
   static const char* function = "float_advance<%1%>(%1%, int)";

   int fpclass = (boost::math::fpclassify)(val);

   if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE))
      return policies::raise_domain_error<T>(
         function,
         "Argument val must be finite, but got %1%", val, pol);

   if(val < 0)
      return -float_advance(-val, -distance, pol);
   if(distance == 0)
      return val;
   if(distance == 1)
      return float_next(val, pol);
   if(distance == -1)
      return float_prior(val, pol);

   if(fabs(val) < detail::get_min_shift_value<T>())
   {
      //
      // Special case: if the value of the least significant bit is a denorm,
      // implement in terms of float_next/float_prior.
      // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
      //
      if(distance > 0)
      {
         do{ val = float_next(val, pol); } while(--distance);
      }
      else
      {
         do{ val = float_prior(val, pol); } while(++distance);
      }
      return val;
   }

   int expon;
   frexp(val, &expon);
   T limit = ldexp((distance < 0 ? T(0.5f) : T(1)), expon);
   if(val <= tools::min_value<T>())
   {
      limit = sign(T(distance)) * tools::min_value<T>();
   }
   T limit_distance = float_distance(val, limit);
   while(fabs(limit_distance) < abs(distance))
   {
      distance -= itrunc(limit_distance);
      val = limit;
      if(distance < 0)
      {
         limit /= 2;
         expon--;
      }
      else
      {
         limit *= 2;
         expon++;
      }
      limit_distance = float_distance(val, limit);
      if(distance && (limit_distance == 0))
      {
         return policies::raise_evaluation_error<T>(function, "Internal logic failed while trying to increment floating point value %1%: most likely your FPU is in non-IEEE conforming mode.", val, pol);
      }
   }
   if((0.5f == frexp(val, &expon)) && (distance < 0))
      --expon;
   T diff = 0;
   if(val != 0)
      diff = distance * ldexp(T(1), expon - tools::digits<T>());
   if(diff == 0)
      diff = distance * detail::get_smallest_value<T>();
   return val += diff;
}

}

template <class T, class Policy>
inline typename tools::promote_args<T>::type float_advance(T val, int distance, const Policy& pol)
{
   typedef typename tools::promote_args<T>::type result_type;
   return detail::float_advance_imp(static_cast<result_type>(val), distance, pol);
}

template <class T>
inline typename tools::promote_args<T>::type float_advance(const T& val, int distance)
{
   return boost::math::float_advance(val, distance, policies::policy<>());
}

}} // namespaces

#endif // BOOST_MATH_SPECIAL_NEXT_HPP


code:
#ifndef BOOST_MATH_NONFINITE_NUM_FACETS_HPP
#define BOOST_MATH_NONFINITE_NUM_FACETS_HPP

// Copyright 2006 Johan Rade
// Copyright 2012 K R Walker
// Copyright 2011, 2012 Paul A. Bristow

// Distributed under 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)

/*
\file

\brief non_finite_num facets for C99 standard output of infinity and NaN.

\details See fuller documentation at Boost.Math Facets
  for Floating-Point Infinities and NaNs.
*/

#include <cstring>
#include <ios>
#include <limits>
#include <locale>

#include <boost/version.hpp>
#include <boost/throw_exception.hpp>

#include <boost/math/special_functions/fpclassify.hpp>
#include <boost/math/special_functions/sign.hpp>

#ifdef _MSC_VER
#  pragma warning(push)
#  pragma warning(disable : 4127) // conditional expression is constant.
#  pragma warning(disable : 4706) // assignment within conditional expression.
#endif

namespace boost {
  namespace math {

    // flags (enums can be ORed together)       -----------------------------------

    const int legacy = 0x1; //!< get facet will recognize most string representations of infinity and NaN.
    const int signed_zero = 0x2; //!< put facet will distinguish between positive and negative zero.
    const int trap_infinity = 0x4; /*!< put facet will throw an exception of type std::ios_base::failure
       when an attempt is made to format positive or negative infinity.
       get will set the fail bit of the stream when an attempt is made
       to parse a string that represents positive or negative sign infinity.
    */
    const int trap_nan = 0x8; /*!< put facet will throw an exception of type std::ios_base::failure
       when an attempt is made to format positive or negative NaN.
       get will set the fail bit of the stream when an attempt is made
       to parse a string that represents positive or negative sign infinity.
       */

    // class nonfinite_num_put -----------------------------------------------------

    template<
      class CharType,
      class OutputIterator = std::ostreambuf_iterator<CharType>
            >
    class nonfinite_num_put : public std::num_put<CharType, OutputIterator>
    {
    public:
      explicit nonfinite_num_put(int flags = 0) : flags_(flags) {}

    protected:
      virtual OutputIterator do_put(
        OutputIterator it, std::ios_base& iosb, CharType fill, double val) const
      {
        put_and_reset_width(it, iosb, fill, val);
        return it;
      }

      virtual OutputIterator do_put(
        OutputIterator it, std::ios_base& iosb,  CharType fill, long double val) const
      {
        put_and_reset_width(it, iosb, fill, val);
        return it;
      }

    private:
      template<class ValType> void put_and_reset_width(
        OutputIterator& it, std::ios_base& iosb,
        CharType fill, ValType val) const
      {
        put_impl(it, iosb, fill, val);
        iosb.width(0);
      }

      template<class ValType> void put_impl(
        OutputIterator& it, std::ios_base& iosb,
        CharType fill, ValType val) const
      {
        static const CharType prefix_plus[2] = { '+', '\0' };
        static const CharType prefix_minus[2] = { '-', '\0' };
        static const CharType body_inf[4] = { 'i', 'n', 'f', '\0' };
        static const CharType body_nan[4] = { 'n', 'a', 'n', '\0' };
        static const CharType* null_string = 0;

        switch((boost::math::fpclassify)(val))
        {

        case FP_INFINITE:
          if(flags_ & trap_infinity)
          {
            BOOST_THROW_EXCEPTION(std::ios_base::failure("Infinity"));
          }
          else if((boost::math::signbit)(val))
          { // negative infinity.
            put_num_and_fill(it, iosb, prefix_minus, body_inf, fill, val);
          }
          else if(iosb.flags() & std::ios_base::showpos)
          { // Explicit "+inf" wanted.
            put_num_and_fill(it, iosb, prefix_plus, body_inf, fill, val);
          }
          else
          { // just "inf" wanted.
            put_num_and_fill(it, iosb, null_string, body_inf, fill, val);
          }
          break;

        case FP_NAN:
          if(flags_ & trap_nan)
          {
            BOOST_THROW_EXCEPTION(std::ios_base::failure("NaN"));
          }
          else if((boost::math::signbit)(val))
          { // negative so "-nan".
            put_num_and_fill(it, iosb, prefix_minus, body_nan, fill, val);
          }
          else if(iosb.flags() & std::ios_base::showpos)
          { // explicit "+nan" wanted.
            put_num_and_fill(it, iosb, prefix_plus, body_nan, fill, val);
          }
          else
          { // Just "nan".
            put_num_and_fill(it, iosb, null_string, body_nan, fill, val);
          }
          break;

        case FP_ZERO:
          if((flags_ & signed_zero) && ((boost::math::signbit)(val)))
          { // Flag set to distinguish between positive and negative zero.
            // But string "0" should have stuff after decimal point if setprecision and/or exp format.

            std::basic_ostringstream<CharType> zeros; // Needs to be CharType version.

            // Copy flags, fill, width and precision.
            zeros.flags(iosb.flags());
            zeros.unsetf(std::ios::showpos); // Ignore showpos because must be negative.
            zeros.precision(iosb.precision());
            //zeros.width is set by put_num_and_fill
            zeros.fill(static_cast<char>(fill));
            zeros << ValType(0);
            put_num_and_fill(it, iosb, prefix_minus, zeros.str().c_str(), fill, val);
          }
          else
          { // Output the platform default for positive and negative zero.
            put_num_and_fill(it, iosb, null_string, null_string, fill, val);
          }
          break;

        default:  // Normal non-zero finite value.
          it = std::num_put<CharType, OutputIterator>::do_put(it, iosb, fill, val);
          break;
        }
      }

      template<class ValType>
      void put_num_and_fill(
        OutputIterator& it, std::ios_base& iosb, const CharType* prefix,
          const CharType* body, CharType fill, ValType val) const
      {
        int prefix_length = prefix ? (int)std::char_traits<CharType>::length(prefix) : 0;
        int body_length = body ? (int)std::char_traits<CharType>::length(body) : 0;
        int width = prefix_length + body_length;
        std::ios_base::fmtflags adjust = iosb.flags() & std::ios_base::adjustfield;
        const std::ctype<CharType>& ct
          = std::use_facet<std::ctype<CharType> >(iosb.getloc());

        if(body || prefix)
        { // adjust == std::ios_base::right, so leading fill needed.
          if(adjust != std::ios_base::internal && adjust != std::ios_base::left)
            put_fill(it, iosb, fill, width);
        }

        if(prefix)
        { // Adjust width for prefix.
          while(*prefix)
            *it = *(prefix++);
          iosb.width( iosb.width() - prefix_length );
          width -= prefix_length;
        }

        if(body)
        { //
          if(adjust == std::ios_base::internal)
          { // Put fill between sign and digits.
            put_fill(it, iosb, fill, width);
          }
          if(iosb.flags() & std::ios_base::uppercase)
          {
              while(*body)
                *it = ct.toupper(*(body++));
          }
          else
          {
            while(*body)
              *it = *(body++);
          }

          if(adjust == std::ios_base::left)
            put_fill(it, iosb, fill, width);
        }
        else
        {
          it = std::num_put<CharType, OutputIterator>::do_put(it, iosb, fill, val);
        }
      }

      void put_fill(
        OutputIterator& it, std::ios_base& iosb, CharType fill, int width) const
      { // Insert fill chars.
        for(std::streamsize i = iosb.width() - static_cast<std::streamsize>(width); i > 0; --i)
          *it = fill;
      }

    private:
      const int flags_;
    };

    // class nonfinite_num_get ------------------------------------------------------

    template<
      class CharType,
      class InputIterator = std::istreambuf_iterator<CharType>
    >
    class nonfinite_num_get : public std::num_get<CharType, InputIterator>
    {

    public:
      explicit nonfinite_num_get(int flags = 0) : flags_(flags)
      {}

    protected:  // float, double and long double versions of do_get.
      virtual InputIterator do_get(
        InputIterator it, InputIterator end, std::ios_base& iosb,
        std::ios_base::iostate& state, float& val) const
      {
        get_and_check_eof(it, end, iosb, state, val);
        return it;
      }

      virtual InputIterator do_get(
        InputIterator it, InputIterator end, std::ios_base& iosb,
        std::ios_base::iostate& state, double& val) const
      {
        get_and_check_eof(it, end, iosb, state, val);
        return it;
      }

      virtual InputIterator do_get(
        InputIterator it, InputIterator end, std::ios_base& iosb,
        std::ios_base::iostate& state, long double& val) const
      {
        get_and_check_eof(it, end, iosb, state, val);
        return it;
      }

      //..............................................................................

    private:
      template<class ValType> static ValType positive_nan()
      {
        // On some platforms quiet_NaN() may be negative.
        return (boost::math::copysign)(
          std::numeric_limits<ValType>::quiet_NaN(), static_cast<ValType>(1)
          );
        // static_cast<ValType>(1) added Paul A. Bristow 5 Apr 11
      }

      template<class ValType> void get_and_check_eof
      (
        InputIterator& it, InputIterator end, std::ios_base& iosb,
        std::ios_base::iostate& state, ValType& val
      ) const
      {
        get_signed(it, end, iosb, state, val);
        if(it == end)
          state |= std::ios_base::eofbit;
      }

      template<class ValType> void get_signed
      (
        InputIterator& it, InputIterator end, std::ios_base& iosb,
        std::ios_base::iostate& state, ValType& val
      ) const
      {
        const std::ctype<CharType>& ct
          = std::use_facet<std::ctype<CharType> >(iosb.getloc());

        char c = peek_char(it, end, ct);

        bool negative = (c == '-');

        if(negative || c == '+')
        {
          ++it;
          c = peek_char(it, end, ct);
          if(c == '-' || c == '+')
          { // Without this check, "++5" etc would be accepted.
            state |= std::ios_base::failbit;
            return;
          }
        }

        get_unsigned(it, end, iosb, ct, state, val);

        if(negative)
        {
          val = (boost::math::changesign)(val);
        }
      } // void get_signed

      template<class ValType> void get_unsigned
      ( //! Get an unsigned floating-point value into val,
        //! but checking for letters indicating non-finites.
        InputIterator& it, InputIterator end, std::ios_base& iosb,
        const std::ctype<CharType>& ct,
        std::ios_base::iostate& state, ValType& val
      ) const
      {
        switch(peek_char(it, end, ct))
        {
        case 'i':
          get_i(it, end, ct, state, val);
          break;

        case 'n':
          get_n(it, end, ct, state, val);
          break;

        case 'q':
        case 's':
          get_q(it, end, ct, state, val);
          break;

        default: // Got a normal floating-point value into val.
          it = std::num_get<CharType, InputIterator>::do_get(
            it, end, iosb, state, val);
          if((flags_ & legacy) && val == static_cast<ValType>(1)
            && peek_char(it, end, ct) == '#')
            get_one_hash(it, end, ct, state, val);
          break;
        }
      } //  get_unsigned

      //..........................................................................

      template<class ValType> void get_i
      ( // Get the rest of all strings starting with 'i', expect "inf", "infinity".
        InputIterator& it, InputIterator end, const std::ctype<CharType>& ct,
        std::ios_base::iostate& state, ValType& val
      ) const
      {
        if(!std::numeric_limits<ValType>::has_infinity
          || (flags_ & trap_infinity))
        {
            state |= std::ios_base::failbit;
            return;
        }

        ++it;
        if(!match_string(it, end, ct, "nf"))
        {
          state |= std::ios_base::failbit;
          return;
        }

        if(peek_char(it, end, ct) != 'i')
        {
          val = std::numeric_limits<ValType>::infinity();  // "inf"
          return;
        }

        ++it;
        if(!match_string(it, end, ct, "nity"))
        { // Expected "infinity"
          state |= std::ios_base::failbit;
          return;
        }

        val = std::numeric_limits<ValType>::infinity(); // "infinity"
      } // void get_i

      template<class ValType> void get_n
      ( // Get expected strings after 'n', "nan", "nanq", "nans", "nan(...)"
        InputIterator& it, InputIterator end, const std::ctype<CharType>& ct,
        std::ios_base::iostate& state, ValType& val
      ) const
      {
        if(!std::numeric_limits<ValType>::has_quiet_NaN
          || (flags_ & trap_nan)) {
            state |= std::ios_base::failbit;
            return;
        }

        ++it;
        if(!match_string(it, end, ct, "an"))
        {
          state |= std::ios_base::failbit;
          return;
        }

        switch(peek_char(it, end, ct)) {
        case 'q':
        case 's':
          if(flags_ && legacy)
            ++it;
          break;  // "nanq", "nans"

        case '(':   // Optional payload field in (...) follows.
         {
            ++it;
            char c;
            while((c = peek_char(it, end, ct))
              && c != ')' && c != ' ' && c != '\n' && c != '\t')
              ++it;
            if(c != ')')
            { // Optional payload field terminator missing!
              state |= std::ios_base::failbit;
              return;
            }
            ++it;
            break;  // "nan(...)"
          }

        default:
          break;  // "nan"
        }

        val = positive_nan<ValType>();
      } // void get_n

      template<class ValType> void get_q
      ( // Get expected rest of string starting with 'q': "qnan".
        InputIterator& it, InputIterator end, const std::ctype<CharType>& ct,
        std::ios_base::iostate& state, ValType& val
      ) const
      {
        if(!std::numeric_limits<ValType>::has_quiet_NaN
          || (flags_ & trap_nan) || !(flags_ & legacy))
        {
          state |= std::ios_base::failbit;
          return;
        }

        ++it;
        if(!match_string(it, end, ct, "nan"))
        {
          state |= std::ios_base::failbit;
          return;
        }

        val = positive_nan<ValType>(); // "QNAN"
      } //  void get_q

      template<class ValType> void get_one_hash
      ( // Get expected string after having read "1.#": "1.#IND", "1.#QNAN", "1.#SNAN".
        InputIterator& it, InputIterator end, const std::ctype<CharType>& ct,
        std::ios_base::iostate& state, ValType& val
      ) const
      {

        ++it;
        switch(peek_char(it, end, ct))
        {
        case 'i': // from IND (indeterminate), considered same a QNAN.
          get_one_hash_i(it, end, ct, state, val); // "1.#IND"
          return;

        case 'q': // from QNAN
        case 's': // from SNAN - treated the same as QNAN.
          if(std::numeric_limits<ValType>::has_quiet_NaN
            && !(flags_ & trap_nan))
          {
            ++it;
            if(match_string(it, end, ct, "nan"))
            { // "1.#QNAN", "1.#SNAN"
 //             ++it; // removed as caused assert() cannot increment iterator).
// (match_string consumes string, so not needed?).
// https://svn.boost.org/trac/boost/ticket/5467
// Change in nonfinite_num_facet.hpp Paul A. Bristow 11 Apr 11 makes legacy_test.cpp work OK.
              val = positive_nan<ValType>(); // "1.#QNAN"
              return;
            }
          }
          break;

        default:
          break;
        }

        state |= std::ios_base::failbit;
      } //  void get_one_hash

      template<class ValType> void get_one_hash_i
      ( // Get expected strings after 'i', "1.#INF", 1.#IND".
        InputIterator& it, InputIterator end, const std::ctype<CharType>& ct,
        std::ios_base::iostate& state, ValType& val
      ) const
      {
        ++it;

        if(peek_char(it, end, ct) == 'n')
        {
          ++it;
          switch(peek_char(it, end, ct))
          {
          case 'f':  // "1.#INF"
            if(std::numeric_limits<ValType>::has_infinity
              && !(flags_ & trap_infinity))
            {
                ++it;
                val = std::numeric_limits<ValType>::infinity();
                return;
            }
            break;

          case 'd':   // 1.#IND"
            if(std::numeric_limits<ValType>::has_quiet_NaN
              && !(flags_ & trap_nan))
            {
                ++it;
                val = positive_nan<ValType>();
                return;
            }
            break;

          default:
            break;
          }
        }

        state |= std::ios_base::failbit;
      } //  void get_one_hash_i

      //..........................................................................

      char peek_char
      ( //! \return next char in the input buffer, ensuring lowercase (but do not 'consume' char).
        InputIterator& it, InputIterator end,
        const std::ctype<CharType>& ct
      ) const
      {
        if(it == end) return 0;
        return ct.narrow(ct.tolower(*it), 0); // Always tolower to ensure case insensitive.
      }

      bool match_string
      ( //! Match remaining chars to expected string (case insensitive),
        //! consuming chars that match OK.
        //! \return true if matched expected string, else false.
        InputIterator& it, InputIterator end,
        const std::ctype<CharType>& ct,
        const char* s
      ) const
      {
        while(it != end && *s && *s == ct.narrow(ct.tolower(*it), 0))
        {
          ++s;
          ++it; //
        }
        return !*s;
      } // bool match_string

    private:
      const int flags_;
    }; //

    //------------------------------------------------------------------------------

  }   // namespace math
}   // namespace boost

#ifdef _MSC_VER
#   pragma warning(pop)
#endif

#endif // BOOST_MATH_NONFINITE_NUM_FACETS_HPP

license: mit


code:
// Copyright Benjamin Sobotta 2012

//  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_OWENS_T_HPP
#define BOOST_OWENS_T_HPP

// Reference:
// Mike Patefield, David Tandy
// FAST AND ACCURATE CALCULATION OF OWEN'S T-FUNCTION
// Journal of Statistical Software, 5 (5), 1-25

#ifdef _MSC_VER
#  pragma once
#endif

#include <boost/math/special_functions/math_fwd.hpp>
#include <boost/config/no_tr1/cmath.hpp>
#include <boost/math/special_functions/erf.hpp>
#include <boost/math/special_functions/expm1.hpp>
#include <boost/throw_exception.hpp>
#include <boost/assert.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/tools/big_constant.hpp>

#include <stdexcept>

#ifdef BOOST_MSVC
#pragma warning(push)
#pragma warning(disable:4127)
#endif

namespace boost
{
   namespace math
   {
      namespace detail
      {
         // owens_t_znorm1(x) = P(-oo<Z<=x)-0.5 with Z being normally distributed.
         template<typename RealType>
         inline RealType owens_t_znorm1(const RealType x)
         {
            using namespace boost::math::constants;
            return erf(x*one_div_root_two<RealType>())*half<RealType>();
         } // RealType owens_t_znorm1(const RealType x)

         // owens_t_znorm2(x) = P(x<=Z<oo) with Z being normally distributed.
         template<typename RealType>
         inline RealType owens_t_znorm2(const RealType x)
         {
            using namespace boost::math::constants;
            return erfc(x*one_div_root_two<RealType>())*half<RealType>();
         } // RealType owens_t_znorm2(const RealType x)

         // Auxiliary function, it computes an array key that is used to determine
         // the specific computation method for Owen's T and the order thereof
         // used in owens_t_dispatch.
         template<typename RealType>
         inline unsigned short owens_t_compute_code(const RealType h, const RealType a)
         {
            static const RealType hrange[] =
            { 0.02f, 0.06f, 0.09f, 0.125f, 0.26f, 0.4f,  0.6f,  1.6f,  1.7f,  2.33f,  2.4f,  3.36f, 3.4f,  4.8f };

            static const RealType arange[] = { 0.025f, 0.09f, 0.15f, 0.36f, 0.5f, 0.9f, 0.99999f };
            /*
            original select array from paper:
            1, 1, 2,13,13,13,13,13,13,13,13,16,16,16, 9
            1, 2, 2, 3, 3, 5, 5,14,14,15,15,16,16,16, 9
            2, 2, 3, 3, 3, 5, 5,15,15,15,15,16,16,16,10
            2, 2, 3, 5, 5, 5, 5, 7, 7,16,16,16,16,16,10
            2, 3, 3, 5, 5, 6, 6, 8, 8,17,17,17,12,12,11
            2, 3, 5, 5, 5, 6, 6, 8, 8,17,17,17,12,12,12
            2, 3, 4, 4, 6, 6, 8, 8,17,17,17,17,17,12,12
            2, 3, 4, 4, 6, 6,18,18,18,18,17,17,17,12,12
            */
            // subtract one because the array is written in FORTRAN in mind - in C arrays start @ zero
            static const unsigned short select[] =
            {
               0,    0 ,   1  , 12   ,12 ,  12  , 12  , 12 ,  12  , 12  , 12  , 15  , 15 ,  15  ,  8,
               0  ,  1  ,  1   , 2 ,   2   , 4  ,  4  , 13 ,  13  , 14  , 14 ,  15  , 15  , 15  ,  8,
               1  ,  1   , 2 ,   2  ,  2  ,  4   , 4  , 14  , 14 ,  14  , 14 ,  15  , 15 ,  15  ,  9,
               1  ,  1   , 2 ,   4  ,  4  ,  4   , 4  ,  6  ,  6 ,  15  , 15 ,  15 ,  15 ,  15  ,  9,
               1  ,  2   , 2  ,  4  ,  4  ,  5   , 5  ,  7  ,  7  , 16   ,16 ,  16 ,  11 ,  11 ,  10,
               1  ,  2   , 4  ,  4   , 4  ,  5   , 5  ,  7  ,  7  , 16  , 16 ,  16 ,  11  , 11 ,  11,
               1  ,  2   , 3  ,  3  ,  5  ,  5   , 7  ,  7  , 16 ,  16  , 16 ,  16 ,  16  , 11 ,  11,
               1  ,  2   , 3   , 3   , 5  ,  5 ,  17  , 17  , 17 ,  17  , 16 ,  16 ,  16 ,  11 ,  11
            };

            unsigned short ihint = 14, iaint = 7;
            for(unsigned short i = 0; i != 14; i++)
            {
               if( h <= hrange[i] )
               {
                  ihint = i;
                  break;
               }
            } // for(unsigned short i = 0; i != 14; i++)

            for(unsigned short i = 0; i != 7; i++)
            {
               if( a <= arange[i] )
               {
                  iaint = i;
                  break;
               }
            } // for(unsigned short i = 0; i != 7; i++)

            // interprete select array as 8x15 matrix
            return select[iaint*15 + ihint];

         } // unsigned short owens_t_compute_code(const RealType h, const RealType a)

         template<typename RealType>
         inline unsigned short owens_t_get_order_imp(const unsigned short icode, RealType, const mpl::int_<53>&)
         {
            static const unsigned short ord[] = {2, 3, 4, 5, 7, 10, 12, 18, 10, 20, 30, 0, 4, 7, 8, 20, 0, 0}; // 18 entries

            BOOST_ASSERT(icode<18);

            return ord[icode];
         } // unsigned short owens_t_get_order(const unsigned short icode, RealType, mpl::int<53> const&)

         template<typename RealType>
         inline unsigned short owens_t_get_order_imp(const unsigned short icode, RealType, const mpl::int_<64>&)
        {
           // method ================>>>       {1, 1, 1, 1, 1,  1,  1,  1,  2,  2,  2,  3, 4,  4,  4,  4,  5, 6}
           static const unsigned short ord[] = {3, 4, 5, 6, 8, 11, 13, 19, 10, 20, 30,  0, 7, 10, 11, 23,  0, 0}; // 18 entries

          BOOST_ASSERT(icode<18);

          return ord[icode];
        } // unsigned short owens_t_get_order(const unsigned short icode, RealType, mpl::int<64> const&)

         template<typename RealType, typename Policy>
         inline unsigned short owens_t_get_order(const unsigned short icode, RealType r, const Policy&)
         {
            typedef typename policies::precision<RealType, Policy>::type precision_type;
            typedef typename mpl::if_<
               mpl::or_<
                  mpl::less_equal<precision_type, mpl::int_<0> >,
                  mpl::greater<precision_type, mpl::int_<53> >
               >,
               mpl::int_<64>,
               mpl::int_<53>
            >::type tag_type;

            return owens_t_get_order_imp(icode, r, tag_type());
         }

         // compute the value of Owen's T function with method T1 from the reference paper
         template<typename RealType, typename Policy>
         inline RealType owens_t_T1(const RealType h, const RealType a, const unsigned short m, const Policy& pol)
         {
            BOOST_MATH_STD_USING
            using namespace boost::math::constants;

            const RealType hs = -h*h*half<RealType>();
            const RealType dhs = exp( hs );
            const RealType as = a*a;

            unsigned short j=1;
            RealType jj = 1;
            RealType aj = a * one_div_two_pi<RealType>();
            RealType dj = boost::math::expm1( hs, pol);
            RealType gj = hs*dhs;

            RealType val = atan( a ) * one_div_two_pi<RealType>();

            while( true )
            {
               val += dj*aj/jj;

               if( m <= j )
                  break;

               j++;
               jj += static_cast<RealType>(2);
               aj *= as;
               dj = gj - dj;
               gj *= hs / static_cast<RealType>(j);
            } // while( true )

            return val;
         } // RealType owens_t_T1(const RealType h, const RealType a, const unsigned short m)

         // compute the value of Owen's T function with method T2 from the reference paper
         template<typename RealType, class Policy>
         inline RealType owens_t_T2(const RealType h, const RealType a, const unsigned short m, const RealType ah, const Policy&, const mpl::false_&)
         {
            BOOST_MATH_STD_USING
            using namespace boost::math::constants;

            const unsigned short maxii = m+m+1;
            const RealType hs = h*h;
            const RealType as = -a*a;
            const RealType y = static_cast<RealType>(1) / hs;

            unsigned short ii = 1;
            RealType val = 0;
            RealType vi = a * exp( -ah*ah*half<RealType>() ) * one_div_root_two_pi<RealType>();
            RealType z = owens_t_znorm1(ah)/h;

            while( true )
            {
               val += z;
               if( maxii <= ii )
               {
                  val *= exp( -hs*half<RealType>() ) * one_div_root_two_pi<RealType>();
                  break;
               } // if( maxii <= ii )
               z = y * ( vi - static_cast<RealType>(ii) * z );
               vi *= as;
               ii += 2;
            } // while( true )

            return val;
         } // RealType owens_t_T2(const RealType h, const RealType a, const unsigned short m, const RealType ah)

         // compute the value of Owen's T function with method T3 from the reference paper
         template<typename RealType>
         inline RealType owens_t_T3_imp(const RealType h, const RealType a, const RealType ah, const mpl::int_<53>&)
         {
            BOOST_MATH_STD_USING
            using namespace boost::math::constants;

      const unsigned short m = 20;

            static const RealType c2[] =
            {
               static_cast<RealType>(0.99999999999999987510),
               static_cast<RealType>(-0.99999999999988796462),      static_cast<RealType>(0.99999999998290743652),
               static_cast<RealType>(-0.99999999896282500134),      static_cast<RealType>(0.99999996660459362918),
               static_cast<RealType>(-0.99999933986272476760),      static_cast<RealType>(0.99999125611136965852),
               static_cast<RealType>(-0.99991777624463387686),      static_cast<RealType>(0.99942835555870132569),
               static_cast<RealType>(-0.99697311720723000295),      static_cast<RealType>(0.98751448037275303682),
               static_cast<RealType>(-0.95915857980572882813),      static_cast<RealType>(0.89246305511006708555),
               static_cast<RealType>(-0.76893425990463999675),      static_cast<RealType>(0.58893528468484693250),
               static_cast<RealType>(-0.38380345160440256652),      static_cast<RealType>(0.20317601701045299653),
               static_cast<RealType>(-0.82813631607004984866E-01),  static_cast<RealType>(0.24167984735759576523E-01),
               static_cast<RealType>(-0.44676566663971825242E-02),  static_cast<RealType>(0.39141169402373836468E-03)
            };

            const RealType as = a*a;
            const RealType hs = h*h;
            const RealType y = static_cast<RealType>(1)/hs;

            RealType ii = 1;
            unsigned short i = 0;
            RealType vi = a * exp( -ah*ah*half<RealType>() ) * one_div_root_two_pi<RealType>();
            RealType zi = owens_t_znorm1(ah)/h;
            RealType val = 0;

            while( true )
            {
               BOOST_ASSERT(i < 21);
               val += zi*c2[i];
               if( m <= i ) // if( m < i+1 )
               {
                  val *= exp( -hs*half<RealType>() ) * one_div_root_two_pi<RealType>();
                  break;
               } // if( m < i )
               zi = y * (ii*zi - vi);
               vi *= as;
               ii += 2;
               i++;
            } // while( true )

            return val;
         } // RealType owens_t_T3(const RealType h, const RealType a, const RealType ah)

        // compute the value of Owen's T function with method T3 from the reference paper
        template<class RealType>
        inline RealType owens_t_T3_imp(const RealType h, const RealType a, const RealType ah, const mpl::int_<64>&)
        {
          BOOST_MATH_STD_USING
          using namespace boost::math::constants;

          const unsigned short m = 30;

          static const RealType c2[] =
          {
             BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.99999999999999999999999729978162447266851932041876728736094298092917625009873),
             BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.99999999999999999999467056379678391810626533251885323416799874878563998732905968),
             BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.99999999999999999824849349313270659391127814689133077036298754586814091034842536),
             BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.9999999999999997703859616213643405880166422891953033591551179153879839440241685),
             BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.99999999999998394883415238173334565554173013941245103172035286759201504179038147),
             BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.9999999999993063616095509371081203145247992197457263066869044528823599399470977),
             BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.9999999999797336340409464429599229870590160411238245275855903767652432017766116267),
             BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.999999999574958412069046680119051639753412378037565521359444170241346845522403274),
             BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.9999999933226234193375324943920160947158239076786103108097456617750134812033362048),
             BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.9999999188923242461073033481053037468263536806742737922476636768006622772762168467),
             BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.9999992195143483674402853783549420883055129680082932629160081128947764415749728967),
             BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.999993935137206712830997921913316971472227199741857386575097250553105958772041501),
             BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.99996135597690552745362392866517133091672395614263398912807169603795088421057688716),
             BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.99979556366513946026406788969630293820987757758641211293079784585126692672425362469),
             BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.999092789629617100153486251423850590051366661947344315423226082520411961968929483),
             BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.996593837411918202119308620432614600338157335862888580671450938858935084316004769854),
             BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.98910017138386127038463510314625339359073956513420458166238478926511821146316469589567),
             BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.970078558040693314521331982203762771512160168582494513347846407314584943870399016019),
             BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.92911438683263187495758525500033707204091967947532160289872782771388170647150321633673),
             BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.8542058695956156057286980736842905011429254735181323743367879525470479126968822863),
             BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.73796526033030091233118357742803709382964420335559408722681794195743240930748630755),
             BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.58523469882837394570128599003785154144164680587615878645171632791404210655891158),
             BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.415997776145676306165661663581868460503874205343014196580122174949645271353372263),
             BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.2588210875241943574388730510317252236407805082485246378222935376279663808416534365),
             BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.1375535825163892648504646951500265585055789019410617565727090346559210218472356689),
             BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.0607952766325955730493900985022020434830339794955745989150270485056436844239206648),
             BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.0216337683299871528059836483840390514275488679530797294557060229266785853764115),
             BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.00593405693455186729876995814181203900550014220428843483927218267309209471516256),
             BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.0011743414818332946510474576182739210553333860106811865963485870668929503649964142),
             BOOST_MATH_BIG_CONSTANT(RealType, 260, -1.489155613350368934073453260689881330166342484405529981510694514036264969925132e-4),
             BOOST_MATH_BIG_CONSTANT(RealType, 260, 9.072354320794357587710929507988814669454281514268844884841547607134260303118208e-6)
          };

          const RealType as = a*a;
          const RealType hs = h*h;
          const RealType y = 1 / hs;

          RealType ii = 1;
          unsigned short i = 0;
          RealType vi = a * exp( -ah*ah*half<RealType>() ) * one_div_root_two_pi<RealType>();
          RealType zi = owens_t_znorm1(ah)/h;
          RealType val = 0;

          while( true )
          {
              BOOST_ASSERT(i < 31);
              val += zi*c2[i];
              if( m <= i ) // if( m < i+1 )
              {
                val *= exp( -hs*half<RealType>() ) * one_div_root_two_pi<RealType>();
                break;
              } // if( m < i )
              zi = y * (ii*zi - vi);
              vi *= as;
              ii += 2;
              i++;
          } // while( true )

          return val;
        } // RealType owens_t_T3(const RealType h, const RealType a, const RealType ah)

        template<class RealType, class Policy>
        inline RealType owens_t_T3(const RealType h, const RealType a, const RealType ah, const Policy&)
        {
            typedef typename policies::precision<RealType, Policy>::type precision_type;
            typedef typename mpl::if_<
               mpl::or_<
                  mpl::less_equal<precision_type, mpl::int_<0> >,
                  mpl::greater<precision_type, mpl::int_<53> >
               >,
               mpl::int_<64>,
               mpl::int_<53>
            >::type tag_type;

            return owens_t_T3_imp(h, a, ah, tag_type());
        }

         // compute the value of Owen's T function with method T4 from the reference paper
         template<typename RealType>
         inline RealType owens_t_T4(const RealType h, const RealType a, const unsigned short m)
         {
            BOOST_MATH_STD_USING
            using namespace boost::math::constants;

            const unsigned short maxii = m+m+1;
            const RealType hs = h*h;
            const RealType as = -a*a;

            unsigned short ii = 1;
            RealType ai = a * exp( -hs*(static_cast<RealType>(1)-as)*half<RealType>() ) * one_div_two_pi<RealType>();
            RealType yi = 1;
            RealType val = 0;

            while( true )
            {
               val += ai*yi;
               if( maxii <= ii )
                  break;
               ii += 2;
               yi = (static_cast<RealType>(1)-hs*yi) / static_cast<RealType>(ii);
               ai *= as;
            } // while( true )

            return val;
         } // RealType owens_t_T4(const RealType h, const RealType a, const unsigned short m)

         // compute the value of Owen's T function with method T5 from the reference paper
         template<typename RealType>
         inline RealType owens_t_T5_imp(const RealType h, const RealType a, const mpl::int_<53>&)
         {
            BOOST_MATH_STD_USING
            /*
               NOTICE:
               - The pts[] array contains the squares (!) of the abscissas, i.e. the roots of the Legendre
                 polynomial P_n(x), instead of the plain roots as required in Gauss-Legendre
                 quadrature, because T5(h,a,m) contains only x^2 terms.
               - The wts[] array contains the weights for Gauss-Legendre quadrature scaled with a factor
                 of 1/(2*pi) according to T5(h,a,m).
             */

            const unsigned short m = 13;
            static const RealType pts[] = {
               static_cast<RealType>(0.35082039676451715489E-02),
               static_cast<RealType>(0.31279042338030753740E-01),  static_cast<RealType>(0.85266826283219451090E-01),
               static_cast<RealType>(0.16245071730812277011),      static_cast<RealType>(0.25851196049125434828),
               static_cast<RealType>(0.36807553840697533536),      static_cast<RealType>(0.48501092905604697475),
               static_cast<RealType>(0.60277514152618576821),      static_cast<RealType>(0.71477884217753226516),
               static_cast<RealType>(0.81475510988760098605),      static_cast<RealType>(0.89711029755948965867),
               static_cast<RealType>(0.95723808085944261843),      static_cast<RealType>(0.99178832974629703586) };
            static const RealType wts[] = {
               static_cast<RealType>(0.18831438115323502887E-01),
               static_cast<RealType>(0.18567086243977649478E-01),  static_cast<RealType>(0.18042093461223385584E-01),
               static_cast<RealType>(0.17263829606398753364E-01),  static_cast<RealType>(0.16243219975989856730E-01),
               static_cast<RealType>(0.14994592034116704829E-01),  static_cast<RealType>(0.13535474469662088392E-01),
               static_cast<RealType>(0.11886351605820165233E-01),  static_cast<RealType>(0.10070377242777431897E-01),
               static_cast<RealType>(0.81130545742299586629E-02),  static_cast<RealType>(0.60419009528470238773E-02),
               static_cast<RealType>(0.38862217010742057883E-02),  static_cast<RealType>(0.16793031084546090448E-02) };

            const RealType as = a*a;
            const RealType hs = -h*h*boost::math::constants::half<RealType>();

            RealType val = 0;
            for(unsigned short i = 0; i < m; ++i)
            {
               BOOST_ASSERT(i < 13);
               const RealType r = static_cast<RealType>(1) + as*pts[i];
               val += wts[i] * exp( hs*r ) / r;
            } // for(unsigned short i = 0; i < m; ++i)

            return val*a;
         } // RealType owens_t_T5(const RealType h, const RealType a)

        // compute the value of Owen's T function with method T5 from the reference paper
        template<typename RealType>
        inline RealType owens_t_T5_imp(const RealType h, const RealType a, const mpl::int_<64>&)
        {
          BOOST_MATH_STD_USING
            /*
              NOTICE:
              - The pts[] array contains the squares (!) of the abscissas, i.e. the roots of the Legendre
              polynomial P_n(x), instead of the plain roots as required in Gauss-Legendre
              quadrature, because T5(h,a,m) contains only x^2 terms.
              - The wts[] array contains the weights for Gauss-Legendre quadrature scaled with a factor
              of 1/(2*pi) according to T5(h,a,m).
            */

          const unsigned short m = 19;
          static const RealType pts[] = {
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0016634282895983227941),
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.014904509242697054183),
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.04103478879005817919),
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.079359853513391511008),
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.1288612130237615133),
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.18822336642448518856),
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.25586876186122962384),
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.32999972011807857222),
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.40864620815774761438),
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.48971819306044782365),
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.57106118513245543894),
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.6505134942981533829),
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.72596367859928091618),
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.79540665919549865924),
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.85699701386308739244),
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.90909804422384697594),
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.95032536436570154409),
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.97958418733152273717),
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.99610366384229088321)
          };
          static const RealType wts[] = {
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012975111395684900835),
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012888764187499150078),
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012716644398857307844),
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012459897461364705691),
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012120231988292330388),
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.011699908404856841158),
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.011201723906897224448),
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.010628993848522759853),
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0099855296835573320047),
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0092756136096132857933),
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0085039700881139589055),
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0076757344408814561254),
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0067964187616556459109),
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.005871875456524750363),
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0049082589542498110071),
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0039119870792519721409),
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0028897090921170700834),
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0018483371329504443947),
               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.00079623320100438873578)
          };

          const RealType as = a*a;
          const RealType hs = -h*h*boost::math::constants::half<RealType>();

          RealType val = 0;
          for(unsigned short i = 0; i < m; ++i)
            {
              BOOST_ASSERT(i < 19);
              const RealType r = 1 + as*pts[i];
              val += wts[i] * exp( hs*r ) / r;
            } // for(unsigned short i = 0; i < m; ++i)

          return val*a;
        } // RealType owens_t_T5(const RealType h, const RealType a)

        template<class RealType, class Policy>
        inline RealType owens_t_T5(const RealType h, const RealType a, const Policy&)
        {
            typedef typename policies::precision<RealType, Policy>::type precision_type;
            typedef typename mpl::if_<
               mpl::or_<
                  mpl::less_equal<precision_type, mpl::int_<0> >,
                  mpl::greater<precision_type, mpl::int_<53> >
               >,
               mpl::int_<64>,
               mpl::int_<53>
            >::type tag_type;

            return owens_t_T5_imp(h, a, tag_type());
        }

         // compute the value of Owen's T function with method T6 from the reference paper
         template<typename RealType>
         inline RealType owens_t_T6(const RealType h, const RealType a)
         {
            BOOST_MATH_STD_USING
            using namespace boost::math::constants;

            const RealType normh = owens_t_znorm2( h );
            const RealType y = static_cast<RealType>(1) - a;
            const RealType r = atan2(y, static_cast<RealType>(1 + a) );

            RealType val = normh * ( static_cast<RealType>(1) - normh ) * half<RealType>();

            if( r != 0 )
               val -= r * exp( -y*h*h*half<RealType>()/r ) * one_div_two_pi<RealType>();

            return val;
         } // RealType owens_t_T6(const RealType h, const RealType a, const unsigned short m)

         template <class T, class Policy>
         std::pair<T, T> owens_t_T1_accelerated(T h, T a, const Policy& pol)
         {
            //
            // This is the same series as T1, but:
            // * The Taylor series for atan has been combined with that for T1,
            //   reducing but not eliminating cancellation error.
            // * The resulting alternating series is then accelerated using method 1
            //   from H. Cohen, F. Rodriguez Villegas, D. Zagier,
            //   "Convergence acceleration of alternating series", Bonn, (1991).
            //
            BOOST_MATH_STD_USING
            static const char* function = "boost::math::owens_t<%1%>(%1%, %1%)";
            T half_h_h = h * h / 2;
            T a_pow = a;
            T aa = a * a;
            T exp_term = exp(-h * h / 2);
            T one_minus_dj_sum = exp_term;
            T sum = a_pow * exp_term;
            T dj_pow = exp_term;
            T term = sum;
            T abs_err;
            int j = 1;

            //
            // Normally with this form of series acceleration we can calculate
            // up front how many terms will be required - based on the assumption
            // that each term decreases in size by a factor of 3.  However,
            // that assumption does not apply here, as the underlying T1 series can
            // go quite strongly divergent in the early terms, before strongly
            // converging later.  Various "guestimates" have been tried to take account
            // of this, but they don't always work.... so instead set "n" to the
            // largest value that won't cause overflow later, and abort iteration
            // when the last accelerated term was small enough...
            //
            int n;
#ifndef BOOST_NO_EXCEPTIONS
            try
            {
#endif
               n = itrunc(T(tools::log_max_value<T>() / 6));
#ifndef BOOST_NO_EXCEPTIONS
            }
            catch(...)
            {
               n = (std::numeric_limits<int>::max)();
            }
#endif
            n = (std::min)(n, 1500);
            T d = pow(3 + sqrt(T(8)), n);
            d = (d + 1 / d) / 2;
            T b = -1;
            T c = -d;
            c = b - c;
            sum *= c;
            b = -n * n * b * 2;
            abs_err = ldexp(fabs(sum), -tools::digits<T>());

            while(j < n)
            {
               a_pow *= aa;
               dj_pow *= half_h_h / j;
               one_minus_dj_sum += dj_pow;
               term = one_minus_dj_sum * a_pow / (2 * j + 1);
               c = b - c;
               sum += c * term;
               abs_err += ldexp((std::max)(T(fabs(sum)), T(fabs(c*term))), -tools::digits<T>());
               b = (j + n) * (j - n) * b / ((j + T(0.5)) * (j + 1));
               ++j;
               //
               // Include an escape route to prevent calculating too many terms:
               //
               if((j > 10) && (fabs(sum * tools::epsilon<T>()) > fabs(c * term)))
                  break;
            }
            abs_err += fabs(c * term);
            if(sum < 0)  // sum must always be positive, if it's negative something really bad has happend:
               policies::raise_evaluation_error(function, 0, T(0), pol);
            return std::pair<T, T>((sum / d) / boost::math::constants::two_pi<T>(), abs_err / sum);
         }

         template<typename RealType, class Policy>
         inline RealType owens_t_T2(const RealType h, const RealType a, const unsigned short m, const RealType ah, const Policy&, const mpl::true_&)
         {
            BOOST_MATH_STD_USING
            using namespace boost::math::constants;

            const unsigned short maxii = m+m+1;
            const RealType hs = h*h;
            const RealType as = -a*a;
            const RealType y = static_cast<RealType>(1) / hs;

            unsigned short ii = 1;
            RealType val = 0;
            RealType vi = a * exp( -ah*ah*half<RealType>() ) / root_two_pi<RealType>();
            RealType z = owens_t_znorm1(ah)/h;
            RealType last_z = fabs(z);
            RealType lim = policies::get_epsilon<RealType, Policy>();

            while( true )
            {
               val += z;
               //
               // This series stops converging after a while, so put a limit
               // on how far we go before returning our best guess:
               //
               if((fabs(lim * val) > fabs(z)) || ((ii > maxii) && (fabs(z) > last_z)) || (z == 0))
               {
                  val *= exp( -hs*half<RealType>() ) / root_two_pi<RealType>();
                  break;
               } // if( maxii <= ii )
               last_z = fabs(z);
               z = y * ( vi - static_cast<RealType>(ii) * z );
               vi *= as;
               ii += 2;
            } // while( true )

            return val;
         } // RealType owens_t_T2(const RealType h, const RealType a, const unsigned short m, const RealType ah)

         template<typename RealType, class Policy>
         inline std::pair<RealType, RealType> owens_t_T2_accelerated(const RealType h, const RealType a, const RealType ah, const Policy&)
         {
            //
            // This is the same series as T2, but with acceleration applied.
            // Note that we have to be *very* careful to check that nothing bad
            // has happened during evaluation - this series will go divergent
            // and/or fail to alternate at a drop of a hat! :-(
            //
            BOOST_MATH_STD_USING
            using namespace boost::math::constants;

            const RealType hs = h*h;
            const RealType as = -a*a;
            const RealType y = static_cast<RealType>(1) / hs;

            unsigned short ii = 1;
            RealType val = 0;
            RealType vi = a * exp( -ah*ah*half<RealType>() ) / root_two_pi<RealType>();
            RealType z = boost::math::detail::owens_t_znorm1(ah)/h;
            RealType last_z = fabs(z);

            //
            // Normally with this form of series acceleration we can calculate
            // up front how many terms will be required - based on the assumption
            // that each term decreases in size by a factor of 3.  However,
            // that assumption does not apply here, as the underlying T1 series can
            // go quite strongly divergent in the early terms, before strongly
            // converging later.  Various "guestimates" have been tried to take account
            // of this, but they don't always work.... so instead set "n" to the
            // largest value that won't cause overflow later, and abort iteration
            // when the last accelerated term was small enough...
            //
            int n;
#ifndef BOOST_NO_EXCEPTIONS
            try
            {
#endif
               n = itrunc(RealType(tools::log_max_value<RealType>() / 6));
#ifndef BOOST_NO_EXCEPTIONS
            }
            catch(...)
            {
               n = (std::numeric_limits<int>::max)();
            }
#endif
            n = (std::min)(n, 1500);
            RealType d = pow(3 + sqrt(RealType(8)), n);
            d = (d + 1 / d) / 2;
            RealType b = -1;
            RealType c = -d;
            int s = 1;

            for(int k = 0; k < n; ++k)
            {
               //
               // Check for both convergence and whether the series has gone bad:
               //
               if(
                  (fabs(z) > last_z)     // Series has gone divergent, abort
                  || (fabs(val) * tools::epsilon<RealType>() > fabs(c * s * z))  // Convergence!
                  || (z * s < 0)         // Series has stopped alternating - all bets are off - abort.
                  )
               {
                  break;
               }
               c = b - c;
               val += c * s * z;
               b = (k + n) * (k - n) * b / ((k + RealType(0.5)) * (k + 1));
               last_z = fabs(z);
               s = -s;
               z = y * ( vi - static_cast<RealType>(ii) * z );
               vi *= as;
               ii += 2;
            } // while( true )
            RealType err = fabs(c * z) / val;
            return std::pair<RealType, RealType>(val * exp( -hs*half<RealType>() ) / (d * root_two_pi<RealType>()), err);
         } // RealType owens_t_T2_accelerated(const RealType h, const RealType a, const RealType ah, const Policy&)

         template<typename RealType, typename Policy>
         inline RealType T4_mp(const RealType h, const RealType a, const Policy& pol)
         {
            BOOST_MATH_STD_USING

            const RealType hs = h*h;
            const RealType as = -a*a;

            unsigned short ii = 1;
            RealType ai = constants::one_div_two_pi<RealType>() * a * exp( -0.5*hs*(1.0-as) );
            RealType yi = 1.0;
            RealType val = 0.0;

            RealType lim = boost::math::policies::get_epsilon<RealType, Policy>();

            while( true )
            {
               RealType term = ai*yi;
               val += term;
               if((yi != 0) && (fabs(val * lim) > fabs(term)))
                  break;
               ii += 2;
               yi = (1.0-hs*yi) / static_cast<RealType>(ii);
               ai *= as;
               if(ii > (std::min)(1500, (int)policies::get_max_series_iterations<Policy>()))
                  policies::raise_evaluation_error("boost::math::owens_t<%1%>", 0, val, pol);
            } // while( true )

            return val;
         } // arg_type owens_t_T4(const arg_type h, const arg_type a, const unsigned short m)

         // This routine dispatches the call to one of six subroutines, depending on the values
         // of h and a.
         // preconditions: h >= 0, 0<=a<=1, ah=a*h
         //
         // Note there are different versions for different precisions....
         template<typename RealType, typename Policy>
         inline RealType owens_t_dispatch(const RealType h, const RealType a, const RealType ah, const Policy& pol, mpl::int_<64> const&)
         {
            // Simple main case for 64-bit precision or less, this is as per the Patefield-Tandy paper:
            BOOST_MATH_STD_USING
            //
            // Handle some special cases first, these are from
            // page 1077 of Owen's original paper:
            //
            if(h == 0)
            {
               return atan(a) * constants::one_div_two_pi<RealType>();
            }
            if(a == 0)
            {
               return 0;
            }
            if(a == 1)
            {
               return owens_t_znorm2(RealType(-h)) * owens_t_znorm2(h) / 2;
            }
            if(a >= tools::max_value<RealType>())
            {
               return owens_t_znorm2(RealType(fabs(h)));
            }
            RealType val = 0; // avoid compiler warnings, 0 will be overwritten in any case
            const unsigned short icode = owens_t_compute_code(h, a);
            const unsigned short m = owens_t_get_order(icode, val /* just a dummy for the type */, pol);
            static const unsigned short meth[] = {1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 4, 4, 4, 4, 5, 6}; // 18 entries

            // determine the appropriate method, T1 ... T6
            switch( meth[icode] )
            {
            case 1: // T1
               val = owens_t_T1(h,a,m,pol);
               break;
            case 2: // T2
               typedef typename policies::precision<RealType, Policy>::type precision_type;
               typedef mpl::bool_<(precision_type::value == 0) || (precision_type::value > 64)> tag_type;
               val = owens_t_T2(h, a, m, ah, pol, tag_type());
               break;
            case 3: // T3
               val = owens_t_T3(h,a,ah, pol);
               break;
            case 4: // T4
               val = owens_t_T4(h,a,m);
               break;
            case 5: // T5
               val = owens_t_T5(h,a, pol);
               break;
            case 6: // T6
               val = owens_t_T6(h,a);
               break;
            default:
               BOOST_THROW_EXCEPTION(std::logic_error("selection routine in Owen's T function failed"));
            }
            return val;
         }

         template<typename RealType, typename Policy>
         inline RealType owens_t_dispatch(const RealType h, const RealType a, const RealType ah, const Policy& pol, const mpl::int_<65>&)
         {
            // Arbitrary precision version:
            BOOST_MATH_STD_USING
            //
            // Handle some special cases first, these are from
            // page 1077 of Owen's original paper:
            //
            if(h == 0)
            {
               return atan(a) * constants::one_div_two_pi<RealType>();
            }
            if(a == 0)
            {
               return 0;
            }
            if(a == 1)
            {
               return owens_t_znorm2(RealType(-h)) * owens_t_znorm2(h) / 2;
            }
            if(a >= tools::max_value<RealType>())
            {
               return owens_t_znorm2(RealType(fabs(h)));
            }
            // Attempt arbitrary precision code, this will throw if it goes wrong:
            typedef typename boost::math::policies::normalise<Policy, boost::math::policies::evaluation_error<> >::type forwarding_policy;
            std::pair<RealType, RealType> p1(0, tools::max_value<RealType>()), p2(0, tools::max_value<RealType>());
            RealType target_precision = policies::get_epsilon<RealType, Policy>() * 1000;
            bool have_t1(false), have_t2(false);
            if(ah < 3)
            {
#ifndef BOOST_NO_EXCEPTIONS
               try
               {
#endif
                  have_t1 = true;
                  p1 = owens_t_T1_accelerated(h, a, forwarding_policy());
                  if(p1.second < target_precision)
                     return p1.first;
#ifndef BOOST_NO_EXCEPTIONS
               }
               catch(const boost::math::evaluation_error&){}  // T1 may fail and throw, that's OK
#endif
            }
            if(ah > 1)
            {
#ifndef BOOST_NO_EXCEPTIONS
               try
               {
#endif
                  have_t2 = true;
                  p2 = owens_t_T2_accelerated(h, a, ah, forwarding_policy());
                  if(p2.second < target_precision)
                     return p2.first;
#ifndef BOOST_NO_EXCEPTIONS
               }
               catch(const boost::math::evaluation_error&){}  // T2 may fail and throw, that's OK
#endif
            }
            //
            // If we haven't tried T1 yet, do it now - sometimes it succeeds and the number of iterations
            // is fairly low compared to T4.
            //
            if(!have_t1)
            {
#ifndef BOOST_NO_EXCEPTIONS
               try
               {
#endif
                  have_t1 = true;
                  p1 = owens_t_T1_accelerated(h, a, forwarding_policy());
                  if(p1.second < target_precision)
                     return p1.first;
#ifndef BOOST_NO_EXCEPTIONS
               }
               catch(const boost::math::evaluation_error&){}  // T1 may fail and throw, that's OK
#endif
            }
            //
            // If we haven't tried T2 yet, do it now - sometimes it succeeds and the number of iterations
            // is fairly low compared to T4.
            //
            if(!have_t2)
            {
#ifndef BOOST_NO_EXCEPTIONS
               try
               {
#endif
                  have_t2 = true;
                  p2 = owens_t_T2_accelerated(h, a, ah, forwarding_policy());
                  if(p2.second < target_precision)
                     return p2.first;
#ifndef BOOST_NO_EXCEPTIONS
               }
               catch(const boost::math::evaluation_error&){}  // T2 may fail and throw, that's OK
#endif
            }
            //
            // OK, nothing left to do but try the most expensive option which is T4,
            // this is often slow to converge, but when it does converge it tends to
            // be accurate:
#ifndef BOOST_NO_EXCEPTIONS
            try
            {
#endif
               return T4_mp(h, a, pol);
#ifndef BOOST_NO_EXCEPTIONS
            }
            catch(const boost::math::evaluation_error&){}  // T4 may fail and throw, that's OK
#endif
            //
            // Now look back at the results from T1 and T2 and see if either gave better
            // results than we could get from the 64-bit precision versions.
            //
            if((std::min)(p1.second, p2.second) < 1e-20)
            {
               return p1.second < p2.second ? p1.first : p2.first;
            }
            //
            // We give up - no arbitrary precision versions succeeded!
            //
            return owens_t_dispatch(h, a, ah, pol, mpl::int_<64>());
         } // RealType owens_t_dispatch(RealType h, RealType a, RealType ah)
         template<typename RealType, typename Policy>
         inline RealType owens_t_dispatch(const RealType h, const RealType a, const RealType ah, const Policy& pol, const mpl::int_<0>&)
         {
            // We don't know what the precision is until runtime:
            if(tools::digits<RealType>() <= 64)
               return owens_t_dispatch(h, a, ah, pol, mpl::int_<64>());
            return owens_t_dispatch(h, a, ah, pol, mpl::int_<65>());
         }
         template<typename RealType, typename Policy>
         inline RealType owens_t_dispatch(const RealType h, const RealType a, const RealType ah, const Policy& pol)
         {
            // Figure out the precision and forward to the correct version:
            typedef typename policies::precision<RealType, Policy>::type precision_type;
            typedef typename mpl::if_c<
               precision_type::value == 0,
               mpl::int_<0>,
               typename mpl::if_c<
                  precision_type::value <= 64,
                  mpl::int_<64>,
                  mpl::int_<65>
               >::type
            >::type tag_type;
            return owens_t_dispatch(h, a, ah, pol, tag_type());
         }
         // compute Owen's T function, T(h,a), for arbitrary values of h and a
         template<typename RealType, class Policy>
         inline RealType owens_t(RealType h, RealType a, const Policy& pol)
         {
            BOOST_MATH_STD_USING
            // exploit that T(-h,a) == T(h,a)
            h = fabs(h);

            // Use equation (2) in the paper to remap the arguments
            // such that h>=0 and 0<=a<=1 for the call of the actual
            // computation routine.

            const RealType fabs_a = fabs(a);
            const RealType fabs_ah = fabs_a*h;

            RealType val = 0.0; // avoid compiler warnings, 0.0 will be overwritten in any case

            if(fabs_a <= 1)
            {
               val = owens_t_dispatch(h, fabs_a, fabs_ah, pol);
            } // if(fabs_a <= 1.0)
            else
            {
               if( h <= 0.67 )
               {
                  const RealType normh = owens_t_znorm1(h);
                  const RealType normah = owens_t_znorm1(fabs_ah);
                  val = static_cast<RealType>(1)/static_cast<RealType>(4) - normh*normah -
                     owens_t_dispatch(fabs_ah, static_cast<RealType>(1 / fabs_a), h, pol);
               } // if( h <= 0.67 )
               else
               {
                  const RealType normh = detail::owens_t_znorm2(h);
                  const RealType normah = detail::owens_t_znorm2(fabs_ah);
                  val = constants::half<RealType>()*(normh+normah) - normh*normah -
                     owens_t_dispatch(fabs_ah, static_cast<RealType>(1 / fabs_a), h, pol);
               } // else [if( h <= 0.67 )]
            } // else [if(fabs_a <= 1)]

            // exploit that T(h,-a) == -T(h,a)
            if(a < 0)
            {
               return -val;
            } // if(a < 0)

            return val;
         } // RealType owens_t(RealType h, RealType a)

         template <class T, class Policy, class tag>
         struct owens_t_initializer
         {
            struct init
            {
               init()
               {
                  do_init(tag());
               }
               template <int N>
               static void do_init(const mpl::int_<N>&){}
               static void do_init(const mpl::int_<64>&)
               {
                  boost::math::owens_t(static_cast<T>(7), static_cast<T>(0.96875), Policy());
                  boost::math::owens_t(static_cast<T>(2), static_cast<T>(0.5), Policy());
               }
               void force_instantiate()const{}
            };
            static const init initializer;
            static void force_instantiate()
            {
               initializer.force_instantiate();
            }
         };

         template <class T, class Policy, class tag>
         const typename owens_t_initializer<T, Policy, tag>::init owens_t_initializer<T, Policy, tag>::initializer;

      } // namespace detail

      template <class T1, class T2, class Policy>
      inline typename tools::promote_args<T1, T2>::type owens_t(T1 h, T2 a, const Policy& pol)
      {
         typedef typename tools::promote_args<T1, T2>::type result_type;
         typedef typename policies::evaluation<result_type, Policy>::type value_type;
         typedef typename policies::precision<value_type, Policy>::type precision_type;
         typedef typename mpl::if_c<
               precision_type::value == 0,
               mpl::int_<0>,
               typename mpl::if_c<
                  precision_type::value <= 64,
                  mpl::int_<64>,
                  mpl::int_<65>
               >::type
            >::type tag_type;

         detail::owens_t_initializer<result_type, Policy, tag_type>::force_instantiate();

         return policies::checked_narrowing_cast<result_type, Policy>(detail::owens_t(static_cast<value_type>(h), static_cast<value_type>(a), pol), "boost::math::owens_t<%1%>(%1%,%1%)");
      }

      template <class T1, class T2>
      inline typename tools::promote_args<T1, T2>::type owens_t(T1 h, T2 a)
      {
         return owens_t(h, a, policies::policy<>());
      }

   } // namespace math
} // namespace boost

#ifdef BOOST_MSVC
#pragma warning(pop)
#endif

#endif
// EOF


code:
///////////////////////////////////////////////////////////////////////////////
//  Copyright 2013 Nikhar Agrawal
//  Copyright 2013 Christopher Kormanyos
//  Copyright 2014 John Maddock
//  Copyright 2013 Paul Bristow
//  Distributed under 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_POLYGAMMA_2013_07_30_HPP_
  #define _BOOST_POLYGAMMA_2013_07_30_HPP_

#include <boost/math/special_functions/factorials.hpp>
#include <boost/math/special_functions/detail/polygamma.hpp>
#include <boost/math/special_functions/trigamma.hpp>

namespace boost { namespace math {

  template<class T, class Policy>
  inline typename tools::promote_args<T>::type polygamma(const int n, T x, const Policy& pol)
  {
     //
     // Filter off special cases right at the start:
     //
     if(n == 0)
        return boost::math::digamma(x, pol);
     if(n == 1)
        return boost::math::trigamma(x, pol);
     //
     // We've found some standard library functions to misbehave if any FPU exception flags
     // are set prior to their call, this code will clear those flags, then reset them
     // on exit:
     //
     BOOST_FPU_EXCEPTION_GUARD
     //
     // The type of the result - the common type of T and U after
     // any integer types have been promoted to double:
     //
     typedef typename tools::promote_args<T>::type result_type;
     //
     // The type used for the calculation.  This may be a wider type than
     // the result in order to ensure full precision:
     //
     typedef typename policies::evaluation<result_type, Policy>::type value_type;
     //
     // The type of the policy to forward to the actual implementation.
     // We disable promotion of float and double as that's [possibly]
     // happened already in the line above.  Also reset to the default
     // any policies we don't use (reduces code bloat if we're called
     // multiple times with differing policies we don't actually use).
     // Also normalise the type, again to reduce code bloat in case we're
     // called multiple times with functionally identical policies that happen
     // to be different types.
     //
     typedef typename policies::normalise<
        Policy,
        policies::promote_float<false>,
        policies::promote_double<false>,
        policies::discrete_quantile<>,
        policies::assert_undefined<> >::type forwarding_policy;
     //
     // Whew.  Now we can make the actual call to the implementation.
     // Arguments are explicitly cast to the evaluation type, and the result
     // passed through checked_narrowing_cast which handles things like overflow
     // according to the policy passed:
     //
     return policies::checked_narrowing_cast<result_type, forwarding_policy>(
        detail::polygamma_imp(n, static_cast<value_type>(x), forwarding_policy()),
        "boost::math::polygamma<%1%>(int, %1%)");
  }

  template<class T>
  inline typename tools::promote_args<T>::type polygamma(const int n, T x)
  {
      return boost::math::polygamma(n, x, policies::policy<>());
  }

} } // namespace boost::math

#endif // _BOOST_BERNOULLI_2013_05_30_HPP_


code:
//   Boost pow.hpp header file
//   Computes a power with exponent known at compile-time

//  (C) Copyright Bruno Lalande 2008.
//  Distributed under 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)

//  See http://www.boost.org for updates, documentation, and revision history.

#ifndef BOOST_MATH_POW_HPP
#define BOOST_MATH_POW_HPP

#include <boost/math/special_functions/math_fwd.hpp>
#include <boost/math/policies/policy.hpp>
#include <boost/math/policies/error_handling.hpp>
#include <boost/math/tools/promotion.hpp>
#include <boost/mpl/greater_equal.hpp>

namespace boost {
namespace math {

#ifdef BOOST_MSVC
#pragma warning(push)
#pragma warning(disable:4702) // Unreachable code, only triggered in release mode and /W4
#endif

namespace detail {

template <int N, int M = N%2>
struct positive_power
{
    template <typename T>
    static T result(T base)
    {
        T power = positive_power<N/2>::result(base);
        return power * power;
    }
};

template <int N>
struct positive_power<N, 1>
{
    template <typename T>
    static T result(T base)
    {
        T power = positive_power<N/2>::result(base);
        return base * power * power;
    }
};

template <>
struct positive_power<1, 1>
{
    template <typename T>
    static T result(T base){ return base; }
};

template <int N, bool>
struct power_if_positive
{
    template <typename T, class Policy>
    static T result(T base, const Policy&)
    { return positive_power<N>::result(base); }
};

template <int N>
struct power_if_positive<N, false>
{
    template <typename T, class Policy>
    static T result(T base, const Policy& policy)
    {
        if (base == 0)
        {
            return policies::raise_overflow_error<T>(
                       "boost::math::pow(%1%)",
                       "Attempted to compute a negative power of 0",
                       policy
                   );
        }

        return T(1) / positive_power<-N>::result(base);
    }
};

template <>
struct power_if_positive<0, true>
{
    template <typename T, class Policy>
    static T result(T base, const Policy& policy)
    {
        if (base == 0)
        {
            return policies::raise_indeterminate_result_error<T>(
                       "boost::math::pow(%1%)",
                       "The result of pow<0>(%1%) is undetermined",
                       base,
                       T(1),
                       policy
                   );
        }

        return T(1);
    }
};

template <int N>
struct select_power_if_positive
{
    typedef typename mpl::greater_equal<
                         mpl::int_<N>,
                         mpl::int_<0>
                     >::type is_positive;

    typedef power_if_positive<N, is_positive::value> type;
};

}  // namespace detail

template <int N, typename T, class Policy>
inline typename tools::promote_args<T>::type pow(T base, const Policy& policy)
{
   typedef typename tools::promote_args<T>::type result_type;
   return detail::select_power_if_positive<N>::type::result(static_cast<result_type>(base), policy);
}

template <int N, typename T>
inline typename tools::promote_args<T>::type pow(T base)
{ return pow<N>(base, policies::policy<>()); }

#ifdef BOOST_MSVC
#pragma warning(pop)
#endif

}  // namespace math
}  // namespace boost

#endif


code:
//  (C) Copyright John Maddock 2006.
//  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_POWM1
#define BOOST_MATH_POWM1

#ifdef _MSC_VER
#pragma once
#pragma warning(push)
#pragma warning(disable:4702) // Unreachable code (release mode only warning)
#endif

#include <boost/math/special_functions/math_fwd.hpp>
#include <boost/math/special_functions/log1p.hpp>
#include <boost/math/special_functions/expm1.hpp>
#include <boost/math/special_functions/trunc.hpp>
#include <boost/assert.hpp>

namespace boost{ namespace math{ namespace detail{

template <class T, class Policy>
inline T powm1_imp(const T x, const T y, const Policy& pol)
{
   BOOST_MATH_STD_USING
   static const char* function = "boost::math::powm1<%1%>(%1%, %1%)";

   if (x > 0)
   {
      if ((fabs(y * (x - 1)) < 0.5) || (fabs(y) < 0.2))
      {
         // We don't have any good/quick approximation for log(x) * y
         // so just try it and see:
         T l = y * log(x);
         if (l < 0.5)
            return boost::math::expm1(l);
         if (l > boost::math::tools::log_max_value<T>())
            return boost::math::policies::raise_overflow_error<T>(function, 0, pol);
         // fall through....
      }
   }
   else
   {
      // y had better be an integer:
      if (boost::math::trunc(y) != y)
         return boost::math::policies::raise_domain_error<T>(function, "For non-integral exponent, expected base > 0 but got %1%", x, pol);
      if (boost::math::trunc(y / 2) == y / 2)
         return powm1_imp(T(-x), y, pol);
   }
   return pow(x, y) - 1;
}

} // detail

template <class T1, class T2>
inline typename tools::promote_args<T1, T2>::type
   powm1(const T1 a, const T2 z)
{
   typedef typename tools::promote_args<T1, T2>::type result_type;
   return detail::powm1_imp(static_cast<result_type>(a), static_cast<result_type>(z), policies::policy<>());
}

template <class T1, class T2, class Policy>
inline typename tools::promote_args<T1, T2>::type
   powm1(const T1 a, const T2 z, const Policy& pol)
{
   typedef typename tools::promote_args<T1, T2>::type result_type;
   return detail::powm1_imp(static_cast<result_type>(a), static_cast<result_type>(z), pol);
}

} // namespace math
} // namespace boost

#ifdef _MSC_VER
#pragma warning(pop)
#endif

#endif // BOOST_MATH_POWM1