code:
//  Copyright (c) 2013 Christopher Kormanyos
//  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)
//
// This work is based on an earlier work:
// "Algorithm 910: A Portable C++ Multiple-Precision System for Special-Function Calculations",
// in ACM TOMS, {VOL 37, ISSUE 4, (February 2011)} (C) ACM, 2011. http://doi.acm.org/10.1145/1916461.1916469
//
// This header contains implementation details for estimating the zeros
// of cylindrical Bessel and Neumann functions on the positive real axis.
// Support is included for both positive as well as negative order.
// Various methods are used to estimate the roots. These include
// empirical curve fitting and McMahon's asymptotic approximation
// for small order, uniform asymptotic expansion for large order,
// and iteration and root interlacing for negative order.
//
#ifndef _BESSEL_JY_ZERO_2013_01_18_HPP_
  #define _BESSEL_JY_ZERO_2013_01_18_HPP_

  #include <algorithm>
  #include <boost/math/constants/constants.hpp>
  #include <boost/math/special_functions/math_fwd.hpp>
  #include <boost/math/special_functions/cbrt.hpp>
  #include <boost/math/special_functions/detail/airy_ai_bi_zero.hpp>

  namespace boost { namespace math {
  namespace detail
  {
    namespace bessel_zero
    {
      template<class T>
      T equation_nist_10_21_19(const T& v, const T& a)
      {
        // Get the initial estimate of the m'th root of Jv or Yv.
        // This subroutine is used for the order m with m > 1.
        // The order m has been used to create the input parameter a.

        // This is Eq. 10.21.19 in the NIST Handbook.
        const T mu                  = (v * v) * 4U;
        const T mu_minus_one        = mu - T(1);
        const T eight_a_inv         = T(1) / (a * 8U);
        const T eight_a_inv_squared = eight_a_inv * eight_a_inv;

        const T term3 = ((mu_minus_one *  4U) *     ((mu *    7U) -     T(31U) )) / 3U;
        const T term5 = ((mu_minus_one * 32U) *   ((((mu *   83U) -    T(982U) ) * mu) +    T(3779U) )) / 15U;
        const T term7 = ((mu_minus_one * 64U) * ((((((mu * 6949U) - T(153855UL)) * mu) + T(1585743UL)) * mu) - T(6277237UL))) / 105U;

        return a + ((((                      - term7
                       * eight_a_inv_squared - term5)
                       * eight_a_inv_squared - term3)
                       * eight_a_inv_squared - mu_minus_one)
                       * eight_a_inv);
      }

      template<typename T>
      class equation_as_9_3_39_and_its_derivative
      {
      public:
        equation_as_9_3_39_and_its_derivative(const T& zt) : zeta(zt) { }

        boost::math::tuple<T, T> operator()(const T& z) const
        {
          BOOST_MATH_STD_USING // ADL of std names, needed for acos, sqrt.

          // Return the function of zeta that is implicitly defined
          // in A&S Eq. 9.3.39 as a function of z. The function is
          // returned along with its derivative with respect to z.

          const T zsq_minus_one_sqrt = sqrt((z * z) - T(1));

          const T the_function(
              zsq_minus_one_sqrt
            - (  acos(T(1) / z) + ((T(2) / 3U) * (zeta * sqrt(zeta)))));

          const T its_derivative(zsq_minus_one_sqrt / z);

          return boost::math::tuple<T, T>(the_function, its_derivative);
        }

      private:
        const equation_as_9_3_39_and_its_derivative& operator=(const equation_as_9_3_39_and_its_derivative&);
        const T zeta;
      };

      template<class T>
      static T equation_as_9_5_26(const T& v, const T& ai_bi_root)
      {
        BOOST_MATH_STD_USING // ADL of std names, needed for log, sqrt.

        // Obtain the estimate of the m'th zero of Jv or Yv.
        // The order m has been used to create the input parameter ai_bi_root.
        // Here, v is larger than about 2.2. The estimate is computed
        // from Abramowitz and Stegun Eqs. 9.5.22 and 9.5.26, page 371.
        //
        // The inversion of z as a function of zeta is mentioned in the text
        // following A&S Eq. 9.5.26. Here, we accomplish the inversion by
        // performing a Taylor expansion of Eq. 9.3.39 for large z to order 2
        // and solving the resulting quadratic equation, thereby taking
        // the positive root of the quadratic.
        // In other words: (2/3)(-zeta)^(3/2) approx = z + 1/(2z) - pi/2.
        // This leads to: z^2 - [(2/3)(-zeta)^(3/2) + pi/2]z + 1/2 = 0.
        //
        // With this initial estimate, Newton-Raphson iteration is used
        // to refine the value of the estimate of the root of z
        // as a function of zeta.

        const T v_pow_third(boost::math::cbrt(v));
        const T v_pow_minus_two_thirds(T(1) / (v_pow_third * v_pow_third));

        // Obtain zeta using the order v combined with the m'th root of
        // an airy function, as shown in  A&S Eq. 9.5.22.
        const T zeta = v_pow_minus_two_thirds * (-ai_bi_root);

        const T zeta_sqrt = sqrt(zeta);

        // Set up a quadratic equation based on the Taylor series
        // expansion mentioned above.
        const T b = -((((zeta * zeta_sqrt) * 2U) / 3U) + boost::math::constants::half_pi<T>());

        // Solve the quadratic equation, taking the positive root.
        const T z_estimate = (-b + sqrt((b * b) - T(2))) / 2U;

        // Establish the range, the digits, and the iteration limit
        // for the upcoming root-finding.
        const T range_zmin = (std::max<T>)(z_estimate - T(1), T(1));
        const T range_zmax = z_estimate + T(1);

        const int my_digits10 = static_cast<int>(static_cast<float>(boost::math::tools::digits<T>() * 0.301F));

        // Select the maximum allowed iterations based on the number
        // of decimal digits in the numeric type T, being at least 12.
        const boost::uintmax_t iterations_allowed = static_cast<boost::uintmax_t>((std::max)(12, my_digits10 * 2));

        boost::uintmax_t iterations_used = iterations_allowed;

        // Calculate the root of z as a function of zeta.
        const T z = boost::math::tools::newton_raphson_iterate(
          boost::math::detail::bessel_zero::equation_as_9_3_39_and_its_derivative<T>(zeta),
          z_estimate,
          range_zmin,
          range_zmax,
          (std::min)(boost::math::tools::digits<T>(), boost::math::tools::digits<float>()),
          iterations_used);

        static_cast<void>(iterations_used);

        // Continue with the implementation of A&S Eq. 9.3.39.
        const T zsq_minus_one      = (z * z) - T(1);
        const T zsq_minus_one_sqrt = sqrt(zsq_minus_one);

        // This is A&S Eq. 9.3.42.
        const T b0_term_5_24 = T(5) / ((zsq_minus_one * zsq_minus_one_sqrt) * 24U);
        const T b0_term_1_8  = T(1) / ( zsq_minus_one_sqrt * 8U);
        const T b0_term_5_48 = T(5) / ((zeta * zeta) * 48U);

        const T b0 = -b0_term_5_48 + ((b0_term_5_24 + b0_term_1_8) / zeta_sqrt);

        // This is the second line of A&S Eq. 9.5.26 for f_k with k = 1.
        const T f1 = ((z * zeta_sqrt) * b0) / zsq_minus_one_sqrt;

        // This is A&S Eq. 9.5.22 expanded to k = 1 (i.e., one term in the series).
        return (v * z) + (f1 / v);
      }

      namespace cyl_bessel_j_zero_detail
      {
        template<class T>
        T equation_nist_10_21_40_a(const T& v)
        {
          const T v_pow_third(boost::math::cbrt(v));
          const T v_pow_minus_two_thirds(T(1) / (v_pow_third * v_pow_third));

          return v * (((((                         + T(0.043)
                          * v_pow_minus_two_thirds - T(0.0908))
                          * v_pow_minus_two_thirds - T(0.00397))
                          * v_pow_minus_two_thirds + T(1.033150))
                          * v_pow_minus_two_thirds + T(1.8557571))
                          * v_pow_minus_two_thirds + T(1));
        }

        template<class T, class Policy>
        class function_object_jv
        {
        public:
          function_object_jv(const T& v,
                             const Policy& pol) : my_v(v),
                                                  my_pol(pol) { }

          T operator()(const T& x) const
          {
            return boost::math::cyl_bessel_j(my_v, x, my_pol);
          }

        private:
          const T my_v;
          const Policy& my_pol;
          const function_object_jv& operator=(const function_object_jv&);
        };

        template<class T, class Policy>
        class function_object_jv_and_jv_prime
        {
        public:
          function_object_jv_and_jv_prime(const T& v,
                                          const bool order_is_zero,
                                          const Policy& pol) : my_v(v),
                                                               my_order_is_zero(order_is_zero),
                                                               my_pol(pol) { }

          boost::math::tuple<T, T> operator()(const T& x) const
          {
            // Obtain Jv(x) and Jv'(x).
            // Chris's original code called the Bessel function implementation layer direct,
            // but that circumvented optimizations for integer-orders.  Call the documented
            // top level functions instead, and let them sort out which implementation to use.
            T j_v;
            T j_v_prime;

            if(my_order_is_zero)
            {
              j_v       =  boost::math::cyl_bessel_j(0, x, my_pol);
              j_v_prime = -boost::math::cyl_bessel_j(1, x, my_pol);
            }
            else
            {
                      j_v       = boost::math::cyl_bessel_j(  my_v,      x, my_pol);
              const T j_v_m1     (boost::math::cyl_bessel_j(T(my_v - 1), x, my_pol));
                      j_v_prime = j_v_m1 - ((my_v * j_v) / x);
            }

            // Return a tuple containing both Jv(x) and Jv'(x).
            return boost::math::make_tuple(j_v, j_v_prime);
          }

        private:
          const T my_v;
          const bool my_order_is_zero;
          const Policy& my_pol;
          const function_object_jv_and_jv_prime& operator=(const function_object_jv_and_jv_prime&);
        };

        template<class T> bool my_bisection_unreachable_tolerance(const T&, const T&) { return false; }

        template<class T, class Policy>
        T initial_guess(const T& v, const int m, const Policy& pol)
        {
          BOOST_MATH_STD_USING // ADL of std names, needed for floor.

          // Compute an estimate of the m'th root of cyl_bessel_j.

          T guess;

          // There is special handling for negative order.
          if(v < 0)
          {
            if((m == 1) && (v > -0.5F))
            {
              // For small, negative v, use the results of empirical curve fitting.
              // Mathematica(R) session for the coefficients:
              //  Table[{n, BesselJZero[n, 1]}, {n, -(1/2), 0, 1/10}]
              //  N[%, 20]
              //  Fit[%, {n^0, n^1, n^2, n^3, n^4, n^5, n^6}, n]
              guess = (((((    - T(0.2321156900729)
                           * v - T(0.1493247777488))
                           * v - T(0.15205419167239))
                           * v + T(0.07814930561249))
                           * v - T(0.17757573537688))
                           * v + T(1.542805677045663))
                           * v + T(2.40482555769577277);

              return guess;
            }

            // Create the positive order and extract its positive floor integer part.
            const T vv(-v);
            const T vv_floor(floor(vv));

            // The to-be-found root is bracketed by the roots of the
            // Bessel function whose reflected, positive integer order
            // is less than, but nearest to vv.

            T root_hi = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess(vv_floor, m, pol);
            T root_lo;

            if(m == 1)
            {
              // The estimate of the first root for negative order is found using
              // an adaptive range-searching algorithm.
              root_lo = T(root_hi - 0.1F);

              const bool hi_end_of_bracket_is_negative = (boost::math::cyl_bessel_j(v, root_hi, pol) < 0);

              while((root_lo > boost::math::tools::epsilon<T>()))
              {
                const bool lo_end_of_bracket_is_negative = (boost::math::cyl_bessel_j(v, root_lo, pol) < 0);

                if(hi_end_of_bracket_is_negative != lo_end_of_bracket_is_negative)
                {
                  break;
                }

                root_hi = root_lo;

                // Decrease the lower end of the bracket using an adaptive algorithm.
                if(root_lo > 0.5F)
                {
                  root_lo -= 0.5F;
                }
                else
                {
                  root_lo *= 0.75F;
                }
              }
            }
            else
            {
              root_lo = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess(vv_floor, m - 1, pol);
            }

            // Perform several steps of bisection iteration to refine the guess.
            boost::uintmax_t number_of_iterations(12U);

            // Do the bisection iteration.
            const boost::math::tuple<T, T> guess_pair =
               boost::math::tools::bisect(
                  boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::function_object_jv<T, Policy>(v, pol),
                  root_lo,
                  root_hi,
                  boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::my_bisection_unreachable_tolerance<T>,
                  number_of_iterations);

            return (boost::math::get<0>(guess_pair) + boost::math::get<1>(guess_pair)) / 2U;
          }

          if(m == 1U)
          {
            // Get the initial estimate of the first root.

            if(v < 2.2F)
            {
              // For small v, use the results of empirical curve fitting.
              // Mathematica(R) session for the coefficients:
              //  Table[{n, BesselJZero[n, 1]}, {n, 0, 22/10, 1/10}]
              //  N[%, 20]
              //  Fit[%, {n^0, n^1, n^2, n^3, n^4, n^5, n^6}, n]
              guess = (((((    - T(0.0008342379046010)
                           * v + T(0.007590035637410))
                           * v - T(0.030640914772013))
                           * v + T(0.078232088020106))
                           * v - T(0.169668712590620))
                           * v + T(1.542187960073750))
                           * v + T(2.4048359915254634);
            }
            else
            {
              // For larger v, use the first line of Eqs. 10.21.40 in the NIST Handbook.
              guess = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::equation_nist_10_21_40_a(v);
            }
          }
          else
          {
            if(v < 2.2F)
            {
              // Use Eq. 10.21.19 in the NIST Handbook.
              const T a(((v + T(m * 2U)) - T(0.5)) * boost::math::constants::half_pi<T>());

              guess = boost::math::detail::bessel_zero::equation_nist_10_21_19(v, a);
            }
            else
            {
              // Get an estimate of the m'th root of airy_ai.
              const T airy_ai_root(boost::math::detail::airy_zero::airy_ai_zero_detail::initial_guess<T>(m));

              // Use Eq. 9.5.26 in the A&S Handbook.
              guess = boost::math::detail::bessel_zero::equation_as_9_5_26(v, airy_ai_root);
            }
          }

          return guess;
        }
      } // namespace cyl_bessel_j_zero_detail

      namespace cyl_neumann_zero_detail
      {
        template<class T>
        T equation_nist_10_21_40_b(const T& v)
        {
          const T v_pow_third(boost::math::cbrt(v));
          const T v_pow_minus_two_thirds(T(1) / (v_pow_third * v_pow_third));

          return v * (((((                         - T(0.001)
                          * v_pow_minus_two_thirds - T(0.0060))
                          * v_pow_minus_two_thirds + T(0.01198))
                          * v_pow_minus_two_thirds + T(0.260351))
                          * v_pow_minus_two_thirds + T(0.9315768))
                          * v_pow_minus_two_thirds + T(1));
        }

        template<class T, class Policy>
        class function_object_yv
        {
        public:
          function_object_yv(const T& v,
                             const Policy& pol) : my_v(v),
                                                  my_pol(pol) { }

          T operator()(const T& x) const
          {
            return boost::math::cyl_neumann(my_v, x, my_pol);
          }

        private:
          const T my_v;
          const Policy& my_pol;
          const function_object_yv& operator=(const function_object_yv&);
        };

        template<class T, class Policy>
        class function_object_yv_and_yv_prime
        {
        public:
          function_object_yv_and_yv_prime(const T& v,
                                          const Policy& pol) : my_v(v),
                                                               my_pol(pol) { }

          boost::math::tuple<T, T> operator()(const T& x) const
          {
            const T half_epsilon(boost::math::tools::epsilon<T>() / 2U);

            const bool order_is_zero = ((my_v > -half_epsilon) && (my_v < +half_epsilon));

            // Obtain Yv(x) and Yv'(x).
            // Chris's original code called the Bessel function implementation layer direct,
            // but that circumvented optimizations for integer-orders.  Call the documented
            // top level functions instead, and let them sort out which implementation to use.
            T y_v;
            T y_v_prime;

            if(order_is_zero)
            {
              y_v       =  boost::math::cyl_neumann(0, x, my_pol);
              y_v_prime = -boost::math::cyl_neumann(1, x, my_pol);
            }
            else
            {
                      y_v       = boost::math::cyl_neumann(  my_v,      x, my_pol);
              const T y_v_m1     (boost::math::cyl_neumann(T(my_v - 1), x, my_pol));
                      y_v_prime = y_v_m1 - ((my_v * y_v) / x);
            }

            // Return a tuple containing both Yv(x) and Yv'(x).
            return boost::math::make_tuple(y_v, y_v_prime);
          }

        private:
          const T my_v;
          const Policy& my_pol;
          const function_object_yv_and_yv_prime& operator=(const function_object_yv_and_yv_prime&);
        };

        template<class T> bool my_bisection_unreachable_tolerance(const T&, const T&) { return false; }

        template<class T, class Policy>
        T initial_guess(const T& v, const int m, const Policy& pol)
        {
          BOOST_MATH_STD_USING // ADL of std names, needed for floor.

          // Compute an estimate of the m'th root of cyl_neumann.

          T guess;

          // There is special handling for negative order.
          if(v < 0)
          {
            // Create the positive order and extract its positive floor and ceiling integer parts.
            const T vv(-v);
            const T vv_floor(floor(vv));

            // The to-be-found root is bracketed by the roots of the
            // Bessel function whose reflected, positive integer order
            // is less than, but nearest to vv.

            // The special case of negative, half-integer order uses
            // the relation between Yv and spherical Bessel functions
            // in order to obtain the bracket for the root.
            // In these special cases, cyl_neumann(-n/2, x) = sph_bessel_j(+n/2, x)
            // for v = -n/2.

            T root_hi;
            T root_lo;

            if(m == 1)
            {
              // The estimate of the first root for negative order is found using
              // an adaptive range-searching algorithm.
              // Take special precautions for the discontinuity at negative,
              // half-integer orders and use different brackets above and below these.
              if(T(vv - vv_floor) < 0.5F)
              {
                root_hi = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::initial_guess(vv_floor, m, pol);
              }
              else
              {
                root_hi = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess(T(vv_floor + 0.5F), m, pol);
              }

              root_lo = T(root_hi - 0.1F);

              const bool hi_end_of_bracket_is_negative = (boost::math::cyl_neumann(v, root_hi, pol) < 0);

              while((root_lo > boost::math::tools::epsilon<T>()))
              {
                const bool lo_end_of_bracket_is_negative = (boost::math::cyl_neumann(v, root_lo, pol) < 0);

                if(hi_end_of_bracket_is_negative != lo_end_of_bracket_is_negative)
                {
                  break;
                }

                root_hi = root_lo;

                // Decrease the lower end of the bracket using an adaptive algorithm.
                if(root_lo > 0.5F)
                {
                  root_lo -= 0.5F;
                }
                else
                {
                  root_lo *= 0.75F;
                }
              }
            }
            else
            {
              if(T(vv - vv_floor) < 0.5F)
              {
                root_lo  = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::initial_guess(vv_floor, m - 1, pol);
                root_hi = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::initial_guess(vv_floor, m, pol);
                root_lo += 0.01F;
                root_hi += 0.01F;
              }
              else
              {
                root_lo = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess(T(vv_floor + 0.5F), m - 1, pol);
                root_hi = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess(T(vv_floor + 0.5F), m, pol);
                root_lo += 0.01F;
                root_hi += 0.01F;
              }
            }

            // Perform several steps of bisection iteration to refine the guess.
            boost::uintmax_t number_of_iterations(12U);

            // Do the bisection iteration.
            const boost::math::tuple<T, T> guess_pair =
               boost::math::tools::bisect(
                  boost::math::detail::bessel_zero::cyl_neumann_zero_detail::function_object_yv<T, Policy>(v, pol),
                  root_lo,
                  root_hi,
                  boost::math::detail::bessel_zero::cyl_neumann_zero_detail::my_bisection_unreachable_tolerance<T>,
                  number_of_iterations);

            return (boost::math::get<0>(guess_pair) + boost::math::get<1>(guess_pair)) / 2U;
          }

          if(m == 1U)
          {
            // Get the initial estimate of the first root.

            if(v < 2.2F)
            {
              // For small v, use the results of empirical curve fitting.
              // Mathematica(R) session for the coefficients:
              //  Table[{n, BesselYZero[n, 1]}, {n, 0, 22/10, 1/10}]
              //  N[%, 20]
              //  Fit[%, {n^0, n^1, n^2, n^3, n^4, n^5, n^6}, n]
              guess = (((((    - T(0.0025095909235652)
                           * v + T(0.021291887049053))
                           * v - T(0.076487785486526))
                           * v + T(0.159110268115362))
                           * v - T(0.241681668765196))
                           * v + T(1.4437846310885244))
                           * v + T(0.89362115190200490);
            }
            else
            {
              // For larger v, use the second line of Eqs. 10.21.40 in the NIST Handbook.
              guess = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::equation_nist_10_21_40_b(v);
            }
          }
          else
          {
            if(v < 2.2F)
            {
              // Use Eq. 10.21.19 in the NIST Handbook.
              const T a(((v + T(m * 2U)) - T(1.5)) * boost::math::constants::half_pi<T>());

              guess = boost::math::detail::bessel_zero::equation_nist_10_21_19(v, a);
            }
            else
            {
              // Get an estimate of the m'th root of airy_bi.
              const T airy_bi_root(boost::math::detail::airy_zero::airy_bi_zero_detail::initial_guess<T>(m));

              // Use Eq. 9.5.26 in the A&S Handbook.
              guess = boost::math::detail::bessel_zero::equation_as_9_5_26(v, airy_bi_root);
            }
          }

          return guess;
        }
      } // namespace cyl_neumann_zero_detail
    } // namespace bessel_zero
  } } } // namespace boost::math::detail

#endif // _BESSEL_JY_ZERO_2013_01_18_HPP_


code:
//  Copyright (c) 2006 Xiaogang Zhang
//  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_BESSEL_K0_HPP
#define BOOST_MATH_BESSEL_K0_HPP

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

#include <boost/math/tools/rational.hpp>
#include <boost/math/tools/big_constant.hpp>
#include <boost/math/policies/error_handling.hpp>
#include <boost/assert.hpp>

// Modified Bessel function of the second kind of order zero
// minimax rational approximations on intervals, see
// Russon and Blair, Chalk River Report AECL-3461, 1969

namespace boost { namespace math { namespace detail{

template <typename T, typename Policy>
T bessel_k0(T x, const Policy&);

template <class T, class Policy>
struct bessel_k0_initializer
{
   struct init
   {
      init()
      {
         do_init();
      }
      static void do_init()
      {
         bessel_k0(T(1), Policy());
      }
      void force_instantiate()const{}
   };
   static const init initializer;
   static void force_instantiate()
   {
      initializer.force_instantiate();
   }
};

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

template <typename T, typename Policy>
T bessel_k0(T x, const Policy& pol)
{
    BOOST_MATH_INSTRUMENT_CODE(x);

    bessel_k0_initializer<T, Policy>::force_instantiate();

    static const T P1[] = {
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.4708152720399552679e+03)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.9169059852270512312e+03)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.6850901201934832188e+02)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.1999463724910714109e+01)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.3166052564989571850e-01)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.8599221412826100000e-04))
    };
    static const T Q1[] = {
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.1312714303849120380e+04)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.4994418972832303646e+02)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0))
    };
    static const T P2[] = {
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.6128136304458193998e+06)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.7333769444840079748e+05)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.7984434409411765813e+04)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.9501657892958843865e+02)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.6414452837299064100e+00))
    };
    static const T Q2[] = {
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.6128136304458193998e+06)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.9865713163054025489e+04)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.5064972445877992730e+02)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0))
    };
    static const T P3[] = {
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.1600249425076035558e+02)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.3444738764199315021e+03)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.8321525870183537725e+04)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 7.1557062783764037541e+04)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.5097646353289914539e+05)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.7398867902565686251e+05)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0577068948034021957e+05)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.1075408980684392399e+04)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.6832589957340267940e+03)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.1394980557384778174e+02))
    };
    static const T Q3[] = {
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 9.2556599177304839811e+01)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.8821890840982713696e+03)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.4847228371802360957e+04)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.8824616785857027752e+04)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.2689839587977598727e+05)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.5144644673520157801e+05)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 9.7418829762268075784e+04)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.1474655750295278825e+04)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.4329628889746408858e+03)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.0013443064949242491e+02)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0))
    };
    T value, factor, r, r1, r2;

    BOOST_MATH_STD_USING
    using namespace boost::math::tools;

    static const char* function = "boost::math::bessel_k0<%1%>(%1%,%1%)";

    if (x < 0)
    {
       return policies::raise_domain_error<T>(function,
            "Got x = %1%, but argument x must be non-negative, complex number result not supported", x, pol);
    }
    if (x == 0)
    {
       return policies::raise_overflow_error<T>(function, 0, pol);
    }
    if (x <= 1)                         // x in (0, 1]
    {
        T y = x * x;
        r1 = evaluate_polynomial(P1, y) / evaluate_polynomial(Q1, y);
        r2 = evaluate_polynomial(P2, y) / evaluate_polynomial(Q2, y);
        factor = log(x);
        value = r1 - factor * r2;
    }
    else                                // x in (1, \infty)
    {
        T y = 1 / x;
        r = evaluate_polynomial(P3, y) / evaluate_polynomial(Q3, y);
        factor = exp(-x) / sqrt(x);
        value = factor * r;
        BOOST_MATH_INSTRUMENT_CODE("y = " << y);
        BOOST_MATH_INSTRUMENT_CODE("r = " << r);
        BOOST_MATH_INSTRUMENT_CODE("factor = " << factor);
        BOOST_MATH_INSTRUMENT_CODE("value = " << value);
    }

    return value;
}

}}} // namespaces

#ifdef _MSC_VER
#pragma warning(pop)
#endif

#endif // BOOST_MATH_BESSEL_K0_HPP

license: mit


code:
//  Copyright (c) 2006 Xiaogang Zhang
//  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_BESSEL_K1_HPP
#define BOOST_MATH_BESSEL_K1_HPP

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

#include <boost/math/tools/rational.hpp>
#include <boost/math/tools/big_constant.hpp>
#include <boost/math/policies/error_handling.hpp>
#include <boost/assert.hpp>

// Modified Bessel function of the second kind of order one
// minimax rational approximations on intervals, see
// Russon and Blair, Chalk River Report AECL-3461, 1969

namespace boost { namespace math { namespace detail{

template <typename T, typename Policy>
T bessel_k1(T x, const Policy&);

template <class T, class Policy>
struct bessel_k1_initializer
{
   struct init
   {
      init()
      {
         do_init();
      }
      static void do_init()
      {
         bessel_k1(T(1), Policy());
      }
      void force_instantiate()const{}
   };
   static const init initializer;
   static void force_instantiate()
   {
      initializer.force_instantiate();
   }
};

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

template <typename T, typename Policy>
T bessel_k1(T x, const Policy& pol)
{
    bessel_k1_initializer<T, Policy>::force_instantiate();

    static const T P1[] = {
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.2149374878243304548e+06)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 7.1938920065420586101e+05)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.7733324035147015630e+05)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 7.1885382604084798576e+03)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 9.9991373567429309922e+01)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.8127070456878442310e-01))
    };
    static const T Q1[] = {
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.2149374878243304548e+06)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.7264298672067697862e+04)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.8143915754538725829e+02)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0))
    };
    static const T P2[] = {
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.0)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.3531161492785421328e+06)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.4758069205414222471e+05)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -4.5051623763436087023e+03)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -5.3103913335180275253e+01)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.2795590826955002390e-01))
    };
    static const T Q2[] = {
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.7062322985570842656e+06)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.3117653211351080007e+04)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.0507151578787595807e+02)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0))
    };
    static const T P3[] = {
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.2196792496874548962e+00)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.4137176114230414036e+01)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.4122953486801312910e+02)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.3319486433183221990e+03)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.8590657697910288226e+03)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.4540675585544584407e+03)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.3123742209168871550e+03)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.1094256146537402173e+02)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.3182609918569941308e+02)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 7.5584584631176030810e+00)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 6.4257745859173138767e-02))
    };
    static const T Q3[] = {
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.7710478032601086579e+00)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.4552228452758912848e+01)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.5951223655579051357e+02)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 9.6929165726802648634e+02)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.9448440788918006154e+03)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.1181000487171943810e+03)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.2082692316002348638e+03)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.3031020088765390854e+02)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.6001069306861518855e+01)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0))
    };
    T value, factor, r, r1, r2;

    BOOST_MATH_STD_USING
    using namespace boost::math::tools;

    static const char* function = "boost::math::bessel_k1<%1%>(%1%,%1%)";

    if (x < 0)
    {
       return policies::raise_domain_error<T>(function,
            "Got x = %1%, but argument x must be non-negative, complex number result not supported.", x, pol);
    }
    if (x == 0)
    {
       return policies::raise_overflow_error<T>(function, 0, pol);
    }
    if (x <= 1)                         // x in (0, 1]
    {
        T y = x * x;
        r1 = evaluate_polynomial(P1, y) /  evaluate_polynomial(Q1, y);
        r2 = evaluate_polynomial(P2, y) /  evaluate_polynomial(Q2, y);
        factor = log(x);
        value = (r1 + factor * r2) / x;
    }
    else                                // x in (1, \infty)
    {
        T y = 1 / x;
        r = evaluate_polynomial(P3, y) /  evaluate_polynomial(Q3, y);
        factor = exp(-x) / sqrt(x);
        value = factor * r;
    }

    return value;
}

}}} // namespaces

#ifdef _MSC_VER
#pragma warning(pop)
#endif

#endif // BOOST_MATH_BESSEL_K1_HPP


code:
//  Copyright (c) 2006 Xiaogang Zhang
//  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_BESSEL_KN_HPP
#define BOOST_MATH_BESSEL_KN_HPP

#ifdef _MSC_VER
#pragma once
#endif

#include <boost/math/special_functions/detail/bessel_k0.hpp>
#include <boost/math/special_functions/detail/bessel_k1.hpp>
#include <boost/math/policies/error_handling.hpp>

// Modified Bessel function of the second kind of integer order
// K_n(z) is the dominant solution, forward recurrence always OK (though unstable)

namespace boost { namespace math { namespace detail{

template <typename T, typename Policy>
T bessel_kn(int n, T x, const Policy& pol)
{
    BOOST_MATH_STD_USING
    T value, current, prev;

    using namespace boost::math::tools;

    static const char* function = "boost::math::bessel_kn<%1%>(%1%,%1%)";

    if (x < 0)
    {
       return policies::raise_domain_error<T>(function,
            "Got x = %1%, but argument x must be non-negative, complex number result not supported.", x, pol);
    }
    if (x == 0)
    {
       return policies::raise_overflow_error<T>(function, 0, pol);
    }

    if (n < 0)
    {
        n = -n;                             // K_{-n}(z) = K_n(z)
    }
    if (n == 0)
    {
        value = bessel_k0(x, pol);
    }
    else if (n == 1)
    {
        value = bessel_k1(x, pol);
    }
    else
    {
       prev = bessel_k0(x, pol);
       current = bessel_k1(x, pol);
       int k = 1;
       BOOST_ASSERT(k < n);
       T scale = 1;
       do
       {
           T fact = 2 * k / x;
           if((tools::max_value<T>() - fabs(prev)) / fact < fabs(current))
           {
              scale /= current;
              prev /= current;
              current = 1;
           }
           value = fact * current + prev;
           prev = current;
           current = value;
           ++k;
       }
       while(k < n);
       if(tools::max_value<T>() * scale < fabs(value))
          return sign(scale) * sign(value) * policies::raise_overflow_error<T>(function, 0, pol);
       value /= scale;
    }
    return value;
}

}}} // namespaces

#endif // BOOST_MATH_BESSEL_KN_HPP


code:
//  Copyright (c) 2006 Xiaogang Zhang
//  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_BESSEL_Y0_HPP
#define BOOST_MATH_BESSEL_Y0_HPP

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

#include <boost/math/special_functions/detail/bessel_j0.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/tools/rational.hpp>
#include <boost/math/tools/big_constant.hpp>
#include <boost/math/policies/error_handling.hpp>
#include <boost/assert.hpp>

// Bessel function of the second kind of order zero
// x <= 8, minimax rational approximations on root-bracketing intervals
// x > 8, Hankel asymptotic expansion in Hart, Computer Approximations, 1968

namespace boost { namespace math { namespace detail{

template <typename T, typename Policy>
T bessel_y0(T x, const Policy&);

template <class T, class Policy>
struct bessel_y0_initializer
{
   struct init
   {
      init()
      {
         do_init();
      }
      static void do_init()
      {
         bessel_y0(T(1), Policy());
      }
      void force_instantiate()const{}
   };
   static const init initializer;
   static void force_instantiate()
   {
      initializer.force_instantiate();
   }
};

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

template <typename T, typename Policy>
T bessel_y0(T x, const Policy& pol)
{
    bessel_y0_initializer<T, Policy>::force_instantiate();

    static const T P1[] = {
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0723538782003176831e+11)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -8.3716255451260504098e+09)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.0422274357376619816e+08)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.1287548474401797963e+06)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0102532948020907590e+04)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.8402381979244993524e+01)),
    };
    static const T Q1[] = {
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.8873865738997033405e+11)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.1617187777290363573e+09)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.5662956624278251596e+07)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.3889393209447253406e+05)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 6.6475986689240190091e+02)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)),
    };
    static const T P2[] = {
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.2213976967566192242e+13)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -5.5107435206722644429e+11)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.3600098638603061642e+10)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -6.9590439394619619534e+08)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.6905288611678631510e+06)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.4566865832663635920e+04)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.7427031242901594547e+01)),
    };
    static const T Q2[] = {
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.3386146580707264428e+14)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.4266824419412347550e+12)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.4015103849971240096e+10)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.3960202770986831075e+08)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.0669982352539552018e+05)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.3030857612070288823e+02)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)),
    };
    static const T P3[] = {
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -8.0728726905150210443e+15)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 6.7016641869173237784e+14)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.2829912364088687306e+11)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.9363051266772083678e+11)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.1958827170518100757e+09)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.0085539923498211426e+07)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.1363534169313901632e+04)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.7439661319197499338e+01)),
    };
    static const T Q3[] = {
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.4563724628846457519e+17)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.9272425569640309819e+15)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.2598377924042897629e+13)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.6926121104209825246e+10)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.4727219475672302327e+08)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.3924739209768057030e+05)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.7903362168128450017e+02)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)),
    };
    static const T PC[] = {
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.2779090197304684302e+04)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.1345386639580765797e+04)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.1170523380864944322e+04)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.4806486443249270347e+03)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.5376201909008354296e+02)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.8961548424210455236e-01)),
    };
    static const T QC[] = {
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.2779090197304684318e+04)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.1370412495510416640e+04)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.1215350561880115730e+04)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.5028735138235608207e+03)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.5711159858080893649e+02)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)),
    };
    static const T PS[] = {
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -8.9226600200800094098e+01)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.8591953644342993800e+02)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.1183429920482737611e+02)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.2300261666214198472e+01)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.2441026745835638459e+00)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -8.8033303048680751817e-03)),
    };
    static const T QS[] = {
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.7105024128512061905e+03)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.1951131543434613647e+04)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 7.2642780169211018836e+03)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.4887231232283756582e+03)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 9.0593769594993125859e+01)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)),
    };
    static const T x1  =  static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.9357696627916752158e-01)),
                   x2  =  static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.9576784193148578684e+00)),
                   x3  =  static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 7.0860510603017726976e+00)),
                   x11 =  static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.280e+02)),
                   x12 =  static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.9519662791675215849e-03)),
                   x21 =  static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0130e+03)),
                   x22 =  static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 6.4716931485786837568e-04)),
                   x31 =  static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.8140e+03)),
                   x32 =  static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.1356030177269762362e-04))
    ;
    T value, factor, r, rc, rs;

    BOOST_MATH_STD_USING
    using namespace boost::math::tools;
    using namespace boost::math::constants;

    static const char* function = "boost::math::bessel_y0<%1%>(%1%,%1%)";

    if (x < 0)
    {
       return policies::raise_domain_error<T>(function,
            "Got x = %1% but x must be non-negative, complex result not supported.", x, pol);
    }
    if (x == 0)
    {
       return -policies::raise_overflow_error<T>(function, 0, pol);
    }
    if (x <= 3)                       // x in (0, 3]
    {
        T y = x * x;
        T z = 2 * log(x/x1) * bessel_j0(x) / pi<T>();
        r = evaluate_rational(P1, Q1, y);
        factor = (x + x1) * ((x - x11/256) - x12);
        value = z + factor * r;
    }
    else if (x <= 5.5f)                  // x in (3, 5.5]
    {
        T y = x * x;
        T z = 2 * log(x/x2) * bessel_j0(x) / pi<T>();
        r = evaluate_rational(P2, Q2, y);
        factor = (x + x2) * ((x - x21/256) - x22);
        value = z + factor * r;
    }
    else if (x <= 8)                  // x in (5.5, 8]
    {
        T y = x * x;
        T z = 2 * log(x/x3) * bessel_j0(x) / pi<T>();
        r = evaluate_rational(P3, Q3, y);
        factor = (x + x3) * ((x - x31/256) - x32);
        value = z + factor * r;
    }
    else                                // x in (8, \infty)
    {
        T y = 8 / x;
        T y2 = y * y;
        rc = evaluate_rational(PC, QC, y2);
        rs = evaluate_rational(PS, QS, y2);
        factor = constants::one_div_root_pi<T>() / sqrt(x);
        //
        // The following code is really just:
        //
        // T z = x - 0.25f * pi<T>();
        // value = factor * (rc * sin(z) + y * rs * cos(z));
        //
        // But using the sin/cos addition formulae and constant values for
        // sin/cos of PI/4 which then cancel part of the "factor" term as they're all
        // 1 / sqrt(2):
        //
        T sx = sin(x);
        T cx = cos(x);
        value = factor * (rc * (sx - cx) + y * rs * (cx + sx));
    }

    return value;
}

}}} // namespaces

#ifdef _MSC_VER
#pragma warning(pop)
#endif

#endif // BOOST_MATH_BESSEL_Y0_HPP


code:
//  Copyright (c) 2006 Xiaogang Zhang
//  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_BESSEL_Y1_HPP
#define BOOST_MATH_BESSEL_Y1_HPP

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

#include <boost/math/special_functions/detail/bessel_j1.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/tools/rational.hpp>
#include <boost/math/tools/big_constant.hpp>
#include <boost/math/policies/error_handling.hpp>
#include <boost/assert.hpp>

// Bessel function of the second kind of order one
// x <= 8, minimax rational approximations on root-bracketing intervals
// x > 8, Hankel asymptotic expansion in Hart, Computer Approximations, 1968

namespace boost { namespace math { namespace detail{

template <typename T, typename Policy>
T bessel_y1(T x, const Policy&);

template <class T, class Policy>
struct bessel_y1_initializer
{
   struct init
   {
      init()
      {
         do_init();
      }
      static void do_init()
      {
         bessel_y1(T(1), Policy());
      }
      void force_instantiate()const{}
   };
   static const init initializer;
   static void force_instantiate()
   {
      initializer.force_instantiate();
   }
};

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

template <typename T, typename Policy>
T bessel_y1(T x, const Policy& pol)
{
    bessel_y1_initializer<T, Policy>::force_instantiate();

    static const T P1[] = {
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.0535726612579544093e+13)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.4708611716525426053e+12)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.7595974497819597599e+11)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 7.2144548214502560419e+09)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -5.9157479997408395984e+07)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.2157953222280260820e+05)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.1714424660046133456e+02)),
    };
    static const T Q1[] = {
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.0737873921079286084e+14)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.1272286200406461981e+12)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.7800352738690585613e+10)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.2250435122182963220e+08)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.8136470753052572164e+05)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.2079908168393867438e+02)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)),
    };
    static const T P2[] = {
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.1514276357909013326e+19)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -5.6808094574724204577e+18)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.3638408497043134724e+16)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.0686275289804744814e+15)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -5.9530713129741981618e+13)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.7453673962438488783e+11)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.1957961912070617006e+09)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.9153806858264202986e+06)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.2337180442012953128e+03)),
    };
    static const T Q2[] = {
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.3321844313316185697e+20)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.6968198822857178911e+18)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.0837179548112881950e+16)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.1187010065856971027e+14)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.0221766852960403645e+11)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 6.3550318087088919566e+08)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0453748201934079734e+06)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.2855164849321609336e+03)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)),
    };
    static const T PC[] = {
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -4.4357578167941278571e+06)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -9.9422465050776411957e+06)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -6.6033732483649391093e+06)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.5235293511811373833e+06)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.0982405543459346727e+05)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.6116166443246101165e+03)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.0)),
    };
    static const T QC[] = {
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -4.4357578167941278568e+06)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -9.9341243899345856590e+06)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -6.5853394797230870728e+06)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.5118095066341608816e+06)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.0726385991103820119e+05)),
        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.4550094401904961825e+03)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)),
    };
    static const T PS[] = {
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.3220913409857223519e+04)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.5145160675335701966e+04)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 6.6178836581270835179e+04)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.8494262873223866797e+04)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.7063754290207680021e+03)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.5265133846636032186e+01)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.0)),
    };
    static const T QS[] = {
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 7.0871281941028743574e+05)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.8194580422439972989e+06)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.4194606696037208929e+06)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.0029443582266975117e+05)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.7890229745772202641e+04)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.6383677696049909675e+02)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)),
    };
    static const T x1  =  static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.1971413260310170351e+00)),
                   x2  =  static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.4296810407941351328e+00)),
                   x11 =  static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.620e+02)),
                   x12 =  static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.8288260310170351490e-03)),
                   x21 =  static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.3900e+03)),
                   x22 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -6.4592058648672279948e-06))
    ;
    T value, factor, r, rc, rs;

    BOOST_MATH_STD_USING
    using namespace boost::math::tools;
    using namespace boost::math::constants;

    if (x <= 0)
    {
       return policies::raise_domain_error<T>("bost::math::bessel_y1<%1%>(%1%,%1%)",
            "Got x == %1%, but x must be > 0, complex result not supported.", x, pol);
    }
    if (x <= 4)                       // x in (0, 4]
    {
        T y = x * x;
        T z = 2 * log(x/x1) * bessel_j1(x) / pi<T>();
        r = evaluate_rational(P1, Q1, y);
        factor = (x + x1) * ((x - x11/256) - x12) / x;
        value = z + factor * r;
    }
    else if (x <= 8)                  // x in (4, 8]
    {
        T y = x * x;
        T z = 2 * log(x/x2) * bessel_j1(x) / pi<T>();
        r = evaluate_rational(P2, Q2, y);
        factor = (x + x2) * ((x - x21/256) - x22) / x;
        value = z + factor * r;
    }
    else                                // x in (8, \infty)
    {
        T y = 8 / x;
        T y2 = y * y;
        rc = evaluate_rational(PC, QC, y2);
        rs = evaluate_rational(PS, QS, y2);
        factor = 1 / (sqrt(x) * root_pi<T>());
        //
        // This code is really just:
        //
        // T z = x - 0.75f * pi<T>();
        // value = factor * (rc * sin(z) + y * rs * cos(z));
        //
        // But using the sin/cos addition rules, plus constants for sin/cos of 3PI/4
        // which then cancel out with corresponding terms in "factor".
        //
        T sx = sin(x);
        T cx = cos(x);
        value = factor * (y * rs * (sx - cx) - rc * (sx + cx));
    }

    return value;
}

}}} // namespaces

#ifdef _MSC_VER
#pragma warning(pop)
#endif

#endif // BOOST_MATH_BESSEL_Y1_HPP


code:
//  Copyright (c) 2006 Xiaogang Zhang
//  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_BESSEL_YN_HPP
#define BOOST_MATH_BESSEL_YN_HPP

#ifdef _MSC_VER
#pragma once
#endif

#include <boost/math/special_functions/detail/bessel_y0.hpp>
#include <boost/math/special_functions/detail/bessel_y1.hpp>
#include <boost/math/special_functions/detail/bessel_jy_series.hpp>
#include <boost/math/policies/error_handling.hpp>

// Bessel function of the second kind of integer order
// Y_n(z) is the dominant solution, forward recurrence always OK (though unstable)

namespace boost { namespace math { namespace detail{

template <typename T, typename Policy>
T bessel_yn(int n, T x, const Policy& pol)
{
    BOOST_MATH_STD_USING
    T value, factor, current, prev;

    using namespace boost::math::tools;

    static const char* function = "boost::math::bessel_yn<%1%>(%1%,%1%)";

    if ((x == 0) && (n == 0))
    {
       return -policies::raise_overflow_error<T>(function, 0, pol);
    }
    if (x <= 0)
    {
       return policies::raise_domain_error<T>(function,
            "Got x = %1%, but x must be > 0, complex result not supported.", x, pol);
    }

    //
    // Reflection comes first:
    //
    if (n < 0)
    {
        factor = static_cast<T>((n & 0x1) ? -1 : 1);  // Y_{-n}(z) = (-1)^n Y_n(z)
        n = -n;
    }
    else
    {
        factor = 1;
    }
    if(x < policies::get_epsilon<T, Policy>())
    {
       T scale = 1;
       value = bessel_yn_small_z(n, x, &scale, pol);
       if(tools::max_value<T>() * fabs(scale) < fabs(value))
          return boost::math::sign(scale) * boost::math::sign(value) * policies::raise_overflow_error<T>(function, 0, pol);
       value /= scale;
    }
    else if(asymptotic_bessel_large_x_limit(n, x))
    {
       value = factor * asymptotic_bessel_y_large_x_2(static_cast<T>(abs(n)), x);
    }
    else if (n == 0)
    {
        value = bessel_y0(x, pol);
    }
    else if (n == 1)
    {
        value = factor * bessel_y1(x, pol);
    }
    else
    {
       prev = bessel_y0(x, pol);
       current = bessel_y1(x, pol);
       int k = 1;
       BOOST_ASSERT(k < n);
       policies::check_series_iterations<T>("boost::math::bessel_y_n<%1%>(%1%,%1%)", n, pol);
       T mult = 2 * k / x;
       value = mult * current - prev;
       prev = current;
       current = value;
       ++k;
       if((mult > 1) && (fabs(current) > 1))
       {
          prev /= current;
          factor /= current;
          value /= current;
          current = 1;
       }
       while(k < n)
       {
           mult = 2 * k / x;
           value = mult * current - prev;
           prev = current;
           current = value;
           ++k;
       }
       if(fabs(tools::max_value<T>() * factor) < fabs(value))
          return sign(value) * sign(factor) * policies::raise_overflow_error<T>(function, 0, pol);
       value /= factor;
    }
    return value;
}

}}} // namespaces

#endif // BOOST_MATH_BESSEL_YN_HPP


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_SF_ERF_INV_HPP
#define BOOST_MATH_SF_ERF_INV_HPP

#ifdef _MSC_VER
#pragma once
#pragma warning(push)
#pragma warning(disable:4127) // Conditional expression is constant
#pragma warning(disable:4702) // Unreachable code: optimization warning
#endif

namespace boost{ namespace math{

namespace detail{
//
// The inverse erf and erfc functions share a common implementation,
// this version is for 80-bit long double's and smaller:
//
template <class T, class Policy>
T erf_inv_imp(const T& p, const T& q, const Policy&, const boost::mpl::int_<64>*)
{
   BOOST_MATH_STD_USING // for ADL of std names.

   T result = 0;

   if(p <= 0.5)
   {
      //
      // Evaluate inverse erf using the rational approximation:
      //
      // x = p(p+10)(Y+R(p))
      //
      // Where Y is a constant, and R(p) is optimised for a low
      // absolute error compared to |Y|.
      //
      // double: Max error found: 2.001849e-18
      // long double: Max error found: 1.017064e-20
      // Maximum Deviation Found (actual error term at infinite precision) 8.030e-21
      //
      static const float Y = 0.0891314744949340820313f;
      static const T P[] = {
         BOOST_MATH_BIG_CONSTANT(T, 64, -0.000508781949658280665617),
         BOOST_MATH_BIG_CONSTANT(T, 64, -0.00836874819741736770379),
         BOOST_MATH_BIG_CONSTANT(T, 64, 0.0334806625409744615033),
         BOOST_MATH_BIG_CONSTANT(T, 64, -0.0126926147662974029034),
         BOOST_MATH_BIG_CONSTANT(T, 64, -0.0365637971411762664006),
         BOOST_MATH_BIG_CONSTANT(T, 64, 0.0219878681111168899165),
         BOOST_MATH_BIG_CONSTANT(T, 64, 0.00822687874676915743155),
         BOOST_MATH_BIG_CONSTANT(T, 64, -0.00538772965071242932965)
      };
      static const T Q[] = {
         BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
         BOOST_MATH_BIG_CONSTANT(T, 64, -0.970005043303290640362),
         BOOST_MATH_BIG_CONSTANT(T, 64, -1.56574558234175846809),
         BOOST_MATH_BIG_CONSTANT(T, 64, 1.56221558398423026363),
         BOOST_MATH_BIG_CONSTANT(T, 64, 0.662328840472002992063),
         BOOST_MATH_BIG_CONSTANT(T, 64, -0.71228902341542847553),
         BOOST_MATH_BIG_CONSTANT(T, 64, -0.0527396382340099713954),
         BOOST_MATH_BIG_CONSTANT(T, 64, 0.0795283687341571680018),
         BOOST_MATH_BIG_CONSTANT(T, 64, -0.00233393759374190016776),
         BOOST_MATH_BIG_CONSTANT(T, 64, 0.000886216390456424707504)
      };
      T g = p * (p + 10);
      T r = tools::evaluate_polynomial(P, p) / tools::evaluate_polynomial(Q, p);
      result = g * Y + g * r;
   }
   else if(q >= 0.25)
   {
      //
      // Rational approximation for 0.5 > q >= 0.25
      //
      // x = sqrt(-2*log(q)) / (Y + R(q))
      //
      // Where Y is a constant, and R(q) is optimised for a low
      // absolute error compared to Y.
      //
      // double : Max error found: 7.403372e-17
      // long double : Max error found: 6.084616e-20
      // Maximum Deviation Found (error term) 4.811e-20
      //
      static const float Y = 2.249481201171875f;
      static const T P[] = {
         BOOST_MATH_BIG_CONSTANT(T, 64, -0.202433508355938759655),
         BOOST_MATH_BIG_CONSTANT(T, 64, 0.105264680699391713268),
         BOOST_MATH_BIG_CONSTANT(T, 64, 8.37050328343119927838),
         BOOST_MATH_BIG_CONSTANT(T, 64, 17.6447298408374015486),
         BOOST_MATH_BIG_CONSTANT(T, 64, -18.8510648058714251895),
         BOOST_MATH_BIG_CONSTANT(T, 64, -44.6382324441786960818),
         BOOST_MATH_BIG_CONSTANT(T, 64, 17.445385985570866523),
         BOOST_MATH_BIG_CONSTANT(T, 64, 21.1294655448340526258),
         BOOST_MATH_BIG_CONSTANT(T, 64, -3.67192254707729348546)
      };
      static const T Q[] = {
         BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
         BOOST_MATH_BIG_CONSTANT(T, 64, 6.24264124854247537712),
         BOOST_MATH_BIG_CONSTANT(T, 64, 3.9713437953343869095),
         BOOST_MATH_BIG_CONSTANT(T, 64, -28.6608180499800029974),
         BOOST_MATH_BIG_CONSTANT(T, 64, -20.1432634680485188801),
         BOOST_MATH_BIG_CONSTANT(T, 64, 48.5609213108739935468),
         BOOST_MATH_BIG_CONSTANT(T, 64, 10.8268667355460159008),
         BOOST_MATH_BIG_CONSTANT(T, 64, -22.6436933413139721736),
         BOOST_MATH_BIG_CONSTANT(T, 64, 1.72114765761200282724)
      };
      T g = sqrt(-2 * log(q));
      T xs = q - 0.25f;
      T r = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
      result = g / (Y + r);
   }
   else
   {
      //
      // For q < 0.25 we have a series of rational approximations all
      // of the general form:
      //
      // let: x = sqrt(-log(q))
      //
      // Then the result is given by:
      //
      // x(Y+R(x-B))
      //
      // where Y is a constant, B is the lowest value of x for which
      // the approximation is valid, and R(x-B) is optimised for a low
      // absolute error compared to Y.
      //
      // Note that almost all code will really go through the first
      // or maybe second approximation.  After than we're dealing with very
      // small input values indeed: 80 and 128 bit long double's go all the
      // way down to ~ 1e-5000 so the "tail" is rather long...
      //
      T x = sqrt(-log(q));
      if(x < 3)
      {
         // Max error found: 1.089051e-20
         static const float Y = 0.807220458984375f;
         static const T P[] = {
            BOOST_MATH_BIG_CONSTANT(T, 64, -0.131102781679951906451),
            BOOST_MATH_BIG_CONSTANT(T, 64, -0.163794047193317060787),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.117030156341995252019),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.387079738972604337464),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.337785538912035898924),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.142869534408157156766),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.0290157910005329060432),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.00214558995388805277169),
            BOOST_MATH_BIG_CONSTANT(T, 64, -0.679465575181126350155e-6),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.285225331782217055858e-7),
            BOOST_MATH_BIG_CONSTANT(T, 64, -0.681149956853776992068e-9)
         };
         static const T Q[] = {
            BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
            BOOST_MATH_BIG_CONSTANT(T, 64, 3.46625407242567245975),
            BOOST_MATH_BIG_CONSTANT(T, 64, 5.38168345707006855425),
            BOOST_MATH_BIG_CONSTANT(T, 64, 4.77846592945843778382),
            BOOST_MATH_BIG_CONSTANT(T, 64, 2.59301921623620271374),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.848854343457902036425),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.152264338295331783612),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.01105924229346489121)
         };
         T xs = x - 1.125f;
         T R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
         result = Y * x + R * x;
      }
      else if(x < 6)
      {
         // Max error found: 8.389174e-21
         static const float Y = 0.93995571136474609375f;
         static const T P[] = {
            BOOST_MATH_BIG_CONSTANT(T, 64, -0.0350353787183177984712),
            BOOST_MATH_BIG_CONSTANT(T, 64, -0.00222426529213447927281),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.0185573306514231072324),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.00950804701325919603619),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.00187123492819559223345),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.000157544617424960554631),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.460469890584317994083e-5),
            BOOST_MATH_BIG_CONSTANT(T, 64, -0.230404776911882601748e-9),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.266339227425782031962e-11)
         };
         static const T Q[] = {
            BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
            BOOST_MATH_BIG_CONSTANT(T, 64, 1.3653349817554063097),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.762059164553623404043),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.220091105764131249824),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.0341589143670947727934),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.00263861676657015992959),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.764675292302794483503e-4)
         };
         T xs = x - 3;
         T R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
         result = Y * x + R * x;
      }
      else if(x < 18)
      {
         // Max error found: 1.481312e-19
         static const float Y = 0.98362827301025390625f;
         static const T P[] = {
            BOOST_MATH_BIG_CONSTANT(T, 64, -0.0167431005076633737133),
            BOOST_MATH_BIG_CONSTANT(T, 64, -0.00112951438745580278863),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.00105628862152492910091),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.000209386317487588078668),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.149624783758342370182e-4),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.449696789927706453732e-6),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.462596163522878599135e-8),
            BOOST_MATH_BIG_CONSTANT(T, 64, -0.281128735628831791805e-13),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.99055709973310326855e-16)
         };
         static const T Q[] = {
            BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.591429344886417493481),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.138151865749083321638),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.0160746087093676504695),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.000964011807005165528527),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.275335474764726041141e-4),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.282243172016108031869e-6)
         };
         T xs = x - 6;
         T R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
         result = Y * x + R * x;
      }
      else if(x < 44)
      {
         // Max error found: 5.697761e-20
         static const float Y = 0.99714565277099609375f;
         static const T P[] = {
            BOOST_MATH_BIG_CONSTANT(T, 64, -0.0024978212791898131227),
            BOOST_MATH_BIG_CONSTANT(T, 64, -0.779190719229053954292e-5),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.254723037413027451751e-4),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.162397777342510920873e-5),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.396341011304801168516e-7),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.411632831190944208473e-9),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.145596286718675035587e-11),
            BOOST_MATH_BIG_CONSTANT(T, 64, -0.116765012397184275695e-17)
         };
         static const T Q[] = {
            BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.207123112214422517181),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.0169410838120975906478),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.000690538265622684595676),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.145007359818232637924e-4),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.144437756628144157666e-6),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.509761276599778486139e-9)
         };
         T xs = x - 18;
         T R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
         result = Y * x + R * x;
      }
      else
      {
         // Max error found: 1.279746e-20
         static const float Y = 0.99941349029541015625f;
         static const T P[] = {
            BOOST_MATH_BIG_CONSTANT(T, 64, -0.000539042911019078575891),
            BOOST_MATH_BIG_CONSTANT(T, 64, -0.28398759004727721098e-6),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.899465114892291446442e-6),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.229345859265920864296e-7),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.225561444863500149219e-9),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.947846627503022684216e-12),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.135880130108924861008e-14),
            BOOST_MATH_BIG_CONSTANT(T, 64, -0.348890393399948882918e-21)
         };
         static const T Q[] = {
            BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.0845746234001899436914),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.00282092984726264681981),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.468292921940894236786e-4),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.399968812193862100054e-6),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.161809290887904476097e-8),
            BOOST_MATH_BIG_CONSTANT(T, 64, 0.231558608310259605225e-11)
         };
         T xs = x - 44;
         T R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
         result = Y * x + R * x;
      }
   }
   return result;
}

template <class T, class Policy>
struct erf_roots
{
   boost::math::tuple<T,T,T> operator()(const T& guess)
   {
      BOOST_MATH_STD_USING
      T derivative = sign * (2 / sqrt(constants::pi<T>())) * exp(-(guess * guess));
      T derivative2 = -2 * guess * derivative;
      return boost::math::make_tuple(((sign > 0) ? static_cast<T>(boost::math::erf(guess, Policy()) - target) : static_cast<T>(boost::math::erfc(guess, Policy())) - target), derivative, derivative2);
   }
   erf_roots(T z, int s) : target(z), sign(s) {}
private:
   T target;
   int sign;
};

template <class T, class Policy>
T erf_inv_imp(const T& p, const T& q, const Policy& pol, const boost::mpl::int_<0>*)
{
   //
   // Generic version, get a guess that's accurate to 64-bits (10^-19)
   //
   T guess = erf_inv_imp(p, q, pol, static_cast<mpl::int_<64> const*>(0));
   T result;
   //
   // If T has more bit's than 64 in it's mantissa then we need to iterate,
   // otherwise we can just return the result:
   //
   if(policies::digits<T, Policy>() > 64)
   {
      boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
      if(p <= 0.5)
      {
         result = tools::halley_iterate(detail::erf_roots<typename remove_cv<T>::type, Policy>(p, 1), guess, static_cast<T>(0), tools::max_value<T>(), (policies::digits<T, Policy>() * 2) / 3, max_iter);
      }
      else
      {
         result = tools::halley_iterate(detail::erf_roots<typename remove_cv<T>::type, Policy>(q, -1), guess, static_cast<T>(0), tools::max_value<T>(), (policies::digits<T, Policy>() * 2) / 3, max_iter);
      }
      policies::check_root_iterations<T>("boost::math::erf_inv<%1%>", max_iter, pol);
   }
   else
   {
      result = guess;
   }
   return result;
}

template <class T, class Policy>
struct erf_inv_initializer
{
   struct init
   {
      init()
      {
         do_init();
      }
      static bool is_value_non_zero(T);
      static void do_init()
      {
         // If std::numeric_limits<T>::digits is zero, we must not call
         // our inituialization code here as the precision presumably
         // varies at runtime, and will not have been set yet.
         if(std::numeric_limits<T>::digits)
         {
            boost::math::erf_inv(static_cast<T>(0.25), Policy());
            boost::math::erf_inv(static_cast<T>(0.55), Policy());
            boost::math::erf_inv(static_cast<T>(0.95), Policy());
            boost::math::erfc_inv(static_cast<T>(1e-15), Policy());
            // These following initializations must not be called if
            // type T can not hold the relevant values without
            // underflow to zero.  We check this at runtime because
            // some tools such as valgrind silently change the precision
            // of T at runtime, and numeric_limits basically lies!
            if(is_value_non_zero(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-130))))
               boost::math::erfc_inv(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-130)), Policy());

            // Some compilers choke on constants that would underflow, even in code that isn't instantiated
            // so try and filter these cases out in the preprocessor:
#if LDBL_MAX_10_EXP >= 800
            if(is_value_non_zero(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-800))))
               boost::math::erfc_inv(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-800)), Policy());
            if(is_value_non_zero(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-900))))
               boost::math::erfc_inv(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-900)), Policy());
#else
            if(is_value_non_zero(static_cast<T>(BOOST_MATH_HUGE_CONSTANT(T, 64, 1e-800))))
               boost::math::erfc_inv(static_cast<T>(BOOST_MATH_HUGE_CONSTANT(T, 64, 1e-800)), Policy());
            if(is_value_non_zero(static_cast<T>(BOOST_MATH_HUGE_CONSTANT(T, 64, 1e-900))))
               boost::math::erfc_inv(static_cast<T>(BOOST_MATH_HUGE_CONSTANT(T, 64, 1e-900)), Policy());
#endif
         }
      }
      void force_instantiate()const{}
   };
   static const init initializer;
   static void force_instantiate()
   {
      initializer.force_instantiate();
   }
};

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

template <class T, class Policy>
bool erf_inv_initializer<T, Policy>::init::is_value_non_zero(T v)
{
   // This needs to be non-inline to detect whether v is non zero at runtime
   // rather than at compile time, only relevant when running under valgrind
   // which changes long double's to double's on the fly.
   return v != 0;
}

} // namespace detail

template <class T, class Policy>
typename tools::promote_args<T>::type erfc_inv(T z, const Policy& pol)
{
   typedef typename tools::promote_args<T>::type result_type;

   //
   // Begin by testing for domain errors, and other special cases:
   //
   static const char* function = "boost::math::erfc_inv<%1%>(%1%, %1%)";
   if((z < 0) || (z > 2))
      return policies::raise_domain_error<result_type>(function, "Argument outside range [0,2] in inverse erfc function (got p=%1%).", z, pol);
   if(z == 0)
      return policies::raise_overflow_error<result_type>(function, 0, pol);
   if(z == 2)
      return -policies::raise_overflow_error<result_type>(function, 0, pol);
   //
   // Normalise the input, so it's in the range [0,1], we will
   // negate the result if z is outside that range.  This is a simple
   // application of the erfc reflection formula: erfc(-z) = 2 - erfc(z)
   //
   result_type p, q, s;
   if(z > 1)
   {
      q = 2 - z;
      p = 1 - q;
      s = -1;
   }
   else
   {
      p = 1 - z;
      q = z;
      s = 1;
   }
   //
   // A bit of meta-programming to figure out which implementation
   // to use, based on the number of bits in the mantissa of T:
   //
   typedef typename policies::precision<result_type, Policy>::type precision_type;
   typedef typename mpl::if_<
      mpl::or_<mpl::less_equal<precision_type, mpl::int_<0> >, mpl::greater<precision_type, mpl::int_<64> > >,
      mpl::int_<0>,
      mpl::int_<64>
   >::type tag_type;
   //
   // Likewise use internal promotion, so we evaluate at a higher
   // precision internally if it's appropriate:
   //
   typedef typename policies::evaluation<result_type, Policy>::type eval_type;
   typedef typename policies::normalise<
      Policy,
      policies::promote_float<false>,
      policies::promote_double<false>,
      policies::discrete_quantile<>,
      policies::assert_undefined<> >::type forwarding_policy;

   detail::erf_inv_initializer<eval_type, forwarding_policy>::force_instantiate();

   //
   // And get the result, negating where required:
   //
   return s * policies::checked_narrowing_cast<result_type, forwarding_policy>(
      detail::erf_inv_imp(static_cast<eval_type>(p), static_cast<eval_type>(q), forwarding_policy(), static_cast<tag_type const*>(0)), function);
}

template <class T, class Policy>
typename tools::promote_args<T>::type erf_inv(T z, const Policy& pol)
{
   typedef typename tools::promote_args<T>::type result_type;

   //
   // Begin by testing for domain errors, and other special cases:
   //
   static const char* function = "boost::math::erf_inv<%1%>(%1%, %1%)";
   if((z < -1) || (z > 1))
      return policies::raise_domain_error<result_type>(function, "Argument outside range [-1, 1] in inverse erf function (got p=%1%).", z, pol);
   if(z == 1)
      return policies::raise_overflow_error<result_type>(function, 0, pol);
   if(z == -1)
      return -policies::raise_overflow_error<result_type>(function, 0, pol);
   if(z == 0)
      return 0;
   //
   // Normalise the input, so it's in the range [0,1], we will
   // negate the result if z is outside that range.  This is a simple
   // application of the erf reflection formula: erf(-z) = -erf(z)
   //
   result_type p, q, s;
   if(z < 0)
   {
      p = -z;
      q = 1 - p;
      s = -1;
   }
   else
   {
      p = z;
      q = 1 - z;
      s = 1;
   }
   //
   // A bit of meta-programming to figure out which implementation
   // to use, based on the number of bits in the mantissa of T:
   //
   typedef typename policies::precision<result_type, Policy>::type precision_type;
   typedef typename mpl::if_<
      mpl::or_<mpl::less_equal<precision_type, mpl::int_<0> >, mpl::greater<precision_type, mpl::int_<64> > >,
      mpl::int_<0>,
      mpl::int_<64>
   >::type tag_type;
   //
   // Likewise use internal promotion, so we evaluate at a higher
   // precision internally if it's appropriate:
   //
   typedef typename policies::evaluation<result_type, Policy>::type eval_type;
   typedef typename policies::normalise<
      Policy,
      policies::promote_float<false>,
      policies::promote_double<false>,
      policies::discrete_quantile<>,
      policies::assert_undefined<> >::type forwarding_policy;
   //
   // Likewise use internal promotion, so we evaluate at a higher
   // precision internally if it's appropriate:
   //
   typedef typename policies::evaluation<result_type, Policy>::type eval_type;

   detail::erf_inv_initializer<eval_type, forwarding_policy>::force_instantiate();
   //
   // And get the result, negating where required:
   //
   return s * policies::checked_narrowing_cast<result_type, forwarding_policy>(
      detail::erf_inv_imp(static_cast<eval_type>(p), static_cast<eval_type>(q), forwarding_policy(), static_cast<tag_type const*>(0)), function);
}

template <class T>
inline typename tools::promote_args<T>::type erfc_inv(T z)
{
   return erfc_inv(z, policies::policy<>());
}

template <class T>
inline typename tools::promote_args<T>::type erf_inv(T z)
{
   return erf_inv(z, policies::policy<>());
}

} // namespace math
} // namespace boost

#ifdef _MSC_VER
#pragma warning(pop)
#endif

#endif // BOOST_MATH_SF_ERF_INV_HPP