Imported Upstream version 1.64.0upstream/1.64.0

Change-Id: Id9212edd016dd55f21172c427aa7894d1d24148b
Signed-off-by: DongHun Kwak <dh0128.kwak@samsung.com>

5 files changed, 1998 insertions, 381 deletions

diff --git a/boost/math/special_functions/detail/bessel_i0.hpp b/boost/math/special_functions/detail/bessel_i0.hpp
index 676eb71511..5896db8f76 100644
--- a/boost/math/special_functions/detail/bessel_i0.hpp
+++ b/boost/math/special_functions/detail/bessel_i0.hpp
@@ -1,4 +1,5 @@
 //  Copyright (c) 2006 Xiaogang Zhang
+//  Copyright (c) 2017 John Maddock
 //  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)
@@ -15,28 +16,45 @@
 #include <boost/assert.hpp>
 
 // Modified Bessel function of the first kind of order zero
-// minimax rational approximations on intervals, see
-// Blair and Edwards, Chalk River Report AECL-4928, 1974
+// we use the approximating forms derived in:
+// "Rational Approximations for the Modified Bessel Function of the First Kind – I0(x) for Computations with Double Precision"
+// by Pavel Holoborodko, 
+// see http://www.advanpix.com/2015/11/11/rational-approximations-for-the-modified-bessel-function-of-the-first-kind-i0-computations-double-precision
+// The actual coefficients used are our own, and extend Pavel's work to precision's other than double.
 
 namespace boost { namespace math { namespace detail{
 
 template <typename T>
-T bessel_i0(T x);
+T bessel_i0(const T& x);
 
-template <class T>
+template <class T, class tag>
 struct bessel_i0_initializer
 {
    struct init
    {
       init()
       {
-         do_init();
+         do_init(tag());
       }
-      static void do_init()
+      static void do_init(const mpl::int_<64>&)
       {
          bessel_i0(T(1));
+         bessel_i0(T(8));
+         bessel_i0(T(12));
+         bessel_i0(T(40));
+         bessel_i0(T(101));
       }
-      void force_instantiate()const{}
+      static void do_init(const mpl::int_<113>&)
+      {
+         bessel_i0(T(1));
+         bessel_i0(T(10));
+         bessel_i0(T(20));
+         bessel_i0(T(40));
+         bessel_i0(T(101));
+      }
+      template <class U>
+      static void do_init(const U&) {}
+      void force_instantiate()const {}
    };
    static const init initializer;
    static void force_instantiate()
@@ -45,82 +63,489 @@ struct bessel_i0_initializer
    }
 };
 
-template <class T>
-const typename bessel_i0_initializer<T>::init bessel_i0_initializer<T>::initializer;
+template <class T, class tag>
+const typename bessel_i0_initializer<T, tag>::init bessel_i0_initializer<T, tag>::initializer;
+
+template <typename T, int N>
+T bessel_i0_imp(const T&, const mpl::int_<N>&)
+{
+   BOOST_ASSERT(0);
+   return 0;
+}
+
+template <typename T>
+T bessel_i0_imp(const T& x, const mpl::int_<24>&)
+{
+   BOOST_MATH_STD_USING
+   if(x < 7.75)
+   {
+      // Max error in interpolated form: 3.929e-08
+      // Max Error found at float precision = Poly: 1.991226e-07
+      static const float P[] = {
+         1.00000003928615375e+00f,
+         2.49999576572179639e-01f,
+         2.77785268558399407e-02f,
+         1.73560257755821695e-03f,
+         6.96166518788906424e-05f,
+         1.89645733877137904e-06f,
+         4.29455004657565361e-08f,
+         3.90565476357034480e-10f,
+         1.48095934745267240e-11f
+      };
+      T a = x * x / 4;
+      return a * boost::math::tools::evaluate_polynomial(P, a) + 1;
+   }
+   else if(x < 50)
+   {
+      // Max error in interpolated form: 5.195e-08
+      // Max Error found at float precision = Poly: 8.502534e-08
+      static const float P[] = {
+         3.98942651588301770e-01f,
+         4.98327234176892844e-02f,
+         2.91866904423115499e-02f,
+         1.35614940793742178e-02f,
+         1.31409251787866793e-01f
+      };
+      return exp(x) * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x);
+   }
+   else
+   {
+      // Max error in interpolated form: 1.782e-09
+      // Max Error found at float precision = Poly: 6.473568e-08
+      static const float P[] = {
+         3.98942391532752700e-01f,
+         4.98455950638200020e-02f,
+         2.94835666900682535e-02f
+      };
+      T ex = exp(x / 2);
+      T result = ex * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x);
+      result *= ex;
+      return result;
+   }
+}
+
+template <typename T>
+T bessel_i0_imp(const T& x, const mpl::int_<53>&)
+{
+   BOOST_MATH_STD_USING
+   if(x < 7.75)
+   {
+      // Bessel I0 over[10 ^ -16, 7.75]
+      // Max error in interpolated form : 3.042e-18
+      // Max Error found at double precision = Poly : 5.106609e-16 Cheb : 5.239199e-16
+      static const double P[] = {
+         1.00000000000000000e+00,
+         2.49999999999999909e-01,
+         2.77777777777782257e-02,
+         1.73611111111023792e-03,
+         6.94444444453352521e-05,
+         1.92901234513219920e-06,
+         3.93675991102510739e-08,
+         6.15118672704439289e-10,
+         7.59407002058973446e-12,
+         7.59389793369836367e-14,
+         6.27767773636292611e-16,
+         4.34709704153272287e-18,
+         2.63417742690109154e-20,
+         1.13943037744822825e-22,
+         9.07926920085624812e-25
+      };
+      T a = x * x / 4;
+      return a * boost::math::tools::evaluate_polynomial(P, a) + 1;
+   }
+   else if(x < 500)
+   {
+      // Max error in interpolated form : 1.685e-16
+      // Max Error found at double precision = Poly : 2.575063e-16 Cheb : 2.247615e+00
+      static const double P[] = {
+         3.98942280401425088e-01,
+         4.98677850604961985e-02,
+         2.80506233928312623e-02,
+         2.92211225166047873e-02,
+         4.44207299493659561e-02,
+         1.30970574605856719e-01,
+         -3.35052280231727022e+00,
+         2.33025711583514727e+02,
+         -1.13366350697172355e+04,
+         4.24057674317867331e+05,
+         -1.23157028595698731e+07,
+         2.80231938155267516e+08,
+         -5.01883999713777929e+09,
+         7.08029243015109113e+10,
+         -7.84261082124811106e+11,
+         6.76825737854096565e+12,
+         -4.49034849696138065e+13,
+         2.24155239966958995e+14,
+         -8.13426467865659318e+14,
+         2.02391097391687777e+15,
+         -3.08675715295370878e+15,
+         2.17587543863819074e+15
+      };
+      return exp(x) * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x);
+   }
+   else
+   {
+      // Max error in interpolated form : 2.437e-18
+      // Max Error found at double precision = Poly : 1.216719e-16
+      static const double P[] = {
+         3.98942280401432905e-01,
+         4.98677850491434560e-02,
+         2.80506308916506102e-02,
+         2.92179096853915176e-02,
+         4.53371208762579442e-02
+      };
+      T ex = exp(x / 2);
+      T result = ex * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x);
+      result *= ex;
+      return result;
+   }
+}
 
 template <typename T>
-T bessel_i0(T x)
+T bessel_i0_imp(const T& x, const mpl::int_<64>&)
 {
-    bessel_i0_initializer<T>::force_instantiate();
-
-    static const T P1[] = {
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.2335582639474375249e+15)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -5.5050369673018427753e+14)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.2940087627407749166e+13)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -8.4925101247114157499e+11)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.1912746104985237192e+10)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.0313066708737980747e+08)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -5.9545626019847898221e+05)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.4125195876041896775e+03)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -7.0935347449210549190e+00)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.5453977791786851041e-02)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.5172644670688975051e-05)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.0517226450451067446e-08)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.6843448573468483278e-11)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.5982226675653184646e-14)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -5.2487866627945699800e-18)),
-    };
-    static const T Q1[] = {
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.2335582639474375245e+15)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 7.8858692566751002988e+12)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.2207067397808979846e+10)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0377081058062166144e+07)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -4.8527560179962773045e+03)),
-        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.2210262233306573296e-04)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.3067392038106924055e-02)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -4.4700805721174453923e-01)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.5674518371240761397e+00)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.3517945679239481621e+01)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.1611322818701131207e+01)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -9.6090021968656180000e+00)),
-    };
-    static const T Q2[] = {
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -5.5194330231005480228e-04)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.2547697594819615062e-02)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.1151759188741312645e+00)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.3982595353892851542e+01)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -6.0228002066743340583e+01)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.5539563258012929600e+01)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.1446690275135491500e+01)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)),
-    };
-    T value, factor, r;
-
-    BOOST_MATH_STD_USING
-    using namespace boost::math::tools;
-
-    BOOST_ASSERT(x >= 0); // negative x is handled before we get here
-    if (x == 0)
-    {
-        return static_cast<T>(1);
-    }
-    if (x <= 15)                        // x in (0, 15]
-    {
-        T y = x * x;
-        value = evaluate_polynomial(P1, y) / evaluate_polynomial(Q1, y);
-    }
-    else                                // x in (15, \infty)
-    {
-        T y = 1 / x - T(1) / 15;
-        r = evaluate_polynomial(P2, y) / evaluate_polynomial(Q2, y);
-        factor = exp(x) / sqrt(x);
-        value = factor * r;
-    }
-
-    return value;
+   BOOST_MATH_STD_USING
+   if(x < 7.75)
+   {
+      // Bessel I0 over[10 ^ -16, 7.75]
+      // Max error in interpolated form : 3.899e-20
+      // Max Error found at float80 precision = Poly : 1.770840e-19
+      static const T P[] = {
+         BOOST_MATH_BIG_CONSTANT(T, 64, 9.99999999999999999961011629e-01),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 2.50000000000000001321873912e-01),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 2.77777777777777703400424216e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 1.73611111111112764793802701e-03),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 6.94444444444251461247253525e-05),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 1.92901234569262206386118739e-06),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 3.93675988851131457141005209e-08),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 6.15118734688297476454205352e-10),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 7.59405797058091016449222685e-12),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 7.59406599631719800679835140e-14),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 6.27598961062070013516660425e-16),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 4.35920318970387940278362992e-18),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 2.57372492687715452949437981e-20),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 1.33908663475949906992942204e-22),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 5.15976668870980234582896010e-25),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 3.46240478946376069211156548e-27)
+      };
+      T a = x * x / 4;
+      return a * boost::math::tools::evaluate_polynomial(P, a) + 1;
+   }
+   else if(x < 10)
+   {
+      // Maximum Deviation Found:                     6.906e-21
+      // Expected Error Term : -6.903e-21
+      // Maximum Relative Change in Control Points : 1.631e-04
+      // Max Error found at float80 precision = Poly : 7.811948e-21
+      static const T Y = 4.051098823547363281250e-01f;
+      static const T P[] = {
+         BOOST_MATH_BIG_CONSTANT(T, 64, -6.158081780620616479492e-03),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 4.883635969834048766148e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 7.892782002476195771920e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 64, -1.478784996478070170327e+00),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 2.988611837308006851257e+01),
+         BOOST_MATH_BIG_CONSTANT(T, 64, -4.140133766747436806179e+02),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 4.117316447921276453271e+03),
+         BOOST_MATH_BIG_CONSTANT(T, 64, -2.942353667455141676001e+04),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 1.493482682461387081534e+05),
+         BOOST_MATH_BIG_CONSTANT(T, 64, -5.228100538921466124653e+05),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 1.195279248600467989454e+06),
+         BOOST_MATH_BIG_CONSTANT(T, 64, -1.601530760654337045917e+06),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 9.504921137873298402679e+05)
+      };
+      return exp(x) * (boost::math::tools::evaluate_polynomial(P, T(1 / x)) + Y) / sqrt(x);
+   }
+   else if(x < 15)
+   {
+      // Maximum Deviation Found:                     4.083e-21
+      // Expected Error Term : -4.025e-21
+      // Maximum Relative Change in Control Points : 1.304e-03
+      // Max Error found at float80 precision = Poly : 2.303527e-20
+      static const T Y = 4.033188819885253906250e-01f;
+      static const T P[] = {
+         BOOST_MATH_BIG_CONSTANT(T, 64, -4.376373876116109401062e-03),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 4.982899138682911273321e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 3.109477529533515397644e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 64, -1.163760580110576407673e-01),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 4.776501832837367371883e+00),
+         BOOST_MATH_BIG_CONSTANT(T, 64, -1.101478069227776656318e+02),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 1.892071912448960299773e+03),
+         BOOST_MATH_BIG_CONSTANT(T, 64, -2.417739279982328117483e+04),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 2.296963447724067390552e+05),
+         BOOST_MATH_BIG_CONSTANT(T, 64, -1.598589306710589358747e+06),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 7.903662411851774878322e+06),
+         BOOST_MATH_BIG_CONSTANT(T, 64, -2.622677059040339516093e+07),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 5.227776578828667629347e+07),
+         BOOST_MATH_BIG_CONSTANT(T, 64, -4.727797957441040896878e+07)
+      };
+      return exp(x) * (boost::math::tools::evaluate_polynomial(P, T(1 / x)) + Y) / sqrt(x);
+   }
+   else if(x < 50)
+   {
+      // Max error in interpolated form: 1.035e-21
+      // Max Error found at float80 precision = Poly: 1.885872e-21
+      static const T Y = 4.011702537536621093750e-01f;
+      static const T P[] = {
+         BOOST_MATH_BIG_CONSTANT(T, 64, -2.227973351806078464328e-03),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 4.986778486088017419036e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 2.805066823812285310011e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 2.921443721160964964623e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 4.517504941996594744052e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 6.316922639868793684401e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 1.535891099168810015433e+00),
+         BOOST_MATH_BIG_CONSTANT(T, 64, -4.706078229522448308087e+01),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 1.351015763079160914632e+03),
+         BOOST_MATH_BIG_CONSTANT(T, 64, -2.948809013999277355098e+04),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 4.967598958582595361757e+05),
+         BOOST_MATH_BIG_CONSTANT(T, 64, -6.346924657995383019558e+06),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 5.998794574259956613472e+07),
+         BOOST_MATH_BIG_CONSTANT(T, 64, -4.016371355801690142095e+08),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 1.768791455631826490838e+09),
+         BOOST_MATH_BIG_CONSTANT(T, 64, -4.441995678177349895640e+09),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 4.482292669974971387738e+09)
+      };
+      return exp(x) * (boost::math::tools::evaluate_polynomial(P, T(1 / x)) + Y) / sqrt(x);
+   }
+   else
+   {
+      // Bessel I0 over[50, INF]
+      // Max error in interpolated form : 5.587e-20
+      // Max Error found at float80 precision = Poly : 8.776852e-20
+      static const T P[] = {
+         BOOST_MATH_BIG_CONSTANT(T, 64, 3.98942280401432677955074061e-01),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 4.98677850501789875615574058e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 2.80506290908675604202206833e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 2.92194052159035901631494784e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 4.47422430732256364094681137e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 9.05971614435738691235525172e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 2.29180522595459823234266708e-01),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 6.15122547776140254569073131e-01),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 7.48491812136365376477357324e+00),
+         BOOST_MATH_BIG_CONSTANT(T, 64, -2.45569740166506688169730713e+02),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 9.66857566379480730407063170e+03),
+         BOOST_MATH_BIG_CONSTANT(T, 64, -2.71924083955641197750323901e+05),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 5.74276685704579268845870586e+06),
+         BOOST_MATH_BIG_CONSTANT(T, 64, -8.89753803265734681907148778e+07),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 9.82590905134996782086242180e+08),
+         BOOST_MATH_BIG_CONSTANT(T, 64, -7.30623197145529889358596301e+09),
+         BOOST_MATH_BIG_CONSTANT(T, 64, 3.27310000726207055200805893e+10),
+         BOOST_MATH_BIG_CONSTANT(T, 64, -6.64365417189215599168817064e+10)
+      };
+      T ex = exp(x / 2);
+      T result = ex * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x);
+      result *= ex;
+      return result;
+   }
+}
+
+template <typename T>
+T bessel_i0_imp(const T& x, const mpl::int_<113>&)
+{
+   BOOST_MATH_STD_USING
+   if(x < 7.75)
+   {
+      // Bessel I0 over[10 ^ -34, 7.75]
+      // Max error in interpolated form : 1.274e-34
+      // Max Error found at float128 precision = Poly : 3.096091e-34
+      static const T P[] = {
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000000000001273856e+00),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 2.4999999999999999999999999999999107477496e-01),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 2.7777777777777777777777777777881795230918e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.7361111111111111111111111106290091648808e-03),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 6.9444444444444444444444445629960334523101e-05),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.9290123456790123456790105563456483249753e-06),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 3.9367598891408415217940836339080514004844e-08),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 6.1511873267825648777900014857992724731476e-10),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 7.5940584281266233066162999610732449709209e-12),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 7.5940584281266232783124723601470051895304e-14),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 6.2760813455591936763439337059117957836078e-16),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 4.3583898233049738471136482147779094353096e-18),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 2.5789288895299965395422423848480340736308e-20),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.3157800456718804437960453545507623434606e-22),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 5.8479113149412360748032684260932041506493e-25),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 2.2843403488398038539283241944594140493394e-27),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 7.9042925594356556196790242908697582021825e-30),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 2.4395919891312152120710245152115597111101e-32),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 6.7580986145276689333214547502373003196707e-35),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.6886514018062348877723837017198859723889e-37),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 3.8540558465757554512570197585002702777999e-40),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 7.4684706070226893763741850944911705726436e-43),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 2.0210715309399646335858150349406935414314e-45)
+      };
+      T a = x * x / 4;
+      return a * boost::math::tools::evaluate_polynomial(P, a) + 1;
+   }
+   else if(x < 15)
+   {
+      // Bessel I0 over[7.75, 15]
+      // Max error in interpolated form : 7.534e-35
+      // Max Error found at float128 precision = Poly : 6.123912e-34
+      static const T P[] = {
+         BOOST_MATH_BIG_CONSTANT(T, 113, 9.9999999999999999992388573069504617493518e-01),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 2.5000000000000000007304739268173096975340e-01),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 2.7777777777777777744261405400543564492074e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.7361111111111111209006987259719750726867e-03),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 6.9444444444444442399703186871329381908321e-05),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.9290123456790126709286741580242189785431e-06),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 3.9367598891408374246503061422528266924389e-08),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 6.1511873267826068395343047827801353170966e-10),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 7.5940584281262673459688011737168286944521e-12),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 7.5940584281291583769928563167645746144508e-14),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 6.2760813455438840231126529638737436950274e-16),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 4.3583898233839583885132809584770578894948e-18),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 2.5789288891798658971960571838369339742994e-20),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.3157800470129311623308216856009970266088e-22),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 5.8479112701534604520063520412207286692581e-25),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 2.2843404822552330714586265081801727491890e-27),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 7.9042888166225242675881424439818162458179e-30),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 2.4396027771820721384198604723320045236973e-32),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 6.7577659910606076328136207973456511895030e-35),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.6896548123724136624716224328803899914646e-37),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 3.8285850162160539150210466453921758781984e-40),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 7.9419071894227736216423562425429524883562e-43),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.4720374049498608905571855665134539425038e-45),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 2.7763533278527958112907118930154738930378e-48),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 3.1213839473168678646697528580511702663617e-51),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.0648035313124146852372607519737686740964e-53),
+         -BOOST_MATH_BIG_CONSTANT(T, 113, 5.1255595184052024349371058585102280860878e-57),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 3.4652470895944157957727948355523715335882e-59)
+      };
+      T a = x * x / 4;
+      return a * boost::math::tools::evaluate_polynomial(P, a) + 1;
+   }
+   else if(x < 30)
+   {
+      // Max error in interpolated form : 1.808e-34
+      // Max Error found at float128 precision = Poly : 2.399403e-34
+      static const T P[] = {
+         BOOST_MATH_BIG_CONSTANT(T, 113, 3.9894228040870793650581242239624530714032e-01),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 4.9867780576714783790784348982178607842250e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 2.8051948347934462928487999569249907599510e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 2.8971143420388958551176254291160976367263e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 7.8197359701715582763961322341827341098897e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -3.3430484862908317377522273217643346601271e+00),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 2.7884507603213662610604413960838990199224e+02),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -1.8304926482356755790062999202373909300514e+04),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 9.8867173178574875515293357145875120137676e+05),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -4.4261178812193528551544261731796888257644e+07),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.6453010340778116475788083817762403540097e+09),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -5.0432401330113978669454035365747869477960e+10),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.2462165331309799059332310595587606836357e+12),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -2.3299800389951335932792950236410844978273e+13),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 2.5748218240248714177527965706790413406639e+14),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.8330014378766930869945511450377736037385e+15),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -1.8494610073827453236940544799030787866218e+17),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 5.7244661371420647691301043350229977856476e+18),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -1.2386378807889388140099109087465781254321e+20),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 2.1104000573102013529518477353943384110982e+21),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -2.9426541092239879262282594572224300191016e+22),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 3.4061439136301913488512592402635688101020e+23),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -3.2836554760521986358980180942859101564671e+24),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 2.6270285589905206294944214795661236766988e+25),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -1.7278631455211972017740134341610659484259e+26),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 9.1971734473772196124736986948034978906801e+26),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -3.8669270707172568763908838463689093500098e+27),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.2368879358870281916900125550129211146626e+28),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -2.8296235063297831758204519071113999839858e+28),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 4.1253861666023020670144616019148954773662e+28),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -2.8809536950051955163648980306847791014734e+28) };
+      return exp(x) * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x);
+   }
+   else if(x < 100)
+   {
+      // Bessel I0 over[30, 100]
+      // Max error in interpolated form : 1.487e-34
+      // Max Error found at float128 precision = Poly : 1.929924e-34
+      static const T P[] = {
+         BOOST_MATH_BIG_CONSTANT(T, 113, 3.9894228040143267793996798658172135362278e-01),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 4.9867785050179084714910130342157246539820e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 2.8050629090725751585266360464766768437048e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 2.9219405302833158254515212437025679637597e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 4.4742214371598631578107310396249912330627e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 9.0602983776478659136184969363625092585520e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 2.2839507231977478205885469900971893734770e-01),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 6.8925739165733823730525449511456529001868e-01),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 2.4238082222874015159424842335385854632223e+00),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 9.6759648427182491050716309699208988458050e+00),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 4.7292246491169360014875196108746167872215e+01),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 3.1001411442786230340015781205680362993575e+01),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 9.8277628835804873490331739499978938078848e+03),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -3.1208326312801432038715638596517882759639e+05),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 9.4813611580683862051838126076298945680803e+06),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -2.1278197693321821164135890132925119054391e+08),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 3.3190303792682886967459489059860595063574e+09),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -2.1580767338646580750893606158043485767644e+10),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -5.0256008808415702780816006134784995506549e+11),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.9044186472918017896554580836514681614475e+13),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -3.2521078890073151875661384381880225635135e+14),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 3.3620352486836976842181057590770636605454e+15),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -2.0375525734060401555856465179734887312420e+16),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 5.6392664899881014534361728644608549445131e+16)
+      };
+      return exp(x) * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x);
+   }
+   else
+   {
+      // Bessel I0 over[100, INF]
+      // Max error in interpolated form : 5.459e-35
+      // Max Error found at float128 precision = Poly : 1.472240e-34
+      static const T P[] = {
+         BOOST_MATH_BIG_CONSTANT(T, 113, 3.9894228040143267793994605993438166526772e-01),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 4.9867785050179084742493257495245185241487e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 2.8050629090725735167652437695397756897920e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 2.9219405302839307466358297347675795965363e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 4.4742214369972689474366968442268908028204e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 9.0602984099194778006610058410222616383078e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 2.2839502241666629677015839125593079416327e-01),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 6.8926354981801627920292655818232972385750e-01),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 2.4231921590621824187100989532173995000655e+00),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 9.7264260959693775207585700654645245723497e+00),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 4.3890136225398811195878046856373030127018e+01),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 2.1999720924619285464910452647408431234369e+02),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.2076909538525038580501368530598517194748e+03),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 7.5684635141332367730007149159063086133399e+03),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 3.5178192543258299267923025833141286569141e+04),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 6.2966297919851965784482163987240461837728e+05) };
+      T ex = exp(x / 2);
+      T result = ex * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x);
+      result *= ex;
+      return result;
+   }
+}
+
+template <typename T>
+T bessel_i0_imp(const T& x, const mpl::int_<0>&)
+{
+   if(boost::math::tools::digits<T>() <= 24)
+      return bessel_i0_imp(x, mpl::int_<24>());
+   else if(boost::math::tools::digits<T>() <= 53)
+      return bessel_i0_imp(x, mpl::int_<53>());
+   else if(boost::math::tools::digits<T>() <= 64)
+      return bessel_i0_imp(x, mpl::int_<64>());
+   else if(boost::math::tools::digits<T>() <= 113)
+      return bessel_i0_imp(x, mpl::int_<113>());
+   BOOST_ASSERT(0);
+   return 0;
+}
+
+template <typename T>
+inline T bessel_i0(const T& x)
+{
+   typedef mpl::int_<
+      std::numeric_limits<T>::digits == 0 ?
+      0 :
+      std::numeric_limits<T>::digits <= 24 ?
+      24 :
+      std::numeric_limits<T>::digits <= 53 ?
+      53 :
+      std::numeric_limits<T>::digits <= 64 ?
+      64 :
+      std::numeric_limits<T>::digits <= 113 ?
+      113 : -1
+   > tag_type;
+
+   bessel_i0_initializer<T, tag_type>::force_instantiate();
+   return bessel_i0_imp(x, tag_type());
 }
 
 }}} // namespaces
diff --git a/boost/math/special_functions/detail/bessel_i1.hpp b/boost/math/special_functions/detail/bessel_i1.hpp
index b85bc67546..0f9ba96bd8 100644
--- a/boost/math/special_functions/detail/bessel_i1.hpp
+++ b/boost/math/special_functions/detail/bessel_i1.hpp
@@ -1,7 +1,9 @@
-//  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)
+// Modified Bessel function of the first kind of order zero
+// we use the approximating forms derived in:
+// "Rational Approximations for the Modified Bessel Function of the First Kind – I1(x) for Computations with Double Precision"
+// by Pavel Holoborodko, 
+// see http://www.advanpix.com/2015/11/12/rational-approximations-for-the-modified-bessel-function-of-the-first-kind-i1-for-computations-with-double-precision/
+// The actual coefficients used are our own, and extend Pavel's work to precision's other than double.
 
 #ifndef BOOST_MATH_BESSEL_I1_HPP
 #define BOOST_MATH_BESSEL_I1_HPP
@@ -21,21 +23,36 @@
 namespace boost { namespace math { namespace detail{
 
 template <typename T>
-T bessel_i1(T x);
+T bessel_i1(const T& x);
 
-template <class T>
+template <class T, class tag>
 struct bessel_i1_initializer
 {
    struct init
    {
       init()
       {
-         do_init();
+         do_init(tag());
       }
-      static void do_init()
+      static void do_init(const mpl::int_<64>&)
       {
          bessel_i1(T(1));
+         bessel_i1(T(15));
+         bessel_i1(T(80));
+         bessel_i1(T(101));
       }
+      static void do_init(const mpl::int_<113>&)
+      {
+         bessel_i1(T(1));
+         bessel_i1(T(10));
+         bessel_i1(T(14));
+         bessel_i1(T(19));
+         bessel_i1(T(34));
+         bessel_i1(T(99));
+         bessel_i1(T(101));
+      }
+      template <class U>
+      static void do_init(const U&) {}
       void force_instantiate()const{}
    };
    static const init initializer;
@@ -45,86 +62,513 @@ struct bessel_i1_initializer
    }
 };
 
-template <class T>
-const typename bessel_i1_initializer<T>::init bessel_i1_initializer<T>::initializer;
+template <class T, class tag>
+const typename bessel_i1_initializer<T, tag>::init bessel_i1_initializer<T, tag>::initializer;
+
+template <typename T, int N>
+T bessel_i1_imp(const T&, const mpl::int_<N>&)
+{
+   BOOST_ASSERT(0);
+   return 0;
+}
+
+template <typename T>
+T bessel_i1_imp(const T& x, const mpl::int_<24>&)
+{
+   BOOST_MATH_STD_USING
+      if(x < 7.75)
+      {
+         //Max error in interpolated form : 1.348e-08
+         // Max Error found at float precision = Poly : 1.469121e-07
+         static const float P[] = {
+            8.333333221e-02f,
+            6.944453712e-03f,
+            3.472097211e-04f,
+            1.158047174e-05f,
+            2.739745142e-07f,
+            5.135884609e-09f,
+            5.262251502e-11f,
+            1.331933703e-12f
+         };
+         T a = x * x / 4;
+         T Q[3] = { 1, 0.5f, boost::math::tools::evaluate_polynomial(P, a) };
+         return x * boost::math::tools::evaluate_polynomial(Q, a) / 2;
+      }
+      else
+      {
+         // Max error in interpolated form: 9.000e-08
+         // Max Error found at float precision = Poly: 1.044345e-07
+
+         static const float P[] = {
+            3.98942115977513013e-01f,
+            -1.49581264836620262e-01f,
+            -4.76475741878486795e-02f,
+            -2.65157315524784407e-02f,
+            -1.47148600683672014e-01f
+         };
+         T ex = exp(x / 2);
+         T result = ex * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x);
+         result *= ex;
+         return result;
+      }
+}
+
+template <typename T>
+T bessel_i1_imp(const T& x, const mpl::int_<53>&)
+{
+   BOOST_MATH_STD_USING
+   if(x < 7.75)
+   {
+      // Bessel I0 over[10 ^ -16, 7.75]
+      // Max error in interpolated form: 5.639e-17
+      // Max Error found at double precision = Poly: 1.795559e-16
+
+      static const double P[] = {
+         8.333333333333333803e-02,
+         6.944444444444341983e-03,
+         3.472222222225921045e-04,
+         1.157407407354987232e-05,
+         2.755731926254790268e-07,
+         4.920949692800671435e-09,
+         6.834657311305621830e-11,
+         7.593969849687574339e-13,
+         6.904822652741917551e-15,
+         5.220157095351373194e-17,
+         3.410720494727771276e-19,
+         1.625212890947171108e-21,
+         1.332898928162290861e-23
+      };
+      T a = x * x / 4;
+      T Q[3] = { 1, 0.5f, boost::math::tools::evaluate_polynomial(P, a) };
+      return x * boost::math::tools::evaluate_polynomial(Q, a) / 2;
+   }
+   else if(x < 500)
+   {
+      // Max error in interpolated form: 1.796e-16
+      // Max Error found at double precision = Poly: 2.898731e-16
+
+      static const double P[] = {
+         3.989422804014406054e-01,
+         -1.496033551613111533e-01,
+         -4.675104253598537322e-02,
+         -4.090895951581637791e-02,
+         -5.719036414430205390e-02,
+         -1.528189554374492735e-01,
+         3.458284470977172076e+00,
+         -2.426181371595021021e+02,
+         1.178785865993440669e+04,
+         -4.404655582443487334e+05,
+         1.277677779341446497e+07,
+         -2.903390398236656519e+08,
+         5.192386898222206474e+09,
+         -7.313784438967834057e+10,
+         8.087824484994859552e+11,
+         -6.967602516005787001e+12,
+         4.614040809616582764e+13,
+         -2.298849639457172489e+14,
+         8.325554073334618015e+14,
+         -2.067285045778906105e+15,
+         3.146401654361325073e+15,
+         -2.213318202179221945e+15
+      };
+      return exp(x) * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x);
+   }
+   else
+   {
+      // Max error in interpolated form: 1.320e-19
+      // Max Error found at double precision = Poly: 7.065357e-17
+      static const double P[] = {
+         3.989422804014314820e-01,
+         -1.496033551467584157e-01,
+         -4.675105322571775911e-02,
+         -4.090421597376992892e-02,
+         -5.843630344778927582e-02
+      };
+      T ex = exp(x / 2);
+      T result = ex * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x);
+      result *= ex;
+      return result;
+   }
+}
+
+template <typename T>
+T bessel_i1_imp(const T& x, const mpl::int_<64>&)
+{
+   BOOST_MATH_STD_USING
+      if(x < 7.75)
+      {
+         // Bessel I0 over[10 ^ -16, 7.75]
+         // Max error in interpolated form: 8.086e-21
+         // Max Error found at float80 precision = Poly: 7.225090e-20
+         static const T P[] = {
+            BOOST_MATH_BIG_CONSTANT(T, 64, 8.33333333333333333340071817e-02),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 6.94444444444444442462728070e-03),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 3.47222222222222318886683883e-04),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 1.15740740740738880709555060e-05),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 2.75573192240046222242685145e-07),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 4.92094986131253986838697503e-09),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 6.83465258979924922633502182e-11),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 7.59405830675154933645967137e-13),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 6.90369179710633344508897178e-15),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 5.23003610041709452814262671e-17),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 3.35291901027762552549170038e-19),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 1.83991379419781823063672109e-21),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 8.87732714140192556332037815e-24),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 3.32120654663773147206454247e-26),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 1.95294659305369207813486871e-28) 
+         };
+         T a = x * x / 4;
+         T Q[3] = { 1, 0.5f, boost::math::tools::evaluate_polynomial(P, a) };
+         return x * boost::math::tools::evaluate_polynomial(Q, a) / 2;
+      }
+      else if(x < 20)
+      {
+         // Max error in interpolated form: 4.258e-20
+         // Max Error found at float80 precision = Poly: 2.851105e-19
+         // Maximum Deviation Found : 3.887e-20
+         // Expected Error Term : 3.887e-20
+         // Maximum Relative Change in Control Points : 1.681e-04
+         static const T P[] = {
+            BOOST_MATH_BIG_CONSTANT(T, 64, 3.98942260530218897338680e-01),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -1.49599542849073670179540e-01),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -4.70492865454119188276875e-02),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -3.12389893307392002405869e-02),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 1.49696126385202602071197e-01),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -3.84206507612717711565967e+01),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 2.14748094784412558689584e+03),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -7.70652726663596993005669e+04),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 2.01659736164815617174439e+06),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -4.04740659606466305607544e+07),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 6.38383394696382837263656e+08),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -8.00779638649147623107378e+09),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 8.02338237858684714480491e+10),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -6.41198553664947312995879e+11),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 4.05915186909564986897554e+12),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -2.00907636964168581116181e+13),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 7.60855263982359981275199e+13),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -2.12901817219239205393806e+14),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 4.14861794397709807823575e+14),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -5.02808138522587680348583e+14),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 2.85505477056514919387171e+14)
+         };
+         return exp(x) * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x);
+      }
+      else if(x < 100)
+      {
+         // Bessel I0 over [15, 50]
+         // Maximum Deviation Found:                     2.444e-20
+         // Expected Error Term : 2.438e-20
+         // Maximum Relative Change in Control Points : 2.101e-03
+         // Max Error found at float80 precision = Poly : 6.029974e-20
+
+         static const T P[] = {
+            BOOST_MATH_BIG_CONSTANT(T, 64, 3.98942280401431675205845e-01),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -1.49603355149968887210170e-01),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -4.67510486284376330257260e-02),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -4.09071458907089270559464e-02),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -5.75278280327696940044714e-02),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -1.10591299500956620739254e-01),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -2.77061766699949309115618e-01),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -5.42683771801837596371638e-01),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -9.17021412070404158464316e+00),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 1.04154379346763380543310e+02),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -1.43462345357478348323006e+03),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 9.98109660274422449523837e+03),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -3.74438822767781410362757e+04)
+         };
+         return exp(x) * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x);
+      }
+      else
+      {
+         // Bessel I0 over[100, INF]
+         // Max error in interpolated form: 2.456e-20
+         // Max Error found at float80 precision = Poly: 5.446356e-20
+         static const T P[] = {
+            BOOST_MATH_BIG_CONSTANT(T, 64, 3.98942280401432677958445e-01),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -1.49603355150537411254359e-01),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -4.67510484842456251368526e-02),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -4.09071676503922479645155e-02),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -5.75256179814881566010606e-02),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -1.10754910257965227825040e-01),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -2.67858639515616079840294e-01),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -9.17266479586791298924367e-01)
+         };
+         T ex = exp(x / 2);
+         T result = ex * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x);
+         result *= ex;
+         return result;
+      }
+}
+
+template <typename T>
+T bessel_i1_imp(const T& x, const mpl::int_<113>&)
+{
+   BOOST_MATH_STD_USING
+   if(x < 7.75)
+   {
+      // Bessel I0 over[10 ^ -34, 7.75]
+      // Max error in interpolated form: 1.835e-35
+      // Max Error found at float128 precision = Poly: 1.645036e-34
+
+      static const T P[] = {
+         BOOST_MATH_BIG_CONSTANT(T, 113, 8.3333333333333333333333333333333331804098e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 6.9444444444444444444444444444445418303082e-03),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 3.4722222222222222222222222222119082346591e-04),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.1574074074074074074074074078415867655987e-05),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 2.7557319223985890652557318255143448192453e-07),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 4.9209498614260519022423916850415000626427e-09),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 6.8346525853139609753354247043900442393686e-11),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 7.5940584281266233060080535940234144302217e-13),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 6.9036894801151120925605467963949641957095e-15),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 5.2300677879659941472662086395055636394839e-17),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 3.3526075563884539394691458717439115962233e-19),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.8420920639497841692288943167036233338434e-21),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 8.7718669711748690065381181691546032291365e-24),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 3.6549445715236427401845636880769861424730e-26),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.3437296196812697924703896979250126739676e-28),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 4.3912734588619073883015937023564978854893e-31),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.2839967682792395867255384448052781306897e-33),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 3.3790094235693528861015312806394354114982e-36),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 8.0423861671932104308662362292359563970482e-39),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.7493858979396446292135661268130281652945e-41),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 3.2786079392547776769387921361408303035537e-44),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 8.2335693685833531118863552173880047183822e-47)
+      };
+      T a = x * x / 4;
+      T Q[3] = { 1, 0.5f, boost::math::tools::evaluate_polynomial(P, a) };
+      return x * boost::math::tools::evaluate_polynomial(Q, a) / 2;
+   }
+   else if(x < 11)
+   {
+      // Max error in interpolated form: 8.574e-36
+      // Maximum Deviation Found : 4.689e-36
+      // Expected Error Term : 3.760e-36
+      // Maximum Relative Change in Control Points : 5.204e-03
+      // Max Error found at float128 precision = Poly : 2.882561e-34
+
+      static const T P[] = {
+         BOOST_MATH_BIG_CONSTANT(T, 113, 8.333333333333333326889717360850080939e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 6.944444444444444511272790848815114507e-03),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 3.472222222222221892451965054394153443e-04),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.157407407407408437378868534321538798e-05),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 2.755731922398566216824909767320161880e-07),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 4.920949861426434829568192525456800388e-09),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 6.834652585308926245465686943255486934e-11),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 7.594058428179852047689599244015979196e-13),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 6.903689479655006062822949671528763738e-15),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 5.230067791254403974475987777406992984e-17),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 3.352607536815161679702105115200693346e-19),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.842092161364672561828681848278567885e-21),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 8.771862912600611801856514076709932773e-24),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 3.654958704184380914803366733193713605e-26),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.343688672071130980471207297730607625e-28),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 4.392252844664709532905868749753463950e-31),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.282086786672692641959912811902298600e-33),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 3.408812012322547015191398229942864809e-36),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 7.681220437734066258673404589233009892e-39),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 2.072417451640733785626701738789290055e-41),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.352218520142636864158849446833681038e-44),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.407918492276267527897751358794783640e-46)
+      };
+      T a = x * x / 4;
+      T Q[3] = { 1, 0.5f, boost::math::tools::evaluate_polynomial(P, a) };
+      return x * boost::math::tools::evaluate_polynomial(Q, a) / 2;
+   }
+   else if(x < 15)
+   {
+      //Max error in interpolated form: 7.599e-36
+      // Maximum Deviation Found : 1.766e-35
+      // Expected Error Term : 1.021e-35
+      // Maximum Relative Change in Control Points : 6.228e-03
+      static const T P[] = {
+         BOOST_MATH_BIG_CONSTANT(T, 113, 8.333333333333255774414858563409941233e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 6.944444444444897867884955912228700291e-03),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 3.472222222220954970397343617150959467e-04),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.157407407409660682751155024932538578e-05),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 2.755731922369973706427272809014190998e-07),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 4.920949861702265600960449699129258153e-09),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 6.834652583208361401197752793379677147e-11),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 7.594058441128280500819776168239988143e-13),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 6.903689413939268702265479276217647209e-15),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 5.230068069012898202890718644753625569e-17),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 3.352606552027491657204243201021677257e-19),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.842095100698532984651921750204843362e-21),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 8.771789051329870174925649852681844169e-24),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 3.655114381199979536997025497438385062e-26),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.343415732516712339472538688374589373e-28),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 4.396177019032432392793591204647901390e-31),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.277563309255167951005939802771456315e-33),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 3.449201419305514579791370198046544736e-36),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 7.415430703400740634202379012388035255e-39),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 2.195458831864936225409005027914934499e-41),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 8.829726762743879793396637797534668039e-45),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.698302711685624490806751012380215488e-46),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -2.062520475425422618494185821587228317e-49),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 6.732372906742845717148185173723304360e-52)
+      };
+      T a = x * x / 4;
+      T Q[3] = { 1, 0.5f, boost::math::tools::evaluate_polynomial(P, a) };
+      return x * boost::math::tools::evaluate_polynomial(Q, a) / 2;
+   }
+   else if(x < 20)
+   {
+      // Max error in interpolated form: 8.864e-36
+      // Max Error found at float128 precision = Poly: 8.522841e-35
+      static const T P[] = {
+         BOOST_MATH_BIG_CONSTANT(T, 113, 3.989422793693152031514179994954750043e-01),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -1.496029423752889591425633234009799670e-01),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -4.682975926820553021482820043377990241e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -3.138871171577224532369979905856458929e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -8.765350219426341341990447005798111212e-01),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 5.321389275507714530941178258122955540e+01),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -2.727748393898888756515271847678850411e+03),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.123040820686242586086564998713862335e+05),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -3.784112378374753535335272752884808068e+06),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.054920416060932189433079126269416563e+08),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -2.450129415468060676827180524327749553e+09),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 4.758831882046487398739784498047935515e+10),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -7.736936520262204842199620784338052937e+11),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.051128683324042629513978256179115439e+13),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -1.188008285959794869092624343537262342e+14),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.108530004906954627420484180793165669e+15),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -8.441516828490144766650287123765318484e+15),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 5.158251664797753450664499268756393535e+16),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -2.467314522709016832128790443932896401e+17),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 8.896222045367960462945885220710294075e+17),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -2.273382139594876997203657902425653079e+18),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 3.669871448568623680543943144842394531e+18),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -2.813923031370708069940575240509912588e+18)
+      };
+      return exp(x) * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x);
+   }
+   else if(x < 35)
+   {
+      // Max error in interpolated form: 6.028e-35
+      // Max Error found at float128 precision = Poly: 1.368313e-34
+
+      static const T P[] = {
+         BOOST_MATH_BIG_CONSTANT(T, 113, 3.989422804012941975429616956496046931e-01),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -1.496033550576049830976679315420681402e-01),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -4.675107835141866009896710750800622147e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -4.090104965125365961928716504473692957e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -5.842241652296980863361375208605487570e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -1.063604828033747303936724279018650633e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -9.113375972811586130949401996332817152e+00),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 6.334748570425075872639817839399823709e+02),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -3.759150758768733692594821032784124765e+04),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.863672813448915255286274382558526321e+06),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -7.798248643371718775489178767529282534e+07),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 2.769963173932801026451013022000669267e+09),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -8.381780137198278741566746511015220011e+10),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 2.163891337116820832871382141011952931e+12),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -4.764325864671438675151635117936912390e+13),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 8.925668307403332887856809510525154955e+14),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -1.416692606589060039334938090985713641e+16),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.892398600219306424294729851605944429e+17),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -2.107232903741874160308537145391245060e+18),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.930223393531877588898224144054112045e+19),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -1.427759576167665663373350433236061007e+20),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 8.306019279465532835530812122374386654e+20),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -3.653753000392125229440044977239174472e+21),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.140760686989511568435076842569804906e+22),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -2.249149337812510200795436107962504749e+22),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 2.101619088427348382058085685849420866e+22)
+      };
+      return exp(x) * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x);
+   }
+   else if(x < 100)
+   {
+      // Max error in interpolated form: 5.494e-35
+      // Max Error found at float128 precision = Poly: 1.214651e-34
+
+      static const T P[] = {
+         BOOST_MATH_BIG_CONSTANT(T, 113, 3.989422804014326779399307367861631577e-01),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -1.496033551505372542086590873271571919e-01),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -4.675104848454290286276466276677172664e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -4.090716742397105403027549796269213215e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -5.752570419098513588311026680089351230e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -1.107369803696534592906420980901195808e-01),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -2.699214194000085622941721628134575121e-01),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -7.953006169077813678478720427604462133e-01),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -2.746618809476524091493444128605380593e+00),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -1.084446249943196826652788161656973391e+01),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -5.020325182518980633783194648285500554e+01),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -1.510195971266257573425196228564489134e+02),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -5.241661863814900938075696173192225056e+03),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.323374362891993686413568398575539777e+05),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -4.112838452096066633754042734723911040e+06),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 9.369270194978310081563767560113534023e+07),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -1.704295412488936504389347368131134993e+09),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 2.320829576277038198439987439508754886e+10),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -2.258818139077875493434420764260185306e+11),
+         BOOST_MATH_BIG_CONSTANT(T, 113, 1.396791306321498426110315039064592443e+12),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -4.217617301585849875301440316301068439e+12)
+      };
+      return exp(x) * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x);
+   }
+   else
+   {
+      // Bessel I0 over[100, INF]
+      // Max error in interpolated form: 6.081e-35
+      // Max Error found at float128 precision = Poly: 1.407151e-34
+      static const T P[] = {
+         BOOST_MATH_BIG_CONSTANT(T, 113, 3.9894228040143267793994605993438200208417e-01),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -1.4960335515053725422747977247811372936584e-01),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -4.6751048484542891946087411826356811991039e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -4.0907167423975030452875828826630006305665e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -5.7525704189964886494791082898669060345483e-02),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -1.1073698056568248642163476807108190176386e-01),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -2.6992139012879749064623499618582631684228e-01),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -7.9530409594026597988098934027440110587905e-01),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -2.7462844478733532517044536719240098183686e+00),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -1.0870711340681926669381449306654104739256e+01),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -4.8510175413216969245241059608553222505228e+01),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -2.4094682286011573747064907919522894740063e+02),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -1.3128845936764406865199641778959502795443e+03),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -8.1655901321962541203257516341266838487359e+03),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -3.8019591025686295090160445920753823994556e+04),
+         BOOST_MATH_BIG_CONSTANT(T, 113, -6.7008089049178178697338128837158732831105e+05)
+      };
+      T ex = exp(x / 2);
+      T result = ex * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x);
+      result *= ex;
+      return result;
+   }
+}
+
+template <typename T>
+T bessel_i1_imp(const T& x, const mpl::int_<0>&)
+{
+   if(boost::math::tools::digits<T>() <= 24)
+      return bessel_i1_imp(x, mpl::int_<24>());
+   else if(boost::math::tools::digits<T>() <= 53)
+      return bessel_i1_imp(x, mpl::int_<53>());
+   else if(boost::math::tools::digits<T>() <= 64)
+      return bessel_i1_imp(x, mpl::int_<64>());
+   else if(boost::math::tools::digits<T>() <= 113)
+      return bessel_i1_imp(x, mpl::int_<113>());
+   BOOST_ASSERT(0);
+   return 0;
+}
 
 template <typename T>
-T bessel_i1(T x)
+inline T bessel_i1(const T& x)
 {
+   typedef mpl::int_<
+      std::numeric_limits<T>::digits == 0 ?
+      0 :
+      std::numeric_limits<T>::digits <= 24 ?
+      24 :
+      std::numeric_limits<T>::digits <= 53 ?
+      53 :
+      std::numeric_limits<T>::digits <= 64 ?
+      64 :
+      std::numeric_limits<T>::digits <= 113 ?
+      113 : -1
+   > tag_type;
 
-    bessel_i1_initializer<T>::force_instantiate();
-
-    static const T P1[] = {
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.4577180278143463643e+15)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.7732037840791591320e+14)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -6.9876779648010090070e+12)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.3357437682275493024e+11)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.4828267606612366099e+09)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.0588550724769347106e+07)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -5.1894091982308017540e+04)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.8225946631657315931e+02)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -4.7207090827310162436e-01)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -9.1746443287817501309e-04)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.3466829827635152875e-06)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.4831904935994647675e-09)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.1928788903603238754e-12)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -6.5245515583151902910e-16)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.9705291802535139930e-19)),
-    };
-    static const T Q1[] = {
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.9154360556286927285e+15)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 9.7887501377547640438e+12)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.4386907088588283434e+10)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.1594225856856884006e+07)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -5.1326864679904189920e+03)),
-        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.4582087408985668208e-05)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -8.9359825138577646443e-04)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.9204895411257790122e-02)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.4198728018058047439e-01)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.3960118277609544334e+00)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.9746376087200685843e+00)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.5591872901933459000e-01)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -6.0437159056137599999e-02)),
-    };
-    static const T Q2[] = {
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.7510433111922824643e-05)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.2835624489492512649e-03)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 7.4212010813186530069e-02)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -8.5017476463217924408e-01)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.2593714889036996297e+00)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.8806586721556593450e+00)),
-        static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)),
-    };
-    T value, factor, r, w;
-
-    BOOST_MATH_STD_USING
-    using namespace boost::math::tools;
-
-    BOOST_ASSERT(x >= 0); // negative x is handled before we get here
-    w = abs(x);
-    if (x == 0)
-    {
-        return static_cast<T>(0);
-    }
-    if (w <= 15)                        // w in (0, 15]
-    {
-        T y = x * x;
-        r = evaluate_polynomial(P1, y) / evaluate_polynomial(Q1, y);
-        factor = w;
-        value = factor * r;
-    }
-    else                                // w in (15, \infty)
-    {
-        T y = 1 / w - T(1) / 15;
-        r = evaluate_polynomial(P2, y) / evaluate_polynomial(Q2, y);
-        factor = exp(w) / sqrt(w);
-        value = factor * r;
-    }
-
-    return value;
+   bessel_i1_initializer<T, tag_type>::force_instantiate();
+   return bessel_i1_imp(x, tag_type());
 }
 
 }}} // namespaces
diff --git a/boost/math/special_functions/detail/bessel_k0.hpp b/boost/math/special_functions/detail/bessel_k0.hpp
index 42041dc673..4f3420ffc0 100644
--- a/boost/math/special_functions/detail/bessel_k0.hpp
+++ b/boost/math/special_functions/detail/bessel_k0.hpp
@@ -1,4 +1,5 @@
 //  Copyright (c) 2006 Xiaogang Zhang
+//  Copyright (c) 2017 John Maddock
 //  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)
@@ -19,26 +20,43 @@
 
 // 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
+// Russon and Blair, Chalk River Report AECL-3461, 1969,
+// as revised by Pavel Holoborodko in "Rational Approximations 
+// for the Modified Bessel Function of the Second Kind - K0(x) 
+// for Computations with Double Precision", see 
+// http://www.advanpix.com/2015/11/25/rational-approximations-for-the-modified-bessel-function-of-the-second-kind-k0-for-computations-with-double-precision/
+//
+// The actual coefficients used are our own derivation (by JM)
+// since we extend to both greater and lesser precision than the
+// references above.  We can also improve performance WRT to
+// Holoborodko without loss of precision.
 
 namespace boost { namespace math { namespace detail{
 
-template <typename T, typename Policy>
-T bessel_k0(T x, const Policy&);
+template <typename T>
+T bessel_k0(const T& x);
 
-template <class T, class Policy>
+template <class T, class tag>
 struct bessel_k0_initializer
 {
    struct init
    {
       init()
       {
-         do_init();
+         do_init(tag());
       }
-      static void do_init()
+      static void do_init(const mpl::int_<113>&)
       {
-         bessel_k0(T(1), Policy());
+         bessel_k0(T(0.5));
+         bessel_k0(T(1.5));
       }
+      static void do_init(const mpl::int_<64>&)
+      {
+         bessel_k0(T(0.5));
+         bessel_k0(T(1.5));
+      }
+      template <class U>
+      static void do_init(const U&){}
       void force_instantiate()const{}
    };
    static const init initializer;
@@ -48,104 +66,437 @@ struct bessel_k0_initializer
    }
 };
 
-template <class T, class Policy>
-const typename bessel_k0_initializer<T, Policy>::init bessel_k0_initializer<T, Policy>::initializer;
+template <class T, class tag>
+const typename bessel_k0_initializer<T, tag>::init bessel_k0_initializer<T, tag>::initializer;
+
+
+template <typename T, int N>
+T bessel_k0_imp(const T& x, const mpl::int_<N>&)
+{
+   BOOST_ASSERT(0);
+   return 0;
+}
+
+template <typename T>
+T bessel_k0_imp(const T& x, const mpl::int_<24>&)
+{
+   BOOST_MATH_STD_USING
+   if(x <= 1)
+   {
+      // Maximum Deviation Found : 2.358e-09
+      // Expected Error Term : -2.358e-09
+      // Maximum Relative Change in Control Points : 9.552e-02
+      // Max Error found at float precision = Poly : 4.448220e-08
+      static const T Y = 1.137250900268554688f;
+      static const T P[] = 
+      {
+         -1.372508979104259711e-01f,
+         2.622545986273687617e-01f,
+         5.047103728247919836e-03f
+      };
+      static const T Q[] = 
+      {
+         1.000000000000000000e+00f,
+         -8.928694018000029415e-02f,
+         2.985980684180969241e-03f
+      };
+      T a = x * x / 4;
+      a = (tools::evaluate_rational(P, Q, a) + Y) * a + 1;
+
+      // Maximum Deviation Found:                     1.346e-09
+      // Expected Error Term : -1.343e-09
+      // Maximum Relative Change in Control Points : 2.405e-02
+      // Max Error found at float precision = Poly : 1.354814e-07
+      static const T P2[] = {
+         1.159315158e-01f,
+         2.789828686e-01f,
+         2.524902861e-02f,
+         8.457241514e-04f,
+         1.530051997e-05f
+      };
+      return tools::evaluate_polynomial(P2, T(x * x)) - log(x) * a;
+   }
+   else
+   {
+      // Maximum Deviation Found:                     1.587e-08
+      // Expected Error Term : 1.531e-08
+      // Maximum Relative Change in Control Points : 9.064e-02
+      // Max Error found at float precision = Poly : 5.065020e-08
+
+      static const T P[] =
+      {
+         2.533141220e-01,
+         5.221502603e-01,
+         6.380180669e-02,
+         -5.934976547e-02
+      };
+      static const T Q[] =
+      {
+         1.000000000e+00,
+         2.679722431e+00,
+         1.561635813e+00,
+         1.573660661e-01
+      };
+      if(x < tools::log_max_value<T>())
+         return ((tools::evaluate_rational(P, Q, T(1 / x)) + 1) * exp(-x) / sqrt(x));
+      else
+      {
+         T ex = exp(-x / 2);
+         return ((tools::evaluate_rational(P, Q, T(1 / x)) + 1) * ex / sqrt(x)) * ex;
+      }
+   }
+}
+
+template <typename T>
+T bessel_k0_imp(const T& x, const mpl::int_<53>&)
+{
+   BOOST_MATH_STD_USING
+   if(x <= 1)
+   {
+      // Maximum Deviation Found:                     6.077e-17
+      // Expected Error Term : -6.077e-17
+      // Maximum Relative Change in Control Points : 7.797e-02
+      // Max Error found at double precision = Poly : 1.003156e-16
+      static const T Y = 1.137250900268554688;
+      static const T P[] =
+      {
+         -1.372509002685546267e-01,
+         2.574916117833312855e-01,
+         1.395474602146869316e-02,
+         5.445476986653926759e-04,
+         7.125159422136622118e-06
+      };
+      static const T Q[] =
+      {
+         1.000000000000000000e+00,
+         -5.458333438017788530e-02,
+         1.291052816975251298e-03,
+         -1.367653946978586591e-05
+      };
+
+      T a = x * x / 4;
+      a = (tools::evaluate_polynomial(P, a) / tools::evaluate_polynomial(Q, a) + Y) * a + 1;
+
+      // Maximum Deviation Found:                     3.429e-18
+      // Expected Error Term : 3.392e-18
+      // Maximum Relative Change in Control Points : 2.041e-02
+      // Max Error found at double precision = Poly : 2.513112e-16
+      static const T P2[] =
+      {
+         1.159315156584124484e-01,
+         2.789828789146031732e-01,
+         2.524892993216121934e-02,
+         8.460350907213637784e-04,
+         1.491471924309617534e-05,
+         1.627106892422088488e-07,
+         1.208266102392756055e-09,
+         6.611686391749704310e-12
+      };
+
+      return tools::evaluate_polynomial(P2, T(x * x)) - log(x) * a;
+   }
+   else
+   {
+      // Maximum Deviation Found:                     4.316e-17
+      // Expected Error Term : 9.570e-18
+      // Maximum Relative Change in Control Points : 2.757e-01
+      // Max Error found at double precision = Poly : 1.001560e-16
+
+      static const T Y = 1;
+      static const T P[] =
+      {
+         2.533141373155002416e-01,
+         3.628342133984595192e+00,
+         1.868441889406606057e+01,
+         4.306243981063412784e+01,
+         4.424116209627428189e+01,
+         1.562095339356220468e+01,
+         -1.810138978229410898e+00,
+         -1.414237994269995877e+00,
+         -9.369168119754924625e-02
+      };
+      static const T Q[] =
+      {
+         1.000000000000000000e+00,
+         1.494194694879908328e+01,
+         8.265296455388554217e+01,
+         2.162779506621866970e+02,
+         2.845145155184222157e+02,
+         1.851714491916334995e+02,
+         5.486540717439723515e+01,
+         6.118075837628957015e+00,
+         1.586261269326235053e-01
+      };
+      if(x < tools::log_max_value<T>())
+         return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * exp(-x) / sqrt(x));
+      else
+      {
+         T ex = exp(-x / 2);
+         return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * ex / sqrt(x)) * ex;
+      }
+   }
+}
 
-template <typename T, typename Policy>
-T bessel_k0(T x, const Policy& pol)
+template <typename T>
+T bessel_k0_imp(const T& x, const mpl::int_<64>&)
 {
-    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;
+   BOOST_MATH_STD_USING
+      if(x <= 1)
+      {
+         // Maximum Deviation Found:                     2.180e-22
+         // Expected Error Term : 2.180e-22
+         // Maximum Relative Change in Control Points : 2.943e-01
+         // Max Error found at float80 precision = Poly : 3.923207e-20
+         static const T Y = 1.137250900268554687500e+00;
+         static const T P[] =
+         {
+            BOOST_MATH_BIG_CONSTANT(T, 64, -1.372509002685546875002e-01),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 2.566481981037407600436e-01),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 1.551881122448948854873e-02),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 6.646112454323276529650e-04),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 1.213747930378196492543e-05),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 9.423709328020389560844e-08)
+         };
+         static const T Q[] =
+         {
+            BOOST_MATH_BIG_CONSTANT(T, 64, 1.000000000000000000000e+00),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -4.843828412587773008342e-02),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 1.088484822515098936140e-03),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -1.374724008530702784829e-05),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 8.452665455952581680339e-08)
+         };
+
+
+         T a = x * x / 4;
+         a = (tools::evaluate_polynomial(P, a) / tools::evaluate_polynomial(Q, a) + Y) * a + 1;
+
+         // Maximum Deviation Found:                     2.440e-21
+         // Expected Error Term : -2.434e-21
+         // Maximum Relative Change in Control Points : 2.459e-02
+         // Max Error found at float80 precision = Poly : 1.482487e-19
+         static const T P2[] =
+         {
+            BOOST_MATH_BIG_CONSTANT(T, 64, 1.159315156584124488110e-01),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 2.764832791416047889734e-01),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 1.926062887220923354112e-02),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 3.660777862036966089410e-04),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 2.094942446930673386849e-06)
+         };
+         static const T Q2[] =
+         {
+            BOOST_MATH_BIG_CONSTANT(T, 64, 1.000000000000000000000e+00),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -2.156100313881251616320e-02),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 2.315993873344905957033e-04),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -1.529444499350703363451e-06),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 5.524988589917857531177e-09)
+         };
+         return tools::evaluate_rational(P2, Q2, T(x * x)) - log(x) * a;
+      }
+      else
+      {
+         // Maximum Deviation Found:                     4.291e-20
+         // Expected Error Term : 2.236e-21
+         // Maximum Relative Change in Control Points : 3.021e-01
+         //Max Error found at float80 precision = Poly : 8.727378e-20
+         static const T Y = 1;
+         static const T P[] =
+         {
+            BOOST_MATH_BIG_CONSTANT(T, 64, 2.533141373155002512056e-01),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 5.417942070721928652715e+00),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 4.477464607463971754433e+01),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 1.838745728725943889876e+02),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 4.009736314927811202517e+02),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 4.557411293123609803452e+02),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 2.360222564015361268955e+02),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 2.385435333168505701022e+01),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -1.750195760942181592050e+01),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -4.059789241612946683713e+00),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -1.612783121537333908889e-01)
+         };
+         static const T Q[] =
+         {
+            BOOST_MATH_BIG_CONSTANT(T, 64, 1.000000000000000000000e+00),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 2.200669254769325861404e+01),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 1.900177593527144126549e+02),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 8.361003989965786932682e+02),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 2.041319870804843395893e+03),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 2.828491555113790345068e+03),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 2.190342229261529076624e+03),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 9.003330795963812219852e+02),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 1.773371397243777891569e+02),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 1.368634935531158398439e+01),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 2.543310879400359967327e-01)
+         };
+         if(x < tools::log_max_value<T>())
+            return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * exp(-x) / sqrt(x));
+         else
+         {
+            T ex = exp(-x / 2);
+            return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * ex / sqrt(x)) * ex;
+         }
+      }
+}
+
+template <typename T>
+T bessel_k0_imp(const T& x, const mpl::int_<113>&)
+{
+   BOOST_MATH_STD_USING
+      if(x <= 1)
+      {
+         // Maximum Deviation Found:                     5.682e-37
+         // Expected Error Term : 5.682e-37
+         // Maximum Relative Change in Control Points : 6.094e-04
+         // Max Error found at float128 precision = Poly : 5.338213e-35
+         static const T Y = 1.137250900268554687500000000000000000e+00f;
+         static const T P[] =
+         {
+            BOOST_MATH_BIG_CONSTANT(T, 113, -1.372509002685546875000000000000000006e-01),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 2.556212905071072782462974351698081303e-01),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.742459135264203478530904179889103929e-02),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 8.077860530453688571555479526961318918e-04),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.868173911669241091399374307788635148e-05),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 2.496405768838992243478709145123306602e-07),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.752489221949580551692915881999762125e-09),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 5.243010555737173524710512824955368526e-12)
+         };
+         static const T Q[] =
+         {
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.000000000000000000000000000000000000e+00),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -4.095631064064621099785696980653193721e-02),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 8.313880983725212151967078809725835532e-04),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -1.095229912293480063501285562382835142e-05),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.022828799511943141130509410251996277e-07),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -6.860874007419812445494782795829046836e-10),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 3.107297802344970725756092082686799037e-12),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -7.460529579244623559164763757787600944e-15)
+         };
+         T a = x * x / 4;
+         a = (tools::evaluate_rational(P, Q, a) + Y) * a + 1;
+
+         // Maximum Deviation Found:                     5.173e-38
+         // Expected Error Term : 5.105e-38
+         // Maximum Relative Change in Control Points : 9.734e-03
+         // Max Error found at float128 precision = Poly : 1.688806e-34
+         static const T P2[] =
+         {
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.159315156584124488107200313757741370e-01),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 2.789828789146031122026800078439435369e-01),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 2.524892993216269451266750049024628432e-02),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 8.460350907082229957222453839935101823e-04),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.491471929926042875260452849503857976e-05),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.627105610481598430816014719558896866e-07),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.208426165007797264194914898538250281e-09),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 6.508697838747354949164182457073784117e-12),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 2.659784680639805301101014383907273109e-14),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 8.531090131964391104248859415958109654e-17),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 2.205195117066478034260323124669936314e-19),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 4.692219280289030165761119775783115426e-22),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 8.362350161092532344171965861545860747e-25),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.277990623924628999539014980773738258e-27)
+         };
+
+         return tools::evaluate_polynomial(P2, T(x * x)) - log(x) * a;
+      }
+      else
+      {
+         // Maximum Deviation Found:                     1.462e-34
+         // Expected Error Term : 4.917e-40
+         // Maximum Relative Change in Control Points : 3.385e-01
+         // Max Error found at float128 precision = Poly : 1.567573e-34
+         static const T Y = 1;
+         static const T P[] =
+         {
+            BOOST_MATH_BIG_CONSTANT(T, 113, 2.533141373155002512078826424055226265e-01),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 2.001949740768235770078339977110749204e+01),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 6.991516715983883248363351472378349986e+02),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.429587951594593159075690819360687720e+04),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.911933815201948768044660065771258450e+05),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.769943016204926614862175317962439875e+06),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.170866154649560750500954150401105606e+07),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 5.634687099724383996792011977705727661e+07),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.989524036456492581597607246664394014e+08),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 5.160394785715328062088529400178080360e+08),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 9.778173054417826368076483100902201433e+08),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.335667778588806892764139643950439733e+09),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.283635100080306980206494425043706838e+09),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 8.300616188213640626577036321085025855e+08),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 3.277591957076162984986406540894621482e+08),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 5.564360536834214058158565361486115932e+07),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -1.043505161612403359098596828115690596e+07),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -7.217035248223503605127967970903027314e+06),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -1.422938158797326748375799596769964430e+06),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -1.229125746200586805278634786674745210e+05),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -4.201632288615609937883545928660649813e+03),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -3.690820607338480548346746717311811406e+01)
+         };
+         static const T Q[] =
+         {
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.000000000000000000000000000000000000e+00),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 7.964877874035741452203497983642653107e+01),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 2.808929943826193766839360018583294769e+03),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 5.814524004679994110944366890912384139e+04),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 7.897794522506725610540209610337355118e+05),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 7.456339470955813675629523617440433672e+06),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 5.057818717813969772198911392875127212e+07),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 2.513821619536852436424913886081133209e+08),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 9.255938846873380596038513316919990776e+08),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 2.537077551699028079347581816919572141e+09),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 5.176769339768120752974843214652367321e+09),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 7.828722317390455845253191337207432060e+09),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 8.698864296569996402006511705803675890e+09),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 7.007803261356636409943826918468544629e+09),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 4.016564631288740308993071395104715469e+09),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.595893010619754750655947035567624730e+09),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 4.241241839120481076862742189989406856e+08),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 7.168778094393076220871007550235840858e+07),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 7.156200301360388147635052029404211109e+06),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 3.752130382550379886741949463587008794e+05),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 8.370574966987293592457152146806662562e+03),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 4.871254714311063594080644835895740323e+01)
+         };
+         if(x < tools::log_max_value<T>())
+            return  ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * exp(-x) / sqrt(x));
+         else
+         {
+            T ex = exp(-x / 2);
+            return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * ex / sqrt(x)) * ex;
+         }
+      }
+}
+
+template <typename T>
+T bessel_k0_imp(const T& x, const mpl::int_<0>&)
+{
+   if(boost::math::tools::digits<T>() <= 24)
+      return bessel_k0_imp(x, mpl::int_<24>());
+   else if(boost::math::tools::digits<T>() <= 53)
+      return bessel_k0_imp(x, mpl::int_<53>());
+   else if(boost::math::tools::digits<T>() <= 64)
+      return bessel_k0_imp(x, mpl::int_<64>());
+   else if(boost::math::tools::digits<T>() <= 113)
+      return bessel_k0_imp(x, mpl::int_<113>());
+   BOOST_ASSERT(0);
+   return 0;
+}
+
+template <typename T>
+inline T bessel_k0(const T& x)
+{
+   typedef mpl::int_<
+      std::numeric_limits<T>::digits == 0 ?
+      0 :
+      std::numeric_limits<T>::digits <= 24 ?
+      24 :
+      std::numeric_limits<T>::digits <= 53 ?
+      53 :
+      std::numeric_limits<T>::digits <= 64 ?
+      64 :
+      std::numeric_limits<T>::digits <= 113 ?
+      113 : -1
+   > tag_type;
+
+   bessel_k0_initializer<T, tag_type>::force_instantiate();
+   return bessel_k0_imp(x, tag_type());
 }
 
 }}} // namespaces
diff --git a/boost/math/special_functions/detail/bessel_k1.hpp b/boost/math/special_functions/detail/bessel_k1.hpp
index c0c5478520..0ae90056cb 100644
--- a/boost/math/special_functions/detail/bessel_k1.hpp
+++ b/boost/math/special_functions/detail/bessel_k1.hpp
@@ -1,4 +1,5 @@
 //  Copyright (c) 2006 Xiaogang Zhang
+//  Copyright (c) 2017 John Maddock 
 //  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)
@@ -17,132 +18,528 @@
 #include <boost/math/policies/error_handling.hpp>
 #include <boost/assert.hpp>
 
-// Modified Bessel function of the second kind of order one
+// 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
+// Russon and Blair, Chalk River Report AECL-3461, 1969,
+// as revised by Pavel Holoborodko in "Rational Approximations 
+// for the Modified Bessel Function of the Second Kind - K0(x) 
+// for Computations with Double Precision", see 
+// http://www.advanpix.com/2016/01/05/rational-approximations-for-the-modified-bessel-function-of-the-second-kind-k1-for-computations-with-double-precision/
+//
+// The actual coefficients used are our own derivation (by JM)
+// since we extend to both greater and lesser precision than the
+// references above.  We can also improve performance WRT to
+// Holoborodko without loss of precision.
 
 namespace boost { namespace math { namespace detail{
 
-template <typename T, typename Policy>
-T bessel_k1(T x, const Policy&);
+   template <typename T>
+   T bessel_k1(const T& x);
 
-template <class T, class Policy>
-struct bessel_k1_initializer
-{
-   struct init
+   template <class T, class tag>
+   struct bessel_k1_initializer
    {
-      init()
+      struct init
       {
-         do_init();
-      }
-      static void do_init()
+         init()
+         {
+            do_init(tag());
+         }
+         static void do_init(const mpl::int_<113>&)
+         {
+            bessel_k1(T(0.5));
+            bessel_k1(T(2));
+            bessel_k1(T(6));
+         }
+         static void do_init(const mpl::int_<64>&)
+         {
+            bessel_k1(T(0.5));
+            bessel_k1(T(6));
+         }
+         template <class U>
+         static void do_init(const U&) {}
+         void force_instantiate()const {}
+      };
+      static const init initializer;
+      static void force_instantiate()
       {
-         bessel_k1(T(1), Policy());
+         initializer.force_instantiate();
       }
-      void force_instantiate()const{}
    };
-   static const init initializer;
-   static void force_instantiate()
+
+   template <class T, class tag>
+   const typename bessel_k1_initializer<T, tag>::init bessel_k1_initializer<T, tag>::initializer;
+
+
+   template <typename T, int N>
+   inline T bessel_k1_imp(const T& x, const mpl::int_<N>&)
    {
-      initializer.force_instantiate();
+      BOOST_ASSERT(0);
+      return 0;
    }
-};
-
-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;
+
+   template <typename T>
+   T bessel_k1_imp(const T& x, const mpl::int_<24>&)
+   {
+      BOOST_MATH_STD_USING
+      if(x <= 1)
+      {
+         // Maximum Deviation Found:                     3.090e-12
+         // Expected Error Term : -3.053e-12
+         // Maximum Relative Change in Control Points : 4.927e-02
+         // Max Error found at float precision = Poly : 7.918347e-10
+         static const T Y = 8.695471287e-02f;
+         static const T P[] =
+         {
+            -3.621379531e-03f,
+            7.131781976e-03f,
+            -1.535278300e-05f
+         };
+         static const T Q[] =
+         {
+            1.000000000e+00f,
+            -5.173102701e-02f,
+            9.203530671e-04f
+         };
+
+         T a = x * x / 4;
+         a = ((tools::evaluate_rational(P, Q, a) + Y) * a * a + a / 2 + 1) * x / 2;
+
+         // Maximum Deviation Found:                     3.556e-08
+         // Expected Error Term : -3.541e-08
+         // Maximum Relative Change in Control Points : 8.203e-02
+         static const T P2[] =
+         {
+            -3.079657469e-01f,
+            -8.537108913e-02f,
+            -4.640275408e-03f,
+            -1.156442414e-04f
+         };
+
+         return tools::evaluate_polynomial(P2, T(x * x)) * x + 1 / x + log(x) * a;
+      }
+      else
+      {
+         // Maximum Deviation Found:                     3.369e-08
+         // Expected Error Term : -3.227e-08
+         // Maximum Relative Change in Control Points : 9.917e-02
+         // Max Error found at float precision = Poly : 6.084411e-08
+         static const T Y = 1.450342178f;
+         static const T P[] =
+         {
+            -1.970280088e-01f,
+            2.188747807e-02f,
+            7.270394756e-01f,
+            2.490678196e-01f
+         };
+         static const T Q[] =
+         {
+            1.000000000e+00f,
+            2.274292882e+00f,
+            9.904984851e-01f,
+            4.585534549e-02f
+         };
+         if(x < tools::log_max_value<T>())
+            return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * exp(-x) / sqrt(x));
+         else
+         {
+            T ex = exp(-x / 2);
+            return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * ex / sqrt(x)) * ex;
+         }
+      }
+   }
+
+   template <typename T>
+   T bessel_k1_imp(const T& x, const mpl::int_<53>&)
+   {
+      BOOST_MATH_STD_USING
+      if(x <= 1)
+      {
+         // Maximum Deviation Found:                     1.922e-17
+         // Expected Error Term : 1.921e-17
+         // Maximum Relative Change in Control Points : 5.287e-03
+         // Max Error found at double precision = Poly : 2.004747e-17
+         static const T Y = 8.69547128677368164e-02f;
+         static const T P[] =
+         {
+            -3.62137953440350228e-03,
+            7.11842087490330300e-03,
+            1.00302560256614306e-05,
+            1.77231085381040811e-06
+         };
+         static const T Q[] =
+         {
+            1.00000000000000000e+00,
+            -4.80414794429043831e-02,
+            9.85972641934416525e-04,
+            -8.91196859397070326e-06
+         };
+
+         T a = x * x / 4;
+         a = ((tools::evaluate_rational(P, Q, a) + Y) * a * a + a / 2 + 1) * x / 2;
+
+         // Maximum Deviation Found:                     4.053e-17
+         // Expected Error Term : -4.053e-17
+         // Maximum Relative Change in Control Points : 3.103e-04
+         // Max Error found at double precision = Poly : 1.246698e-16
+
+         static const T P2[] =
+         {
+            -3.07965757829206184e-01,
+            -7.80929703673074907e-02,
+            -2.70619343754051620e-03,
+            -2.49549522229072008e-05
+         };
+         static const T Q2[] = 
+         {
+            1.00000000000000000e+00,
+            -2.36316836412163098e-02,
+            2.64524577525962719e-04,
+            -1.49749618004162787e-06
+         };
+
+         return tools::evaluate_rational(P2, Q2, T(x * x)) * x + 1 / x + log(x) * a;
+      }
+      else
+      {
+         // Maximum Deviation Found:                     8.883e-17
+         // Expected Error Term : -1.641e-17
+         // Maximum Relative Change in Control Points : 2.786e-01
+         // Max Error found at double precision = Poly : 1.258798e-16
+
+         static const T Y = 1.45034217834472656f;
+         static const T P[] =
+         {
+            -1.97028041029226295e-01,
+            -2.32408961548087617e+00,
+            -7.98269784507699938e+00,
+            -2.39968410774221632e+00,
+            3.28314043780858713e+01,
+            5.67713761158496058e+01,
+            3.30907788466509823e+01,
+            6.62582288933739787e+00,
+            3.08851840645286691e-01
+         };
+         static const T Q[] =
+         {
+            1.00000000000000000e+00,
+            1.41811409298826118e+01,
+            7.35979466317556420e+01,
+            1.77821793937080859e+02,
+            2.11014501598705982e+02,
+            1.19425262951064454e+02,
+            2.88448064302447607e+01,
+            2.27912927104139732e+00,
+            2.50358186953478678e-02
+         };
+         if(x < tools::log_max_value<T>())
+            return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * exp(-x) / sqrt(x));
+         else
+         {
+            T ex = exp(-x / 2);
+            return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * ex / sqrt(x)) * ex;
+         }
+      }
+   }
+
+   template <typename T>
+   T bessel_k1_imp(const T& x, const mpl::int_<64>&)
+   {
+      BOOST_MATH_STD_USING
+      if(x <= 1)
+      {
+         // Maximum Deviation Found:                     5.549e-23
+         // Expected Error Term : -5.548e-23
+         // Maximum Relative Change in Control Points : 2.002e-03
+         // Max Error found at float80 precision = Poly : 9.352785e-22
+         static const T Y = 8.695471286773681640625e-02f;
+         static const T P[] =
+         {
+            BOOST_MATH_BIG_CONSTANT(T, 64, -3.621379534403483072861e-03),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 7.102135866103952705932e-03),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 4.167545240236717601167e-05),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 2.537484002571894870830e-06),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 6.603228256820000135990e-09)
+         };
+         static const T Q[] =
+         {
+            BOOST_MATH_BIG_CONSTANT(T, 64, 1.000000000000000000000e+00),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -4.354457194045068370363e-02),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 8.709137201220209072820e-04),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -9.676151796359590545143e-06),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 5.162715192766245311659e-08)
+         };
+
+         T a = x * x / 4;
+         a = ((tools::evaluate_rational(P, Q, a) + Y) * a * a + a / 2 + 1) * x / 2;
+
+         // Maximum Deviation Found:                     1.995e-23
+         // Expected Error Term : 1.995e-23
+         // Maximum Relative Change in Control Points : 8.174e-04
+         // Max Error found at float80 precision = Poly : 4.137325e-20
+         static const T P2[] =
+         {
+            BOOST_MATH_BIG_CONSTANT(T, 64, -3.079657578292062244054e-01),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -7.963049154965966503231e-02),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -3.103277523735639924895e-03),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -4.023052834702215699504e-05),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -1.719459155018493821839e-07)
+         };
+         static const T Q2[] = 
+         {
+            BOOST_MATH_BIG_CONSTANT(T, 64, 1.000000000000000000000e+00),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -1.863917670410152669768e-02),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 1.699367098849735298090e-04),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -9.309358790546076298429e-07),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 2.708893480271612711933e-09)
+         };
+
+         return tools::evaluate_rational(P2, Q2, T(x * x)) * x + 1 / x + log(x) * a;
+      }
+      else
+      {
+         // Maximum Deviation Found:                     9.785e-20
+         // Expected Error Term : -3.302e-21
+         // Maximum Relative Change in Control Points : 3.432e-01
+         // Max Error found at float80 precision = Poly : 1.083755e-19
+         static const T Y = 1.450342178344726562500e+00f;
+         static const T P[] =
+         {
+            BOOST_MATH_BIG_CONSTANT(T, 64, -1.970280410292263112917e-01),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -4.058564803062959169322e+00),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -3.036658174194917777473e+01),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -9.576825392332820142173e+01),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -6.706969489248020941949e+01),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 3.264572499406168221382e+02),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 8.584972047303151034100e+02),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 8.422082733280017909550e+02),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 3.738005441471368178383e+02),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 7.016938390144121276609e+01),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 4.319614662598089438939e+00),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 3.710715864316521856193e-02)
+         };
+         static const T Q[] =
+         {
+            BOOST_MATH_BIG_CONSTANT(T, 64, 1.000000000000000000000e+00),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 2.298433045824439052398e+01),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 2.082047745067709230037e+02),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 9.662367854250262046592e+02),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 2.504148628460454004686e+03),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 3.712730364911389908905e+03),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 3.108002081150068641112e+03),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 1.400149940532448553143e+03),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 3.083303048095846226299e+02),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 2.748706060530351833346e+01),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 6.321900849331506946977e-01),
+         };
+         if(x < tools::log_max_value<T>())
+            return ((tools::evaluate_polynomial(P, T(1 / x)) / tools::evaluate_polynomial(Q, T(1 / x)) + Y) * exp(-x) / sqrt(x));
+         else
+         {
+            T ex = exp(-x / 2);
+            return ((tools::evaluate_polynomial(P, T(1 / x)) / tools::evaluate_polynomial(Q, T(1 / x)) + Y) * ex / sqrt(x)) * ex;
+         }
+      }
+   }
+
+   template <typename T>
+   T bessel_k1_imp(const T& x, const mpl::int_<113>&)
+   {
+      BOOST_MATH_STD_USING
+      if(x <= 1)
+      {
+         // Maximum Deviation Found:                     7.120e-35
+         // Expected Error Term : -7.119e-35
+         // Maximum Relative Change in Control Points : 1.207e-03
+         // Max Error found at float128 precision = Poly : 7.143688e-35
+         static const T Y = 8.695471286773681640625000000000000000e-02f;
+         static const T P[] =
+         {
+            BOOST_MATH_BIG_CONSTANT(T, 113, -3.621379534403483072916666666666595475e-03),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 7.074117676930975433219826471336547627e-03),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 9.631337631362776369069668419033041661e-05),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 3.468935967870048731821071646104412775e-06),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 2.956705020559599861444492614737168261e-08),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 2.347140307321161346703214099534250263e-10),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 5.569608494081482873946791086435679661e-13)
+         };
+         static const T Q[] =
+         {
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.000000000000000000000000000000000000e+00),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -3.580768910152105375615558920428350204e-02),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 6.197467671701485365363068445534557369e-04),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -6.707466533308630411966030561446666237e-06),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 4.846687802282250112624373388491123527e-08),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -2.248493131151981569517383040323900343e-10),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 5.319279786372775264555728921709381080e-13)
+         };
+
+         T a = x * x / 4;
+         a = ((tools::evaluate_rational(P, Q, a) + Y) * a * a + a / 2 + 1) * x / 2;
+
+         // Maximum Deviation Found:                     4.473e-37
+         // Expected Error Term : 4.473e-37
+         // Maximum Relative Change in Control Points : 8.550e-04
+         // Max Error found at float128 precision = Poly : 8.167701e-35
+         static const T P2[] =
+         {
+            BOOST_MATH_BIG_CONSTANT(T, 113, -3.079657578292062244053600156878870690e-01),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -8.133183745732467770755578848987414875e-02),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -3.548968792764174773125420229299431951e-03),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -5.886125468718182876076972186152445490e-05),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -4.506712111733707245745396404449639865e-07),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -1.632502325880313239698965376754406011e-09),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -2.311973065898784812266544485665624227e-12)
+         };
+         static const T Q2[] = 
+         {
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.000000000000000000000000000000000000e+00),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -1.311471216733781016657962995723287450e-02),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 8.571876054797365417068164018709472969e-05),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -3.630181215268238731442496851497901293e-07),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.070176111227805048604885986867484807e-09),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -2.129046580769872602793220056461084761e-12),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 2.294906469421390890762001971790074432e-15)
+         };
+
+         return tools::evaluate_rational(P2, Q2, T(x * x)) * x + 1 / x + log(x) * a;
+      }
+      else if(x < 4)
+      {
+         // Max error in interpolated form: 5.307e-37
+         // Max Error found at float128 precision = Poly: 7.087862e-35
+         static const T Y = 1.5023040771484375f;
+         static const T P[] =
+         {
+            BOOST_MATH_BIG_CONSTANT(T, 113, -2.489899398329369710528254347931380044e-01),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -6.819080211203854781858815596508456873e+00),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -7.599915699069767382647695624952723034e+01),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -4.450211910821295507926582231071300718e+02),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -1.451374687870925175794150513723956533e+03),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -2.405805746895098802803503988539098226e+03),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -5.638808326778389656403861103277220518e+02),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 5.513958744081268456191778822780865708e+03),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.121301640926540743072258116122834804e+04),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.080094900175649541266613109971296190e+04),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 5.896531083639613332407534434915552429e+03),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.856602122319645694042555107114028437e+03),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 3.237121918853145421414003823957537419e+02),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 2.842072954561323076230238664623893504e+01),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.039705646510167437971862966128055524e+00),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.008418100718254816100425022904039530e-02)
+         };
+         static const T Q[] =
+         {
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.000000000000000000000000000000000000e+00),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 2.927456835239137986889227412815459529e+01),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 3.598985593265577043711382994516531273e+02),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 2.449897377085510281395819892689690579e+03),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.025555887684561913263090023158085327e+04),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 2.774140447181062463181892531100679195e+04),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 4.962055507843204417243602332246120418e+04),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 5.908269326976180183216954452196772931e+04),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 4.655160454422016855911700790722577942e+04),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 2.383586885019548163464418964577684608e+04),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 7.679920375586960324298491662159976419e+03),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.478586421028842906987799049804565008e+03),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.565384974896746094224942654383537090e+02),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 7.902617937084010911005732488607114511e+00),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.429293010387921526110949911029094926e-01),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 3.880342607911083143560111853491047663e-04)
+         };
+         return ((tools::evaluate_polynomial(P, T(1 / x)) / tools::evaluate_polynomial(Q, T(1 / x)) + Y) * exp(-x) / sqrt(x));
+      }
+      else
+      {
+         // Maximum Deviation Found:                     4.359e-37
+         // Expected Error Term : -6.565e-40
+         // Maximum Relative Change in Control Points : 1.880e-01
+         // Max Error found at float128 precision = Poly : 2.943572e-35
+         static const T Y = 1.308816909790039062500000000000000000f;
+         static const T P[] =
+         {
+            BOOST_MATH_BIG_CONSTANT(T, 113, -5.550277247453881129211735759447737350e-02),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -3.485883080219574328217554864956175929e+00),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -8.903760658131484239300875153154881958e+01),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -1.144813672213626237418235110712293337e+03),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -6.498400501156131446691826557494158173e+03),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.573531831870363502604119835922166116e+04),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 5.417416550054632009958262596048841154e+05),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 4.271266450613557412825896604269130661e+06),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.898386013314389952534433455681107783e+07),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 5.353798784656436259250791761023512750e+07),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 9.839619195427352438957774052763490067e+07),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.169246368651532232388152442538005637e+08),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 8.696368884166831199967845883371116431e+07),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 3.810226630422736458064005843327500169e+07),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 8.854996610560406127438950635716757614e+06),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 8.981057433937398731355768088809437625e+05),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 2.519440069856232098711793483639792952e+04)
+         };
+         static const T Q[] =
+         {
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.000000000000000000000000000000000000e+00),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 7.127348248283623146544565916604103560e+01),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 2.205092684176906740104488180754982065e+03),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 3.911249195069050636298346469740075758e+04),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 4.426103406579046249654548481377792614e+05),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 3.365861555422488771286500241966208541e+06),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.765377714160383676864913709252529840e+07),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 6.453822726931857253365138260720815246e+07),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.643207885048369990391975749439783892e+08),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 2.882540678243694621895816336640877878e+08),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 3.410120808992380266174106812005338148e+08),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 2.628138016559335882019310900426773027e+08),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.250794693811010646965360198541047961e+08),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 3.378723408195485594610593014072950078e+07),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 4.488253856312453816451380319061865560e+06),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 2.202167197882689873967723350537104582e+05),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.673233230356966539460728211412989843e+03)
+         };
+         if(x < tools::log_max_value<T>())
+            return ((tools::evaluate_polynomial(P, T(1 / x)) / tools::evaluate_polynomial(Q, T(1 / x)) + Y) * exp(-x) / sqrt(x));
+         else
+         {
+            T ex = exp(-x / 2);
+            return ((tools::evaluate_polynomial(P, T(1 / x)) / tools::evaluate_polynomial(Q, T(1 / x)) + Y) * ex / sqrt(x)) * ex;
+         }
+      }
     }
-    else                                // x in (1, \infty)
+
+    template <typename T>
+    T bessel_k1_imp(const T& x, const mpl::int_<0>&)
     {
-        T y = 1 / x;
-        r = evaluate_polynomial(P3, y) /  evaluate_polynomial(Q3, y);
-        factor = exp(-x) / sqrt(x);
-        value = factor * r;
+       if(boost::math::tools::digits<T>() <= 24)
+          return bessel_k1_imp(x, mpl::int_<24>());
+       else if(boost::math::tools::digits<T>() <= 53)
+          return bessel_k1_imp(x, mpl::int_<53>());
+       else if(boost::math::tools::digits<T>() <= 64)
+          return bessel_k1_imp(x, mpl::int_<64>());
+       else if(boost::math::tools::digits<T>() <= 113)
+          return bessel_k1_imp(x, mpl::int_<113>());
+       BOOST_ASSERT(0);
+       return 0;
     }
 
-    return value;
-}
+    template <typename T>
+   inline T bessel_k1(const T& x)
+   {
+      typedef mpl::int_<
+         std::numeric_limits<T>::digits == 0 ?
+         0 :
+         std::numeric_limits<T>::digits <= 24 ?
+         24 :
+         std::numeric_limits<T>::digits <= 53 ?
+         53 :
+         std::numeric_limits<T>::digits <= 64 ?
+         64 :
+         std::numeric_limits<T>::digits <= 113 ?
+         113 : -1
+      > tag_type;
+
+      bessel_k1_initializer<T, tag_type>::force_instantiate();
+      return bessel_k1_imp(x, tag_type());
+   }
 
 }}} // namespaces
 
diff --git a/boost/math/special_functions/detail/bessel_kn.hpp b/boost/math/special_functions/detail/bessel_kn.hpp
index e3a5023c63..54c4a1cfa7 100644
--- a/boost/math/special_functions/detail/bessel_kn.hpp
+++ b/boost/math/special_functions/detail/bessel_kn.hpp
@@ -45,16 +45,16 @@ T bessel_kn(int n, T x, const Policy& pol)
     }
     if (n == 0)
     {
-        value = bessel_k0(x, pol);
+        value = bessel_k0(x);
     }
     else if (n == 1)
     {
-        value = bessel_k1(x, pol);
+        value = bessel_k1(x);
     }
     else
     {
-       prev = bessel_k0(x, pol);
-       current = bessel_k1(x, pol);
+       prev = bessel_k0(x);
+       current = bessel_k1(x);
        int k = 1;
        BOOST_ASSERT(k < n);
        T scale = 1;