summaryrefslogtreecommitdiff
path: root/boost/math/special_functions/detail
diff options
context:
space:
mode:
Diffstat (limited to 'boost/math/special_functions/detail')
-rw-r--r--boost/math/special_functions/detail/bessel_i0.hpp30
-rw-r--r--boost/math/special_functions/detail/bessel_i1.hpp31
-rw-r--r--boost/math/special_functions/detail/bessel_j0.hpp30
-rw-r--r--boost/math/special_functions/detail/bessel_j1.hpp30
-rw-r--r--boost/math/special_functions/detail/bessel_jy.hpp16
-rw-r--r--boost/math/special_functions/detail/bessel_k0.hpp30
-rw-r--r--boost/math/special_functions/detail/bessel_k1.hpp30
-rw-r--r--boost/math/special_functions/detail/bessel_y0.hpp30
-rw-r--r--boost/math/special_functions/detail/bessel_y1.hpp30
-rw-r--r--boost/math/special_functions/detail/erf_inv.hpp40
-rw-r--r--boost/math/special_functions/detail/iconv.hpp2
-rw-r--r--boost/math/special_functions/detail/igamma_large.hpp1
-rw-r--r--boost/math/special_functions/detail/lanczos_sse2.hpp5
-rw-r--r--boost/math/special_functions/detail/lgamma_small.hpp8
14 files changed, 304 insertions, 9 deletions
diff --git a/boost/math/special_functions/detail/bessel_i0.hpp b/boost/math/special_functions/detail/bessel_i0.hpp
index 2c129facc7..7dc65d1a1b 100644
--- a/boost/math/special_functions/detail/bessel_i0.hpp
+++ b/boost/math/special_functions/detail/bessel_i0.hpp
@@ -21,8 +21,38 @@
namespace boost { namespace math { namespace detail{
template <typename T>
+T bessel_i0(T x);
+
+template <class T>
+struct bessel_i0_initializer
+{
+ struct init
+ {
+ init()
+ {
+ do_init();
+ }
+ static void do_init()
+ {
+ bessel_i0(T(1));
+ }
+ void force_instantiate()const{}
+ };
+ static const init initializer;
+ static void force_instantiate()
+ {
+ initializer.force_instantiate();
+ }
+};
+
+template <class T>
+const typename bessel_i0_initializer<T>::init bessel_i0_initializer<T>::initializer;
+
+template <typename T>
T bessel_i0(T x)
{
+ 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)),
diff --git a/boost/math/special_functions/detail/bessel_i1.hpp b/boost/math/special_functions/detail/bessel_i1.hpp
index aa4596cfd4..47f1b79883 100644
--- a/boost/math/special_functions/detail/bessel_i1.hpp
+++ b/boost/math/special_functions/detail/bessel_i1.hpp
@@ -21,8 +21,39 @@
namespace boost { namespace math { namespace detail{
template <typename T>
+T bessel_i1(T x);
+
+template <class T>
+struct bessel_i1_initializer
+{
+ struct init
+ {
+ init()
+ {
+ do_init();
+ }
+ static void do_init()
+ {
+ bessel_i1(T(1));
+ }
+ void force_instantiate()const{}
+ };
+ static const init initializer;
+ static void force_instantiate()
+ {
+ initializer.force_instantiate();
+ }
+};
+
+template <class T>
+const typename bessel_i1_initializer<T>::init bessel_i1_initializer<T>::initializer;
+
+template <typename T>
T bessel_i1(T x)
{
+
+ 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)),
diff --git a/boost/math/special_functions/detail/bessel_j0.hpp b/boost/math/special_functions/detail/bessel_j0.hpp
index ee25d46f61..a07052d73e 100644
--- a/boost/math/special_functions/detail/bessel_j0.hpp
+++ b/boost/math/special_functions/detail/bessel_j0.hpp
@@ -22,8 +22,38 @@
namespace boost { namespace math { namespace detail{
template <typename T>
+T bessel_j0(T x);
+
+template <class T>
+struct bessel_j0_initializer
+{
+ struct init
+ {
+ init()
+ {
+ do_init();
+ }
+ static void do_init()
+ {
+ bessel_j0(T(1));
+ }
+ void force_instantiate()const{}
+ };
+ static const init initializer;
+ static void force_instantiate()
+ {
+ initializer.force_instantiate();
+ }
+};
+
+template <class T>
+const typename bessel_j0_initializer<T>::init bessel_j0_initializer<T>::initializer;
+
+template <typename T>
T bessel_j0(T x)
{
+ bessel_j0_initializer<T>::force_instantiate();
+
static const T P1[] = {
static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -4.1298668500990866786e+11)),
static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.7282507878605942706e+10)),
diff --git a/boost/math/special_functions/detail/bessel_j1.hpp b/boost/math/special_functions/detail/bessel_j1.hpp
index 3db2503ff6..09d862c240 100644
--- a/boost/math/special_functions/detail/bessel_j1.hpp
+++ b/boost/math/special_functions/detail/bessel_j1.hpp
@@ -22,8 +22,38 @@
namespace boost { namespace math{ namespace detail{
template <typename T>
+T bessel_j1(T x);
+
+template <class T>
+struct bessel_j1_initializer
+{
+ struct init
+ {
+ init()
+ {
+ do_init();
+ }
+ static void do_init()
+ {
+ bessel_j1(T(1));
+ }
+ void force_instantiate()const{}
+ };
+ static const init initializer;
+ static void force_instantiate()
+ {
+ initializer.force_instantiate();
+ }
+};
+
+template <class T>
+const typename bessel_j1_initializer<T>::init bessel_j1_initializer<T>::initializer;
+
+template <typename T>
T bessel_j1(T x)
{
+ bessel_j1_initializer<T>::force_instantiate();
+
static const T P1[] = {
static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.4258509801366645672e+11)),
static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 6.6781041261492395835e+09)),
diff --git a/boost/math/special_functions/detail/bessel_jy.hpp b/boost/math/special_functions/detail/bessel_jy.hpp
index 19f951ab70..d60dda2d41 100644
--- a/boost/math/special_functions/detail/bessel_jy.hpp
+++ b/boost/math/special_functions/detail/bessel_jy.hpp
@@ -215,8 +215,6 @@ int CF2_jy(T v, T x, T* p, T* q, const Policy& pol)
Cr = br + a / temp;
Dr = br;
Di = bi;
- //std::cout << "C = " << Cr << " " << Ci << std::endl;
- //std::cout << "D = " << Dr << " " << Di << std::endl;
if (fabs(Cr) + fabs(Ci) < tiny) { Cr = tiny; }
if (fabs(Dr) + fabs(Di) < tiny) { Dr = tiny; }
temp = Dr * Dr + Di * Di;
@@ -227,7 +225,6 @@ int CF2_jy(T v, T x, T* p, T* q, const Policy& pol)
temp = fr;
fr = temp * delta_r - fi * delta_i;
fi = temp * delta_i + fi * delta_r;
- //std::cout << fr << " " << fi << std::endl;
for (k = 2; k < policies::get_max_series_iterations<Policy>(); k++)
{
a = k - 0.5f;
@@ -239,8 +236,6 @@ int CF2_jy(T v, T x, T* p, T* q, const Policy& pol)
Ci = bi - a * Ci / temp;
Dr = br + a * Dr;
Di = bi + a * Di;
- //std::cout << "C = " << Cr << " " << Ci << std::endl;
- //std::cout << "D = " << Dr << " " << Di << std::endl;
if (fabs(Cr) + fabs(Ci) < tiny) { Cr = tiny; }
if (fabs(Dr) + fabs(Di) < tiny) { Dr = tiny; }
temp = Dr * Dr + Di * Di;
@@ -253,7 +248,6 @@ int CF2_jy(T v, T x, T* p, T* q, const Policy& pol)
fi = temp * delta_i + fi * delta_r;
if (fabs(delta_r - 1) + fabs(delta_i) < tolerance)
break;
- //std::cout << fr << " " << fi << std::endl;
}
policies::check_series_iterations<T>("boost::math::bessel_jy<%1%>(%1%,%1%) in CF2_jy", k, pol);
*p = fr;
@@ -491,6 +485,16 @@ int bessel_jy(T v, T x, T* J, T* Y, int kind, const Policy& pol)
CF2_jy(u, x, &p, &q, pol); // continued fraction CF2_jy
T t = u / x - fu; // t = J'/J
gamma = (p - t) / q;
+ //
+ // We can't allow gamma to cancel out to zero competely as it messes up
+ // the subsequent logic. So pretend that one bit didn't cancel out
+ // and set to a suitably small value. The only test case we've been able to
+ // find for this, is when v = 8.5 and x = 4*PI.
+ //
+ if(gamma == 0)
+ {
+ gamma = u * tools::epsilon<T>() / x;
+ }
Ju = sign(current) * sqrt(W / (q + gamma * (p - t)));
Jv = Ju * ratio; // normalization
diff --git a/boost/math/special_functions/detail/bessel_k0.hpp b/boost/math/special_functions/detail/bessel_k0.hpp
index 81407dab10..e209168e87 100644
--- a/boost/math/special_functions/detail/bessel_k0.hpp
+++ b/boost/math/special_functions/detail/bessel_k0.hpp
@@ -22,10 +22,40 @@
namespace boost { namespace math { namespace detail{
template <typename T, typename Policy>
+T bessel_k0(T x, const Policy&);
+
+template <class T, class Policy>
+struct bessel_k0_initializer
+{
+ struct init
+ {
+ init()
+ {
+ do_init();
+ }
+ static void do_init()
+ {
+ bessel_k0(T(1), Policy());
+ }
+ void force_instantiate()const{}
+ };
+ static const init initializer;
+ static void force_instantiate()
+ {
+ initializer.force_instantiate();
+ }
+};
+
+template <class T, class Policy>
+const typename bessel_k0_initializer<T, Policy>::init bessel_k0_initializer<T, Policy>::initializer;
+
+template <typename T, typename Policy>
T bessel_k0(T x, const Policy& pol)
{
BOOST_MATH_INSTRUMENT_CODE(x);
+ bessel_k0_initializer<T, Policy>::force_instantiate();
+
static const T P1[] = {
static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.4708152720399552679e+03)),
static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.9169059852270512312e+03)),
diff --git a/boost/math/special_functions/detail/bessel_k1.hpp b/boost/math/special_functions/detail/bessel_k1.hpp
index 225209f7ba..0d17cd3057 100644
--- a/boost/math/special_functions/detail/bessel_k1.hpp
+++ b/boost/math/special_functions/detail/bessel_k1.hpp
@@ -22,8 +22,38 @@
namespace boost { namespace math { namespace detail{
template <typename T, typename Policy>
+T bessel_k1(T x, const Policy&);
+
+template <class T, class Policy>
+struct bessel_k1_initializer
+{
+ struct init
+ {
+ init()
+ {
+ do_init();
+ }
+ static void do_init()
+ {
+ bessel_k1(T(1), Policy());
+ }
+ void force_instantiate()const{}
+ };
+ static const init initializer;
+ static void force_instantiate()
+ {
+ initializer.force_instantiate();
+ }
+};
+
+template <class T, class Policy>
+const typename bessel_k1_initializer<T, Policy>::init bessel_k1_initializer<T, Policy>::initializer;
+
+template <typename T, typename Policy>
T bessel_k1(T x, const Policy& pol)
{
+ bessel_k1_initializer<T, Policy>::force_instantiate();
+
static const T P1[] = {
static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.2149374878243304548e+06)),
static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 7.1938920065420586101e+05)),
diff --git a/boost/math/special_functions/detail/bessel_y0.hpp b/boost/math/special_functions/detail/bessel_y0.hpp
index e23f861bf0..289bda5f18 100644
--- a/boost/math/special_functions/detail/bessel_y0.hpp
+++ b/boost/math/special_functions/detail/bessel_y0.hpp
@@ -24,8 +24,38 @@
namespace boost { namespace math { namespace detail{
template <typename T, typename Policy>
+T bessel_y0(T x, const Policy&);
+
+template <class T, class Policy>
+struct bessel_y0_initializer
+{
+ struct init
+ {
+ init()
+ {
+ do_init();
+ }
+ static void do_init()
+ {
+ bessel_y0(T(1), Policy());
+ }
+ void force_instantiate()const{}
+ };
+ static const init initializer;
+ static void force_instantiate()
+ {
+ initializer.force_instantiate();
+ }
+};
+
+template <class T, class Policy>
+const typename bessel_y0_initializer<T, Policy>::init bessel_y0_initializer<T, Policy>::initializer;
+
+template <typename T, typename Policy>
T bessel_y0(T x, const Policy& pol)
{
+ bessel_y0_initializer<T, Policy>::force_instantiate();
+
static const T P1[] = {
static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0723538782003176831e+11)),
static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -8.3716255451260504098e+09)),
diff --git a/boost/math/special_functions/detail/bessel_y1.hpp b/boost/math/special_functions/detail/bessel_y1.hpp
index b85e7011ea..caf09ffd26 100644
--- a/boost/math/special_functions/detail/bessel_y1.hpp
+++ b/boost/math/special_functions/detail/bessel_y1.hpp
@@ -24,8 +24,38 @@
namespace boost { namespace math { namespace detail{
template <typename T, typename Policy>
+T bessel_y1(T x, const Policy&);
+
+template <class T, class Policy>
+struct bessel_y1_initializer
+{
+ struct init
+ {
+ init()
+ {
+ do_init();
+ }
+ static void do_init()
+ {
+ bessel_y1(T(1), Policy());
+ }
+ void force_instantiate()const{}
+ };
+ static const init initializer;
+ static void force_instantiate()
+ {
+ initializer.force_instantiate();
+ }
+};
+
+template <class T, class Policy>
+const typename bessel_y1_initializer<T, Policy>::init bessel_y1_initializer<T, Policy>::initializer;
+
+template <typename T, typename Policy>
T bessel_y1(T x, const Policy& pol)
{
+ bessel_y1_initializer<T, Policy>::force_instantiate();
+
static const T P1[] = {
static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.0535726612579544093e+13)),
static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.4708611716525426053e+12)),
diff --git a/boost/math/special_functions/detail/erf_inv.hpp b/boost/math/special_functions/detail/erf_inv.hpp
index f2f625f991..d51db9d52f 100644
--- a/boost/math/special_functions/detail/erf_inv.hpp
+++ b/boost/math/special_functions/detail/erf_inv.hpp
@@ -322,12 +322,47 @@ T erf_inv_imp(const T& p, const T& q, const Policy& pol, const boost::mpl::int_<
return result;
}
+template <class T, class Policy>
+struct erf_inv_initializer
+{
+ struct init
+ {
+ init()
+ {
+ do_init();
+ }
+ static void do_init()
+ {
+ boost::math::erf_inv(static_cast<T>(0.25), Policy());
+ boost::math::erf_inv(static_cast<T>(0.55), Policy());
+ boost::math::erf_inv(static_cast<T>(0.95), Policy());
+ boost::math::erfc_inv(static_cast<T>(1e-15), Policy());
+ if(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-130)) != 0)
+ boost::math::erfc_inv(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-130)), Policy());
+ if(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-800)) != 0)
+ boost::math::erfc_inv(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-800)), Policy());
+ if(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-900)) != 0)
+ boost::math::erfc_inv(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-900)), Policy());
+ }
+ void force_instantiate()const{}
+ };
+ static const init initializer;
+ static void force_instantiate()
+ {
+ initializer.force_instantiate();
+ }
+};
+
+template <class T, class Policy>
+const typename erf_inv_initializer<T, Policy>::init erf_inv_initializer<T, Policy>::initializer;
+
} // namespace detail
template <class T, class Policy>
typename tools::promote_args<T>::type erfc_inv(T z, const Policy& pol)
{
typedef typename tools::promote_args<T>::type result_type;
+
//
// Begin by testing for domain errors, and other special cases:
//
@@ -378,6 +413,8 @@ typename tools::promote_args<T>::type erfc_inv(T z, const Policy& pol)
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
+ detail::erf_inv_initializer<eval_type, forwarding_policy>::force_instantiate();
+
//
// And get the result, negating where required:
//
@@ -389,6 +426,7 @@ template <class T, class Policy>
typename tools::promote_args<T>::type erf_inv(T z, const Policy& pol)
{
typedef typename tools::promote_args<T>::type result_type;
+
//
// Begin by testing for domain errors, and other special cases:
//
@@ -445,6 +483,8 @@ typename tools::promote_args<T>::type erf_inv(T z, const Policy& pol)
// precision internally if it's appropriate:
//
typedef typename policies::evaluation<result_type, Policy>::type eval_type;
+
+ detail::erf_inv_initializer<eval_type, forwarding_policy>::force_instantiate();
//
// And get the result, negating where required:
//
diff --git a/boost/math/special_functions/detail/iconv.hpp b/boost/math/special_functions/detail/iconv.hpp
index 8916eaed1d..4256ffcc88 100644
--- a/boost/math/special_functions/detail/iconv.hpp
+++ b/boost/math/special_functions/detail/iconv.hpp
@@ -25,7 +25,7 @@ template <class T, class Policy>
inline int iconv_imp(T v, Policy const& pol, mpl::false_ const&)
{
BOOST_MATH_STD_USING
- return iround(v);
+ return iround(v, pol);
}
template <class T, class Policy>
diff --git a/boost/math/special_functions/detail/igamma_large.hpp b/boost/math/special_functions/detail/igamma_large.hpp
index f9a810c489..eb3d4ba93e 100644
--- a/boost/math/special_functions/detail/igamma_large.hpp
+++ b/boost/math/special_functions/detail/igamma_large.hpp
@@ -759,7 +759,6 @@ T igamma_temme_large(T a, T x, const Policy& pol, mpl::int_<113> const *)
return result;
}
-
} // namespace detail
} // namespace math
} // namespace math
diff --git a/boost/math/special_functions/detail/lanczos_sse2.hpp b/boost/math/special_functions/detail/lanczos_sse2.hpp
index 6a3f3e5347..f8846bf376 100644
--- a/boost/math/special_functions/detail/lanczos_sse2.hpp
+++ b/boost/math/special_functions/detail/lanczos_sse2.hpp
@@ -13,7 +13,7 @@
#include <emmintrin.h>
#if defined(__GNUC__) || defined(__PGI)
-#define ALIGN16 __attribute__((aligned(16)))
+#define ALIGN16 __attribute__((__aligned__(16)))
#else
#define ALIGN16 __declspec(align(16))
#endif
@@ -194,8 +194,11 @@ inline double lanczos13m53::lanczos_sum_expG_scaled<double>(const double& x)
} // namespace math
} // namespace boost
+#undef ALIGN16
+
#endif // BOOST_MATH_SPECIAL_FUNCTIONS_LANCZOS
+
diff --git a/boost/math/special_functions/detail/lgamma_small.hpp b/boost/math/special_functions/detail/lgamma_small.hpp
index 526a573583..ec28ed2adf 100644
--- a/boost/math/special_functions/detail/lgamma_small.hpp
+++ b/boost/math/special_functions/detail/lgamma_small.hpp
@@ -15,6 +15,14 @@
namespace boost{ namespace math{ namespace detail{
//
+// These need forward declaring to keep GCC happy:
+//
+template <class T, class Policy, class Lanczos>
+T gamma_imp(T z, const Policy& pol, const Lanczos& l);
+template <class T, class Policy>
+T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos& l);
+
+//
// lgamma for small arguments:
//
template <class T, class Policy, class Lanczos>
generated by cgit v1.2.3 (git 2.25.1) at 2024-04-28 05:36:38 +0000