mirror of
https://github.com/boostorg/math.git
synced 2026-01-19 04:22:09 +00:00
Add improved approximations for K0 and K1.
Based on http://www.advanpix.com/2016/01/05/rational-approximations-for-the-modified-bessel-function-of-the-second-kind-k1-for-computations-with-double-precision/
This commit is contained in:
jzmaddock
@@ -572,7 +572,7 @@ inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_k(
|
||||
{
|
||||
BOOST_FPU_EXCEPTION_GUARD
|
||||
typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
|
||||
typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type;
|
||||
typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag128 tag_type;
|
||||
typedef typename policies::evaluation<result_type, Policy>::type value_type;
|
||||
typedef typename policies::normalise<
|
||||
Policy,
|
||||
|
||||
@@ -24,21 +24,29 @@
|
||||
namespace boost { namespace math { namespace detail{
|
||||
|
||||
template <typename T, typename Policy>
|
||||
T bessel_k0(T x, const Policy&);
|
||||
T bessel_k0(const T& x, const Policy&);
|
||||
|
||||
template <class T, class Policy>
|
||||
template <class T, class Policy, 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), Policy());
|
||||
bessel_k0(T(1.5), Policy());
|
||||
}
|
||||
static void do_init(const mpl::int_<64>&)
|
||||
{
|
||||
bessel_k0(T(0.5), Policy());
|
||||
bessel_k0(T(1.5), Policy());
|
||||
}
|
||||
template <class U>
|
||||
static void do_init(const U&){}
|
||||
void force_instantiate()const{}
|
||||
};
|
||||
static const init initializer;
|
||||
@@ -48,104 +56,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 Policy, class tag>
|
||||
const typename bessel_k0_initializer<T, Policy, tag>::init bessel_k0_initializer<T, Policy, tag>::initializer;
|
||||
|
||||
template <typename T, typename Policy>
|
||||
T bessel_k0(T x, const Policy& pol)
|
||||
|
||||
template <typename T, class Policy, int N>
|
||||
T bessel_k0_imp(const T& x, const Policy&, const mpl::int_<N>&)
|
||||
{
|
||||
BOOST_MATH_INSTRUMENT_CODE(x);
|
||||
BOOST_ASSERT(0);
|
||||
return 0;
|
||||
}
|
||||
|
||||
bessel_k0_initializer<T, Policy>::force_instantiate();
|
||||
template <typename T, class Policy>
|
||||
T bessel_k0_imp(const T& x, const Policy& pol, const mpl::int_<0>&)
|
||||
{
|
||||
if(boost::math::tools::digits<T>() <= 24)
|
||||
return bessel_k0_imp(x, pol, mpl::int_<24>());
|
||||
else if(boost::math::tools::digits<T>() <= 53)
|
||||
return bessel_k0_imp(x, pol, mpl::int_<53>());
|
||||
else if(boost::math::tools::digits<T>() <= 64)
|
||||
return bessel_k0_imp(x, pol, mpl::int_<64>());
|
||||
else if(boost::math::tools::digits<T>() <= 113)
|
||||
return bessel_k0_imp(x, pol, mpl::int_<113>());
|
||||
BOOST_ASSERT(0);
|
||||
return 0;
|
||||
}
|
||||
|
||||
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;
|
||||
template <typename T, class Policy>
|
||||
T bessel_k0_imp(const T& x, const Policy& pol, 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;
|
||||
|
||||
BOOST_MATH_STD_USING
|
||||
using namespace boost::math::tools;
|
||||
// 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, 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 char* function = "boost::math::bessel_k0<%1%>(%1%,%1%)";
|
||||
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;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
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);
|
||||
}
|
||||
template <typename T, class Policy>
|
||||
T bessel_k0_imp(const T& x, const Policy& pol, 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
|
||||
};
|
||||
|
||||
return value;
|
||||
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, 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, class Policy>
|
||||
T bessel_k0_imp(const T& x, const Policy& pol, const mpl::int_<64>&)
|
||||
{
|
||||
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, 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, class Policy>
|
||||
T bessel_k0_imp(const T& x, const Policy& pol, 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, 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, class Policy>
|
||||
inline T bessel_k0(const T& x, const Policy& pol)
|
||||
{
|
||||
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, Policy, tag_type>::force_instantiate();
|
||||
return bessel_k0_imp(x, pol, tag_type());
|
||||
}
|
||||
|
||||
}}} // namespaces
|
||||
|
||||
@@ -23,126 +23,513 @@
|
||||
|
||||
namespace boost { namespace math { namespace detail{
|
||||
|
||||
template <typename T, typename Policy>
|
||||
T bessel_k1(T x, const Policy&);
|
||||
template <typename T, typename Policy>
|
||||
T bessel_k1(const T& x, const Policy&);
|
||||
|
||||
template <class T, class Policy>
|
||||
struct bessel_k1_initializer
|
||||
{
|
||||
struct init
|
||||
template <class T, class Policy, 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), Policy());
|
||||
bessel_k1(T(2), Policy());
|
||||
bessel_k1(T(6), Policy());
|
||||
}
|
||||
static void do_init(const mpl::int_<64>&)
|
||||
{
|
||||
bessel_k1(T(0.5), Policy());
|
||||
bessel_k1(T(6), Policy());
|
||||
}
|
||||
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 Policy, class tag>
|
||||
const typename bessel_k1_initializer<T, Policy, tag>::init bessel_k1_initializer<T, Policy, tag>::initializer;
|
||||
|
||||
|
||||
template <typename T, class Policy, int N>
|
||||
inline T bessel_k1_imp(const T& x, const Policy&, 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, class Policy>
|
||||
T bessel_k1_imp(const T& x, const Policy& pol, const mpl::int_<0>&)
|
||||
{
|
||||
if(boost::math::tools::digits<T>() <= 24)
|
||||
return bessel_k1_imp(x, pol, mpl::int_<24>());
|
||||
else if(boost::math::tools::digits<T>() <= 53)
|
||||
return bessel_k1_imp(x, pol, mpl::int_<53>());
|
||||
else if(boost::math::tools::digits<T>() <= 64)
|
||||
return bessel_k1_imp(x, pol, mpl::int_<64>());
|
||||
else if(boost::math::tools::digits<T>() <= 113)
|
||||
return bessel_k1_imp(x, pol, mpl::int_<113>());
|
||||
BOOST_ASSERT(0);
|
||||
return 0;
|
||||
}
|
||||
|
||||
template <typename T, typename Policy>
|
||||
T bessel_k1(T x, const Policy& pol)
|
||||
{
|
||||
bessel_k1_initializer<T, Policy>::force_instantiate();
|
||||
template <typename T, class Policy>
|
||||
T bessel_k1_imp(const T& x, const Policy& pol, 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
|
||||
};
|
||||
|
||||
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;
|
||||
T a = x * x / 4;
|
||||
a = ((tools::evaluate_rational(P, Q, a) + Y) * a * a + a / 2 + 1) * x / 2;
|
||||
|
||||
BOOST_MATH_STD_USING
|
||||
using namespace boost::math::tools;
|
||||
// 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
|
||||
};
|
||||
|
||||
static const char* function = "boost::math::bessel_k1<%1%>(%1%,%1%)";
|
||||
return tools::evaluate_polynomial(P2, 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;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (x < 0)
|
||||
{
|
||||
return policies::raise_domain_error<T>(function,
|
||||
"Got x = %1%, but argument x must be non-negative, complex number result not supported.", x, pol);
|
||||
}
|
||||
if (x == 0)
|
||||
{
|
||||
return policies::raise_overflow_error<T>(function, 0, pol);
|
||||
}
|
||||
if (x <= 1) // x in (0, 1]
|
||||
{
|
||||
T y = x * x;
|
||||
r1 = evaluate_polynomial(P1, y) / evaluate_polynomial(Q1, y);
|
||||
r2 = evaluate_polynomial(P2, y) / evaluate_polynomial(Q2, y);
|
||||
factor = log(x);
|
||||
value = (r1 + factor * r2) / x;
|
||||
}
|
||||
else // x in (1, \infty)
|
||||
{
|
||||
T y = 1 / x;
|
||||
r = evaluate_polynomial(P3, y) / evaluate_polynomial(Q3, y);
|
||||
factor = exp(-x) / sqrt(x);
|
||||
value = factor * r;
|
||||
template <typename T, class Policy>
|
||||
T bessel_k1_imp(const T& x, const Policy& pol, 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, 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, class Policy>
|
||||
T bessel_k1_imp(const T& x, const Policy& pol, 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, 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, class Policy>
|
||||
T bessel_k1_imp(const T& x, const Policy& pol, 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, 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;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return value;
|
||||
}
|
||||
template <typename T, class Policy>
|
||||
inline T bessel_k1(const T& x, const Policy& pol)
|
||||
{
|
||||
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, Policy, tag_type>::force_instantiate();
|
||||
return bessel_k1_imp(x, pol, tag_type());
|
||||
}
|
||||
|
||||
}}} // namespaces
|
||||
|
||||
|
||||
@@ -626,6 +626,17 @@ namespace boost
|
||||
bessel_maybe_int_tag
|
||||
>::type
|
||||
>::type optimisation_tag;
|
||||
typedef typename mpl::if_<
|
||||
mpl::or_<
|
||||
mpl::less_equal<precision_type, mpl::int_<0> >,
|
||||
mpl::greater<precision_type, mpl::int_<113> > >,
|
||||
bessel_no_int_tag,
|
||||
typename mpl::if_<
|
||||
is_integral<T1>,
|
||||
bessel_int_tag,
|
||||
bessel_maybe_int_tag
|
||||
>::type
|
||||
>::type optimisation_tag128;
|
||||
};
|
||||
} // detail
|
||||
|
||||
|
||||
@@ -334,6 +334,30 @@ mp_type f(const mp_type& x, int variant)
|
||||
mp_type xx = 1 / x;
|
||||
return boost::math::cyl_bessel_i(1, xx) * sqrt(xx) * exp(-xx);
|
||||
}
|
||||
case 36:
|
||||
{
|
||||
// K0 over [0, 1]
|
||||
mp_type xx = sqrt(x);
|
||||
return boost::math::cyl_bessel_k(0, xx) + log(xx) * boost::math::cyl_bessel_i(0, xx);
|
||||
}
|
||||
case 37:
|
||||
{
|
||||
// K0 over [1, INF]
|
||||
mp_type xx = 1 / x;
|
||||
return boost::math::cyl_bessel_k(0, xx) * exp(xx) * sqrt(xx);
|
||||
}
|
||||
case 38:
|
||||
{
|
||||
// K1 over [0, 1]
|
||||
mp_type xx = sqrt(x);
|
||||
return (boost::math::cyl_bessel_k(1, xx) - log(xx) * boost::math::cyl_bessel_i(1, xx) - 1 / xx) / xx;
|
||||
}
|
||||
case 39:
|
||||
{
|
||||
// K1 over [1, INF]
|
||||
mp_type xx = 1 / x;
|
||||
return boost::math::cyl_bessel_k(1, xx) * sqrt(xx) * exp(xx);
|
||||
}
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user