summaryrefslogtreecommitdiff
path: root/boost/math/special_functions/detail
diff options
context:
space:
mode:
authorDongHun Kwak <dh0128.kwak@samsung.com>2016-10-06 01:33:54 (GMT)
committerDongHun Kwak <dh0128.kwak@samsung.com>2016-10-06 01:36:09 (GMT)
commitd9ec475d945d3035377a0d89ed42e382d8988891 (patch)
tree34aff2cee4b209906243ab5499d61f3edee2982f /boost/math/special_functions/detail
parent71d216b90256936a9638f325af9bc69d720e75de (diff)
downloadboost-d9ec475d945d3035377a0d89ed42e382d8988891.zip
boost-d9ec475d945d3035377a0d89ed42e382d8988891.tar.gz
boost-d9ec475d945d3035377a0d89ed42e382d8988891.tar.bz2
Imported Upstream version 1.60.0refs/changes/07/91107/1
Change-Id: Ie709530d6d5841088ceaba025cbe175a4ef43050 Signed-off-by: DongHun Kwak <dh0128.kwak@samsung.com>
Diffstat (limited to 'boost/math/special_functions/detail')
-rw-r--r--boost/math/special_functions/detail/bernoulli_details.hpp18
-rw-r--r--boost/math/special_functions/detail/bessel_derivatives_linear.hpp16
-rw-r--r--boost/math/special_functions/detail/bessel_ik.hpp6
-rw-r--r--boost/math/special_functions/detail/bessel_jn.hpp7
-rw-r--r--boost/math/special_functions/detail/bessel_jy_asym.hpp18
-rw-r--r--boost/math/special_functions/detail/bessel_jy_derivatives_series.hpp5
-rw-r--r--boost/math/special_functions/detail/bessel_yn.hpp34
-rw-r--r--boost/math/special_functions/detail/lanczos_sse2.hpp2
-rw-r--r--boost/math/special_functions/detail/polygamma.hpp8
-rw-r--r--boost/math/special_functions/detail/t_distribution_inv.hpp40
10 files changed, 100 insertions, 54 deletions
diff --git a/boost/math/special_functions/detail/bernoulli_details.hpp b/boost/math/special_functions/detail/bernoulli_details.hpp
index 525c1fc..c12a172 100644
--- a/boost/math/special_functions/detail/bernoulli_details.hpp
+++ b/boost/math/special_functions/detail/bernoulli_details.hpp
@@ -202,9 +202,13 @@ struct bernoulli_initializer
// initialize our dymanic table:
//
boost::math::bernoulli_b2n<T>(2, Policy());
+#ifndef BOOST_NO_EXCEPTIONS
try{
+#endif
boost::math::bernoulli_b2n<T>(max_bernoulli_b2n<T>::value + 1, Policy());
+#ifndef BOOST_NO_EXCEPTIONS
} catch(const std::overflow_error&){}
+#endif
boost::math::tangent_t2n<T>(2, Policy());
}
void force_instantiate()const{}
@@ -254,7 +258,9 @@ struct fixed_vector : private std::allocator<T>
void resize(unsigned n, const T& val)
{
if(n > m_capacity)
- throw std::runtime_error("Exhausted storage for Bernoulli numbers.");
+ {
+ BOOST_THROW_EXCEPTION(std::runtime_error("Exhausted storage for Bernoulli numbers."));
+ }
for(unsigned i = m_used; i < n; ++i)
new (m_data + i) T(val);
m_used = n;
@@ -432,11 +438,11 @@ public:
//
if(start + n >= bn.size())
{
- std::size_t new_size = (std::min)((std::max)((std::max)(start + n, std::size_t(bn.size() + 20)), std::size_t(50)), std::size_t(bn.capacity()));
+ std::size_t new_size = (std::min)((std::max)((std::max)(std::size_t(start + n), std::size_t(bn.size() + 20)), std::size_t(50)), std::size_t(bn.capacity()));
tangent_numbers_series(new_size);
}
- for(std::size_t i = (std::max)(max_bernoulli_b2n<T>::value + 1, start); i < start + n; ++i)
+ for(std::size_t i = (std::max)(std::size_t(max_bernoulli_b2n<T>::value + 1), start); i < start + n; ++i)
{
*out = (i >= m_overflow_limit) ? policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol) : bn[i];
++out;
@@ -448,11 +454,11 @@ public:
boost::detail::lightweight_mutex::scoped_lock l(m_mutex);
if(start + n >= bn.size())
{
- std::size_t new_size = (std::min)((std::max)((std::max)(start + n, std::size_t(bn.size() + 20)), std::size_t(50)), std::size_t(bn.capacity()));
+ std::size_t new_size = (std::min)((std::max)((std::max)(std::size_t(start + n), std::size_t(bn.size() + 20)), std::size_t(50)), std::size_t(bn.capacity()));
tangent_numbers_series(new_size);
}
- for(std::size_t i = (std::max)(max_bernoulli_b2n<T>::value + 1, start); i < start + n; ++i)
+ for(std::size_t i = (std::max)(std::size_t(max_bernoulli_b2n<T>::value + 1), start); i < start + n; ++i)
{
*out = (i >= m_overflow_limit) ? policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol) : bn[i];
++out;
@@ -473,7 +479,7 @@ public:
{
if(start + n >= bn.size())
{
- std::size_t new_size = (std::min)((std::max)((std::max)(start + n, std::size_t(bn.size() + 20)), std::size_t(50)), std::size_t(bn.capacity()));
+ std::size_t new_size = (std::min)((std::max)((std::max)(std::size_t(start + n), std::size_t(bn.size() + 20)), std::size_t(50)), std::size_t(bn.capacity()));
tangent_numbers_series(new_size);
}
m_counter.store(static_cast<atomic_integer_type>(bn.size()), BOOST_MATH_ATOMIC_NS::memory_order_release);
diff --git a/boost/math/special_functions/detail/bessel_derivatives_linear.hpp b/boost/math/special_functions/detail/bessel_derivatives_linear.hpp
index 2ee86a0..72fddc6 100644
--- a/boost/math/special_functions/detail/bessel_derivatives_linear.hpp
+++ b/boost/math/special_functions/detail/bessel_derivatives_linear.hpp
@@ -37,13 +37,25 @@ inline T sph_bessel_j_derivative_linear(unsigned v, T x, Policy pol)
template <class T, class Policy>
inline T bessel_i_derivative_linear(T v, T x, Policy pol)
{
- return (boost::math::detail::cyl_bessel_i_imp<T>(v-1, x, pol) + boost::math::detail::cyl_bessel_i_imp<T>(v+1, x, pol)) / 2;
+ T result = boost::math::detail::cyl_bessel_i_imp<T>(v - 1, x, pol);
+ if(result >= tools::max_value<T>())
+ return result; // result is infinite
+ T result2 = boost::math::detail::cyl_bessel_i_imp<T>(v + 1, x, pol);
+ if(result2 >= tools::max_value<T>() - result)
+ return result2; // result is infinite
+ return (result + result2) / 2;
}
template <class T, class Tag, class Policy>
inline T bessel_k_derivative_linear(T v, T x, Tag tag, Policy pol)
{
- return (boost::math::detail::cyl_bessel_k_imp<T>(v-1, x, tag, pol) + boost::math::detail::cyl_bessel_k_imp<T>(v+1, x, tag, pol)) / -2;
+ T result = boost::math::detail::cyl_bessel_k_imp<T>(v - 1, x, tag, pol);
+ if(result >= tools::max_value<T>())
+ return -result; // result is infinite
+ T result2 = boost::math::detail::cyl_bessel_k_imp<T>(v + 1, x, tag, pol);
+ if(result2 >= tools::max_value<T>() - result)
+ return -result2; // result is infinite
+ return (result + result2) / -2;
}
template <class T, class Policy>
diff --git a/boost/math/special_functions/detail/bessel_ik.hpp b/boost/math/special_functions/detail/bessel_ik.hpp
index 10118d9..aac1781 100644
--- a/boost/math/special_functions/detail/bessel_ik.hpp
+++ b/boost/math/special_functions/detail/bessel_ik.hpp
@@ -374,6 +374,7 @@ int bessel_ik(T v, T x, T* I, T* K, int kind, const Policy& pol)
prev = Ku;
current = Ku1;
T scale = 1;
+ T scale_sign = 1;
for (k = 1; k <= n; k++) // forward recurrence for K
{
T fact = 2 * (u + k) / x;
@@ -381,6 +382,7 @@ int bessel_ik(T v, T x, T* I, T* K, int kind, const Policy& pol)
{
prev /= current;
scale /= current;
+ scale_sign *= boost::math::sign(current);
current = 1;
}
next = fact * current + prev;
@@ -426,7 +428,7 @@ int bessel_ik(T v, T x, T* I, T* K, int kind, const Policy& pol)
if(fact == 0)
*I = Iv;
else if(tools::max_value<T>() * scale < fact)
- *I = (org_kind & need_i) ? T(sign(fact) * sign(scale) * policies::raise_overflow_error<T>(function, 0, pol)) : T(0);
+ *I = (org_kind & need_i) ? T(sign(fact) * scale_sign * policies::raise_overflow_error<T>(function, 0, pol)) : T(0);
else
*I = Iv + fact / scale; // reflection formula
}
@@ -435,7 +437,7 @@ int bessel_ik(T v, T x, T* I, T* K, int kind, const Policy& pol)
*I = Iv;
}
if(tools::max_value<T>() * scale < Kv)
- *K = (org_kind & need_k) ? T(sign(Kv) * sign(scale) * policies::raise_overflow_error<T>(function, 0, pol)) : T(0);
+ *K = (org_kind & need_k) ? T(sign(Kv) * scale_sign * policies::raise_overflow_error<T>(function, 0, pol)) : T(0);
else
*K = Kv / scale;
BOOST_MATH_INSTRUMENT_VARIABLE(*I);
diff --git a/boost/math/special_functions/detail/bessel_jn.hpp b/boost/math/special_functions/detail/bessel_jn.hpp
index 3f15f9c..2413630 100644
--- a/boost/math/special_functions/detail/bessel_jn.hpp
+++ b/boost/math/special_functions/detail/bessel_jn.hpp
@@ -35,7 +35,7 @@ T bessel_jn(int n, T x, const Policy& pol)
//
if (n < 0)
{
- factor = (n & 0x1) ? -1 : 1; // J_{-n}(z) = (-1)^n J_n(z)
+ factor = static_cast<T>((n & 0x1) ? -1 : 1); // J_{-n}(z) = (-1)^n J_n(z)
n = -n;
}
else
@@ -50,6 +50,8 @@ T bessel_jn(int n, T x, const Policy& pol)
//
// Special cases:
//
+ if(asymptotic_bessel_large_x_limit(T(n), x))
+ return factor * asymptotic_bessel_j_large_x_2<T>(T(n), x);
if (n == 0)
{
return factor * bessel_j0(x);
@@ -64,9 +66,6 @@ T bessel_jn(int n, T x, const Policy& pol)
return static_cast<T>(0);
}
- if(asymptotic_bessel_large_x_limit(T(n), x))
- return factor * asymptotic_bessel_j_large_x_2<T>(n, x);
-
BOOST_ASSERT(n > 1);
T scale = 1;
if (n < abs(x)) // forward recurrence
diff --git a/boost/math/special_functions/detail/bessel_jy_asym.hpp b/boost/math/special_functions/detail/bessel_jy_asym.hpp
index 81f6238..4d7ac48 100644
--- a/boost/math/special_functions/detail/bessel_jy_asym.hpp
+++ b/boost/math/special_functions/detail/bessel_jy_asym.hpp
@@ -120,6 +120,24 @@ inline T asymptotic_bessel_j_large_x_2(T v, T x)
}
template <class T>
+inline bool asymptotic_bessel_large_x_limit(int v, const T& x)
+{
+ BOOST_MATH_STD_USING
+ //
+ // Determines if x is large enough compared to v to take the asymptotic
+ // forms above. From A&S 9.2.28 we require:
+ // v < x * eps^1/8
+ // and from A&S 9.2.29 we require:
+ // v^12/10 < 1.5 * x * eps^1/10
+ // using the former seems to work OK in practice with broadly similar
+ // error rates either side of the divide for v < 10000.
+ // At double precision eps^1/8 ~= 0.01.
+ //
+ BOOST_ASSERT(v >= 0);
+ return (v ? v : 1) < x * 0.004f;
+}
+
+template <class T>
inline bool asymptotic_bessel_large_x_limit(const T& v, const T& x)
{
BOOST_MATH_STD_USING
diff --git a/boost/math/special_functions/detail/bessel_jy_derivatives_series.hpp b/boost/math/special_functions/detail/bessel_jy_derivatives_series.hpp
index 0dc68fc..5425caf 100644
--- a/boost/math/special_functions/detail/bessel_jy_derivatives_series.hpp
+++ b/boost/math/special_functions/detail/bessel_jy_derivatives_series.hpp
@@ -157,7 +157,8 @@ inline T bessel_y_derivative_small_z_series(T v, T x, const Policy& pol)
gam = 1;
if (boost::math::tools::max_value<T>() * p < gam)
{
- return -boost::math::policies::raise_overflow_error<T>(function, 0, pol);
+ // This term will overflow to -INF, when combined with the series below it becomes +INF:
+ return boost::math::policies::raise_overflow_error<T>(function, 0, pol);
}
}
prefix = -gam / (boost::math::constants::pi<T>() * p);
@@ -173,7 +174,7 @@ inline T bessel_y_derivative_small_z_series(T v, T x, const Policy& pol)
scale /= (boost::math::tools::max_value<T>() / 4);
if (boost::math::tools::log_max_value<T>() < prefix)
{
- return -boost::math::policies::raise_overflow_error<T>(function, 0, pol);
+ return boost::math::policies::raise_overflow_error<T>(function, 0, pol);
}
}
prefix = -exp(prefix);
diff --git a/boost/math/special_functions/detail/bessel_yn.hpp b/boost/math/special_functions/detail/bessel_yn.hpp
index 0509062..62d7377 100644
--- a/boost/math/special_functions/detail/bessel_yn.hpp
+++ b/boost/math/special_functions/detail/bessel_yn.hpp
@@ -45,14 +45,13 @@ T bessel_yn(int n, T x, const Policy& pol)
//
if (n < 0)
{
- factor = (n & 0x1) ? -1 : 1; // Y_{-n}(z) = (-1)^n Y_n(z)
+ factor = static_cast<T>((n & 0x1) ? -1 : 1); // Y_{-n}(z) = (-1)^n Y_n(z)
n = -n;
}
else
{
factor = 1;
}
-
if(x < policies::get_epsilon<T, Policy>())
{
T scale = 1;
@@ -61,6 +60,10 @@ T bessel_yn(int n, T x, const Policy& pol)
return boost::math::sign(scale) * boost::math::sign(value) * policies::raise_overflow_error<T>(function, 0, pol);
value /= scale;
}
+ else if(asymptotic_bessel_large_x_limit(n, x))
+ {
+ value = factor * asymptotic_bessel_y_large_x_2(static_cast<T>(abs(n)), x);
+ }
else if (n == 0)
{
value = bessel_y0(x, pol);
@@ -76,23 +79,28 @@ T bessel_yn(int n, T x, const Policy& pol)
int k = 1;
BOOST_ASSERT(k < n);
policies::check_series_iterations<T>("boost::math::bessel_y_n<%1%>(%1%,%1%)", n, pol);
- do
+ T mult = 2 * k / x;
+ value = mult * current - prev;
+ prev = current;
+ current = value;
+ ++k;
+ if((mult > 1) && (fabs(current) > 1))
+ {
+ prev /= current;
+ factor /= current;
+ value /= current;
+ current = 1;
+ }
+ while(k < n)
{
- T fact = 2 * k / x;
- if((fact > 1) && ((tools::max_value<T>() - fabs(prev)) / fact < fabs(current)))
- {
- prev /= current;
- factor /= current;
- current = 1;
- }
- value = fact * current - prev;
+ mult = 2 * k / x;
+ value = mult * current - prev;
prev = current;
current = value;
++k;
}
- while(k < n);
if(fabs(tools::max_value<T>() * factor) < fabs(value))
- return sign(value) * sign(value) * policies::raise_overflow_error<T>(function, 0, pol);
+ return sign(value) * sign(factor) * policies::raise_overflow_error<T>(function, 0, pol);
value /= factor;
}
return value;
diff --git a/boost/math/special_functions/detail/lanczos_sse2.hpp b/boost/math/special_functions/detail/lanczos_sse2.hpp
index edef3a0..b2dd7ea 100644
--- a/boost/math/special_functions/detail/lanczos_sse2.hpp
+++ b/boost/math/special_functions/detail/lanczos_sse2.hpp
@@ -12,7 +12,7 @@
#include <emmintrin.h>
-#if defined(__GNUC__) || defined(__PGI)
+#if defined(__GNUC__) || defined(__PGI) || defined(__SUNPRO_CC)
#define ALIGN16 __attribute__((__aligned__(16)))
#else
#define ALIGN16 __declspec(align(16))
diff --git a/boost/math/special_functions/detail/polygamma.hpp b/boost/math/special_functions/detail/polygamma.hpp
index 4ef503b..0267bd3 100644
--- a/boost/math/special_functions/detail/polygamma.hpp
+++ b/boost/math/special_functions/detail/polygamma.hpp
@@ -44,7 +44,7 @@
// x is crazy large, just concentrate on the first part of the expression and use logs:
if(n == 1) return 1 / x;
T nlx = n * log(x);
- if((nlx < tools::log_max_value<T>()) && (n < max_factorial<T>::value))
+ if((nlx < tools::log_max_value<T>()) && (n < (int)max_factorial<T>::value))
return ((n & 1) ? 1 : -1) * boost::math::factorial<T>(n - 1) * pow(x, -n);
else
return ((n & 1) ? 1 : -1) * exp(boost::math::lgamma(T(n), pol) - n * log(x));
@@ -71,13 +71,13 @@
// or the power term underflows, this just gets set to 0 and then we
// know that we have to use logs for the initial terms:
//
- part_term = ((n > boost::math::max_factorial<T>::value) && (T(n) * n > tools::log_max_value<T>()))
+ part_term = ((n > (int)boost::math::max_factorial<T>::value) && (T(n) * n > tools::log_max_value<T>()))
? T(0) : static_cast<T>(boost::math::factorial<T>(n - 1, pol) * pow(x, -n - 1));
if(part_term == 0)
{
// Either n is very large, or the power term underflows,
// set the initial values of part_term, term and sum via logs:
- part_term = boost::math::lgamma(n, pol) - (n + 1) * log(x);
+ part_term = static_cast<T>(boost::math::lgamma(n, pol) - (n + 1) * log(x));
sum = exp(part_term + log(n + 2 * x) - boost::math::constants::ln_two<T>());
part_term += log(T(n) * (n + 1)) - boost::math::constants::ln_two<T>() - log(x);
part_term = exp(part_term);
@@ -505,7 +505,7 @@
// would mean setting the limit to ~ 1 / n,
// but we can tolerate a small amount of divergence:
//
- T small_x_limit = std::min(T(T(5) / n), T(0.25f));
+ T small_x_limit = (std::min)(T(T(5) / n), T(0.25f));
if(x < small_x_limit)
{
return polygamma_nearzero(n, x, pol, function);
diff --git a/boost/math/special_functions/detail/t_distribution_inv.hpp b/boost/math/special_functions/detail/t_distribution_inv.hpp
index 72f6f0c..e8bd0db 100644
--- a/boost/math/special_functions/detail/t_distribution_inv.hpp
+++ b/boost/math/special_functions/detail/t_distribution_inv.hpp
@@ -56,9 +56,9 @@ T inverse_students_t_hill(T ndf, T u, const Policy& pol)
}
else
{
- y = ((1 / (((ndf + 6) / (ndf * y) - 0.089f * d - 0.822f)
+ y = static_cast<T>(((1 / (((ndf + 6) / (ndf * y) - 0.089f * d - 0.822f)
* (ndf + 2) * 3) + 0.5 / (ndf + 4)) * y - 1)
- * (ndf + 1) / (ndf + 2) + 1 / y;
+ * (ndf + 1) / (ndf + 2) + 1 / y);
}
q = sqrt(ndf * y);
@@ -144,41 +144,41 @@ T inverse_students_t_body_series(T df, T u, const Policy& pol)
// only on the degrees of freedom (Eq 57 of Shaw):
//
T in = 1 / df;
- c[2] = 0.16666666666666666667 + 0.16666666666666666667 * in;
- c[3] = (0.0083333333333333333333 * in
+ c[2] = static_cast<T>(0.16666666666666666667 + 0.16666666666666666667 * in);
+ c[3] = static_cast<T>((0.0083333333333333333333 * in
+ 0.066666666666666666667) * in
- + 0.058333333333333333333;
- c[4] = ((0.00019841269841269841270 * in
+ + 0.058333333333333333333);
+ c[4] = static_cast<T>(((0.00019841269841269841270 * in
+ 0.0017857142857142857143) * in
+ 0.026785714285714285714) * in
- + 0.025198412698412698413;
- c[5] = (((2.7557319223985890653e-6 * in
+ + 0.025198412698412698413);
+ c[5] = static_cast<T>((((2.7557319223985890653e-6 * in
+ 0.00037477954144620811287) * in
- 0.0011078042328042328042) * in
+ 0.010559964726631393298) * in
- + 0.012039792768959435626;
- c[6] = ((((2.5052108385441718775e-8 * in
+ + 0.012039792768959435626);
+ c[6] = static_cast<T>(((((2.5052108385441718775e-8 * in
- 0.000062705427288760622094) * in
+ 0.00059458674042007375341) * in
- 0.0016095979637646304313) * in
+ 0.0061039211560044893378) * in
- + 0.0038370059724226390893;
- c[7] = (((((1.6059043836821614599e-10 * in
+ + 0.0038370059724226390893);
+ c[7] = static_cast<T>((((((1.6059043836821614599e-10 * in
+ 0.000015401265401265401265) * in
- 0.00016376804137220803887) * in
+ 0.00069084207973096861986) * in
- 0.0012579159844784844785) * in
+ 0.0010898206731540064873) * in
- + 0.0032177478835464946576;
- c[8] = ((((((7.6471637318198164759e-13 * in
+ + 0.0032177478835464946576);
+ c[8] = static_cast<T>(((((((7.6471637318198164759e-13 * in
- 3.9851014346715404916e-6) * in
+ 0.000049255746366361445727) * in
- 0.00024947258047043099953) * in
+ 0.00064513046951456342991) * in
- 0.00076245135440323932387) * in
+ 0.000033530976880017885309) * in
- + 0.0017438262298340009980;
- c[9] = (((((((2.8114572543455207632e-15 * in
+ + 0.0017438262298340009980);
+ c[9] = static_cast<T>((((((((2.8114572543455207632e-15 * in
+ 1.0914179173496789432e-6) * in
- 0.000015303004486655377567) * in
+ 0.000090867107935219902229) * in
@@ -186,8 +186,8 @@ T inverse_students_t_body_series(T df, T u, const Policy& pol)
+ 0.00051406605788341121363) * in
- 0.00036307660358786885787) * in
- 0.00031101086326318780412) * in
- + 0.00096472747321388644237;
- c[10] = ((((((((8.2206352466243297170e-18 * in
+ + 0.00096472747321388644237);
+ c[10] = static_cast<T>(((((((((8.2206352466243297170e-18 * in
- 3.1239569599829868045e-7) * in
+ 4.8903045291975346210e-6) * in
- 0.000033202652391372058698) * in
@@ -196,7 +196,7 @@ T inverse_students_t_body_series(T df, T u, const Policy& pol)
+ 0.00035764655430568632777) * in
- 0.00010230378073700412687) * in
- 0.00036942667800009661203) * in
- + 0.00054229262813129686486;
+ + 0.00054229262813129686486);
//
// The result is then a polynomial in v (see Eq 56 of Shaw):
//
@@ -285,7 +285,7 @@ T inverse_students_t(T df, T u, T v, const Policy& pol, bool* pexact = 0)
//
T a = 4 * (u - u * u);//1 - 4 * (u - 0.5f) * (u - 0.5f);
T b = boost::math::cbrt(a);
- static const T c = 0.85498797333834849467655443627193;
+ static const T c = static_cast<T>(0.85498797333834849467655443627193L);
T p = 6 * (1 + c * (1 / b - 1));
T p0;
do{