summaryrefslogtreecommitdiff
path: root/boost/math/special_functions
diff options
context:
space:
mode:
Diffstat (limited to 'boost/math/special_functions')
-rw-r--r--boost/math/special_functions/detail/airy_ai_bi_zero.hpp4
-rw-r--r--boost/math/special_functions/detail/bessel_i0.hpp2
-rw-r--r--boost/math/special_functions/detail/bessel_i1.hpp2
-rw-r--r--boost/math/special_functions/detail/bessel_k0.hpp2
-rw-r--r--boost/math/special_functions/detail/bessel_k1.hpp2
-rw-r--r--boost/math/special_functions/legendre.hpp191
-rw-r--r--boost/math/special_functions/legendre_stieltjes.hpp235
-rw-r--r--boost/math/special_functions/math_fwd.hpp21
-rw-r--r--boost/math/special_functions/next.hpp384
-rw-r--r--boost/math/special_functions/sign.hpp2
-rw-r--r--boost/math/special_functions/ulp.hpp38
11 files changed, 853 insertions, 30 deletions
diff --git a/boost/math/special_functions/detail/airy_ai_bi_zero.hpp b/boost/math/special_functions/detail/airy_ai_bi_zero.hpp
index dbb7388dda..b5a71c78eb 100644
--- a/boost/math/special_functions/detail/airy_ai_bi_zero.hpp
+++ b/boost/math/special_functions/detail/airy_ai_bi_zero.hpp
@@ -86,7 +86,7 @@
class function_object_ai_and_ai_prime
{
public:
- function_object_ai_and_ai_prime(const Policy pol) : my_pol(pol) { }
+ function_object_ai_and_ai_prime(const Policy& pol) : my_pol(pol) { }
boost::math::tuple<T, T> operator()(const T& x) const
{
@@ -137,7 +137,7 @@
class function_object_bi_and_bi_prime
{
public:
- function_object_bi_and_bi_prime(const Policy pol) : my_pol(pol) { }
+ function_object_bi_and_bi_prime(const Policy& pol) : my_pol(pol) { }
boost::math::tuple<T, T> operator()(const T& x) const
{
diff --git a/boost/math/special_functions/detail/bessel_i0.hpp b/boost/math/special_functions/detail/bessel_i0.hpp
index 5896db8f76..7738229454 100644
--- a/boost/math/special_functions/detail/bessel_i0.hpp
+++ b/boost/math/special_functions/detail/bessel_i0.hpp
@@ -532,7 +532,7 @@ template <typename T>
inline T bessel_i0(const T& x)
{
typedef mpl::int_<
- std::numeric_limits<T>::digits == 0 ?
+ ((std::numeric_limits<T>::digits == 0) || (std::numeric_limits<T>::radix != 2)) ?
0 :
std::numeric_limits<T>::digits <= 24 ?
24 :
diff --git a/boost/math/special_functions/detail/bessel_i1.hpp b/boost/math/special_functions/detail/bessel_i1.hpp
index 0f9ba96bd8..a6364a6cc9 100644
--- a/boost/math/special_functions/detail/bessel_i1.hpp
+++ b/boost/math/special_functions/detail/bessel_i1.hpp
@@ -555,7 +555,7 @@ template <typename T>
inline T bessel_i1(const T& x)
{
typedef mpl::int_<
- std::numeric_limits<T>::digits == 0 ?
+ ((std::numeric_limits<T>::digits == 0) || (std::numeric_limits<T>::radix != 2)) ?
0 :
std::numeric_limits<T>::digits <= 24 ?
24 :
diff --git a/boost/math/special_functions/detail/bessel_k0.hpp b/boost/math/special_functions/detail/bessel_k0.hpp
index 4f3420ffc0..74f4014bd9 100644
--- a/boost/math/special_functions/detail/bessel_k0.hpp
+++ b/boost/math/special_functions/detail/bessel_k0.hpp
@@ -483,7 +483,7 @@ template <typename T>
inline T bessel_k0(const T& x)
{
typedef mpl::int_<
- std::numeric_limits<T>::digits == 0 ?
+ ((std::numeric_limits<T>::digits == 0) || (std::numeric_limits<T>::radix != 2)) ?
0 :
std::numeric_limits<T>::digits <= 24 ?
24 :
diff --git a/boost/math/special_functions/detail/bessel_k1.hpp b/boost/math/special_functions/detail/bessel_k1.hpp
index 0ae90056cb..2ab191fb49 100644
--- a/boost/math/special_functions/detail/bessel_k1.hpp
+++ b/boost/math/special_functions/detail/bessel_k1.hpp
@@ -525,7 +525,7 @@ namespace boost { namespace math { namespace detail{
inline T bessel_k1(const T& x)
{
typedef mpl::int_<
- std::numeric_limits<T>::digits == 0 ?
+ ((std::numeric_limits<T>::digits == 0) || (std::numeric_limits<T>::radix != 2)) ?
0 :
std::numeric_limits<T>::digits <= 24 ?
24 :
diff --git a/boost/math/special_functions/legendre.hpp b/boost/math/special_functions/legendre.hpp
index 1a2ef5d615..6028b377d3 100644
--- a/boost/math/special_functions/legendre.hpp
+++ b/boost/math/special_functions/legendre.hpp
@@ -11,8 +11,11 @@
#pragma once
#endif
+#include <utility>
+#include <vector>
#include <boost/math/special_functions/math_fwd.hpp>
#include <boost/math/special_functions/factorials.hpp>
+#include <boost/math/tools/roots.hpp>
#include <boost/math/tools/config.hpp>
namespace boost{
@@ -68,6 +71,149 @@ T legendre_imp(unsigned l, T x, const Policy& pol, bool second = false)
return p1;
}
+template <class T, class Policy>
+T legendre_p_prime_imp(unsigned l, T x, const Policy& pol, T* Pn
+#ifdef BOOST_NO_CXX11_NULLPTR
+ = 0
+#else
+ = nullptr
+#endif
+)
+{
+ static const char* function = "boost::math::legrendre_p_prime<%1%>(unsigned, %1%)";
+ // Error handling:
+ if ((x < -1) || (x > 1))
+ return policies::raise_domain_error<T>(
+ function,
+ "The Legendre Polynomial is defined for"
+ " -1 <= x <= 1, but got x = %1%.", x, pol);
+
+ if (l == 0)
+ {
+ if (Pn)
+ {
+ *Pn = 1;
+ }
+ return 0;
+ }
+ T p0 = 1;
+ T p1 = x;
+ T p_prime;
+ bool odd = l & 1;
+ // If the order is odd, we sum all the even polynomials:
+ if (odd)
+ {
+ p_prime = p0;
+ }
+ else // Otherwise we sum the odd polynomials * (2n+1)
+ {
+ p_prime = 3*p1;
+ }
+
+ unsigned n = 1;
+ while(n < l - 1)
+ {
+ std::swap(p0, p1);
+ p1 = boost::math::legendre_next(n, x, p0, p1);
+ ++n;
+ if (odd)
+ {
+ p_prime += (2*n+1)*p1;
+ odd = false;
+ }
+ else
+ {
+ odd = true;
+ }
+ }
+ // This allows us to evaluate the derivative and the function for the same cost.
+ if (Pn)
+ {
+ std::swap(p0, p1);
+ *Pn = boost::math::legendre_next(n, x, p0, p1);
+ }
+ return p_prime;
+}
+
+template <class T, class Policy>
+struct legendre_p_zero_func
+{
+ int n;
+ const Policy& pol;
+
+ legendre_p_zero_func(int n_, const Policy& p) : n(n_), pol(p) {}
+
+ std::pair<T, T> operator()(T x) const
+ {
+ T Pn;
+ T Pn_prime = detail::legendre_p_prime_imp(n, x, pol, &Pn);
+ return std::pair<T, T>(Pn, Pn_prime);
+ };
+};
+
+template <class T, class Policy>
+std::vector<T> legendre_p_zeros_imp(int n, const Policy& pol)
+{
+ using std::cos;
+ using std::sin;
+ using std::ceil;
+ using std::sqrt;
+ using boost::math::constants::pi;
+ using boost::math::constants::half;
+ using boost::math::tools::newton_raphson_iterate;
+
+ BOOST_ASSERT(n >= 0);
+ std::vector<T> zeros;
+ if (n == 0)
+ {
+ // There are no zeros of P_0(x) = 1.
+ return zeros;
+ }
+ int k;
+ if (n & 1)
+ {
+ zeros.resize((n-1)/2 + 1, std::numeric_limits<T>::quiet_NaN());
+ zeros[0] = 0;
+ k = 1;
+ }
+ else
+ {
+ zeros.resize(n/2, std::numeric_limits<T>::quiet_NaN());
+ k = 0;
+ }
+ T half_n = ceil(n*half<T>());
+
+ while (k < (int)zeros.size())
+ {
+ // Bracket the root: Szego:
+ // Gabriel Szego, Inequalities for the Zeros of Legendre Polynomials and Related Functions, Transactions of the American Mathematical Society, Vol. 39, No. 1 (1936)
+ T theta_nk = ((half_n - half<T>()*half<T>() - static_cast<T>(k))*pi<T>())/(static_cast<T>(n)+half<T>());
+ T lower_bound = cos( (half_n - static_cast<T>(k))*pi<T>()/static_cast<T>(n + 1));
+ T cos_nk = cos(theta_nk);
+ T upper_bound = cos_nk;
+ // First guess follows from:
+ // F. G. Tricomi, Sugli zeri dei polinomi sferici ed ultrasferici, Ann. Mat. Pura Appl., 31 (1950), pp. 93–97;
+ T inv_n_sq = 1/static_cast<T>(n*n);
+ T sin_nk = sin(theta_nk);
+ T x_nk_guess = (1 - inv_n_sq/static_cast<T>(8) + inv_n_sq /static_cast<T>(8*n) - (inv_n_sq*inv_n_sq/384)*(39 - 28 / (sin_nk*sin_nk) ) )*cos_nk;
+
+ boost::uintmax_t number_of_iterations = policies::get_max_root_iterations<Policy>();
+
+ legendre_p_zero_func<T, Policy> f(n, pol);
+
+ const T x_nk = newton_raphson_iterate(f, x_nk_guess,
+ lower_bound, upper_bound,
+ policies::digits<T, Policy>(),
+ number_of_iterations);
+
+ BOOST_ASSERT(lower_bound < x_nk);
+ BOOST_ASSERT(upper_bound > x_nk);
+ zeros[k] = x_nk;
+ ++k;
+ }
+ return zeros;
+}
+
} // namespace detail
template <class T, class Policy>
@@ -82,13 +228,49 @@ inline typename boost::enable_if_c<policies::is_policy<Policy>::value, typename
return policies::checked_narrowing_cast<result_type, Policy>(detail::legendre_imp(l, static_cast<value_type>(x), pol, false), function);
}
+
+template <class T, class Policy>
+inline typename boost::enable_if_c<policies::is_policy<Policy>::value, typename tools::promote_args<T>::type>::type
+ legendre_p_prime(int l, T x, const Policy& pol)
+{
+ typedef typename tools::promote_args<T>::type result_type;
+ typedef typename policies::evaluation<result_type, Policy>::type value_type;
+ static const char* function = "boost::math::legendre_p_prime<%1%>(unsigned, %1%)";
+ if(l < 0)
+ return policies::checked_narrowing_cast<result_type, Policy>(detail::legendre_p_prime_imp(-l-1, static_cast<value_type>(x), pol), function);
+ return policies::checked_narrowing_cast<result_type, Policy>(detail::legendre_p_prime_imp(l, static_cast<value_type>(x), pol), function);
+}
+
template <class T>
-inline typename tools::promote_args<T>::type
+inline typename tools::promote_args<T>::type
legendre_p(int l, T x)
{
return boost::math::legendre_p(l, x, policies::policy<>());
}
+template <class T>
+inline typename tools::promote_args<T>::type
+ legendre_p_prime(int l, T x)
+{
+ return boost::math::legendre_p_prime(l, x, policies::policy<>());
+}
+
+template <class T, class Policy>
+inline std::vector<T> legendre_p_zeros(int l, const Policy& pol)
+{
+ if(l < 0)
+ return detail::legendre_p_zeros_imp<T>(-l-1, pol);
+
+ return detail::legendre_p_zeros_imp<T>(l, pol);
+}
+
+
+template <class T>
+inline std::vector<T> legendre_p_zeros(int l)
+{
+ return boost::math::legendre_p_zeros<T>(l, policies::policy<>());
+}
+
template <class T, class Policy>
inline typename boost::enable_if_c<policies::is_policy<Policy>::value, typename tools::promote_args<T>::type>::type
legendre_q(unsigned l, T x, const Policy& pol)
@@ -99,7 +281,7 @@ inline typename boost::enable_if_c<policies::is_policy<Policy>::value, typename
}
template <class T>
-inline typename tools::promote_args<T>::type
+inline typename tools::promote_args<T>::type
legendre_q(unsigned l, T x)
{
return boost::math::legendre_q(l, x, policies::policy<>());
@@ -107,7 +289,7 @@ inline typename tools::promote_args<T>::type
// Recurrence for associated polynomials:
template <class T1, class T2, class T3>
-inline typename tools::promote_args<T1, T2, T3>::type
+inline typename tools::promote_args<T1, T2, T3>::type
legendre_next(unsigned l, unsigned m, T1 x, T2 Pl, T3 Plm1)
{
typedef typename tools::promote_args<T1, T2, T3>::type result_type;
@@ -189,6 +371,3 @@ inline typename tools::promote_args<T>::type
} // namespace boost
#endif // BOOST_MATH_SPECIAL_LEGENDRE_HPP
-
-
-
diff --git a/boost/math/special_functions/legendre_stieltjes.hpp b/boost/math/special_functions/legendre_stieltjes.hpp
new file mode 100644
index 0000000000..5d68198870
--- /dev/null
+++ b/boost/math/special_functions/legendre_stieltjes.hpp
@@ -0,0 +1,235 @@
+// Copyright Nick Thompson 2017.
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_MATH_SPECIAL_LEGENDRE_STIELTJES_HPP
+#define BOOST_MATH_SPECIAL_LEGENDRE_STIELTJES_HPP
+
+/*
+ * Constructs the Legendre-Stieltjes polynomial of degree m.
+ * The Legendre-Stieltjes polynomials are used to create extensions for Gaussian quadratures,
+ * commonly called "Gauss-Konrod" quadratures.
+ *
+ * References:
+ * Patterson, TNL. "The optimum addition of points to quadrature formulae." Mathematics of Computation 22.104 (1968): 847-856.
+ */
+
+#include <iostream>
+#include <vector>
+#include <boost/math/tools/roots.hpp>
+#include <boost/math/special_functions/legendre.hpp>
+
+namespace boost{
+namespace math{
+
+template<class Real>
+class legendre_stieltjes
+{
+public:
+ legendre_stieltjes(size_t m)
+ {
+ if (m == 0)
+ {
+ throw std::domain_error("The Legendre-Stieltjes polynomial is defined for order m > 0.\n");
+ }
+ m_m = m;
+ std::ptrdiff_t n = m - 1;
+ std::ptrdiff_t q;
+ std::ptrdiff_t r;
+ bool odd = n & 1;
+ if (odd)
+ {
+ q = 1;
+ r = (n-1)/2 + 2;
+ }
+ else
+ {
+ q = 0;
+ r = n/2 + 1;
+ }
+ m_a.resize(r + 1);
+ // We'll keep the ones-based indexing at the cost of storing a superfluous element
+ // so that we can follow Patterson's notation exactly.
+ m_a[r] = static_cast<Real>(1);
+ // Make sure using the zero index is a bug:
+ m_a[0] = std::numeric_limits<Real>::quiet_NaN();
+
+ for (std::ptrdiff_t k = 1; k < r; ++k)
+ {
+ Real ratio = 1;
+ m_a[r - k] = 0;
+ for (std::ptrdiff_t i = r + 1 - k; i <= r; ++i)
+ {
+ // See Patterson, equation 12
+ std::ptrdiff_t num = (n - q + 2*(i + k - 1))*(n + q + 2*(k - i + 1))*(n-1-q+2*(i-k))*(2*(k+i-1) -1 -q -n);
+ std::ptrdiff_t den = (n - q + 2*(i - k))*(2*(k + i - 1) - q - n)*(n + 1 + q + 2*(k - i))*(n - 1 - q + 2*(i + k));
+ ratio *= static_cast<Real>(num)/static_cast<Real>(den);
+ m_a[r - k] -= ratio*m_a[i];
+ }
+ }
+ }
+
+
+ Real norm_sq() const
+ {
+ Real t = 0;
+ bool odd = m_m & 1;
+ for (size_t i = 1; i < m_a.size(); ++i)
+ {
+ if(odd)
+ {
+ t += 2*m_a[i]*m_a[i]/static_cast<Real>(4*i-1);
+ }
+ else
+ {
+ t += 2*m_a[i]*m_a[i]/static_cast<Real>(4*i-3);
+ }
+ }
+ return t;
+ }
+
+
+ Real operator()(Real x) const
+ {
+ // Trivial implementation:
+ // Em += m_a[i]*legendre_p(2*i - 1, x); m odd
+ // Em += m_a[i]*legendre_p(2*i - 2, x); m even
+ size_t r = m_a.size() - 1;
+ Real p0 = 1;
+ Real p1 = x;
+
+ Real Em;
+ bool odd = m_m & 1;
+ if (odd)
+ {
+ Em = m_a[1]*p1;
+ }
+ else
+ {
+ Em = m_a[1]*p0;
+ }
+
+ unsigned n = 1;
+ for (size_t i = 2; i <= r; ++i)
+ {
+ std::swap(p0, p1);
+ p1 = boost::math::legendre_next(n, x, p0, p1);
+ ++n;
+ if (!odd)
+ {
+ Em += m_a[i]*p1;
+ }
+ std::swap(p0, p1);
+ p1 = boost::math::legendre_next(n, x, p0, p1);
+ ++n;
+ if(odd)
+ {
+ Em += m_a[i]*p1;
+ }
+ }
+ return Em;
+ }
+
+
+ Real prime(Real x) const
+ {
+ Real Em_prime = 0;
+
+ for (size_t i = 1; i < m_a.size(); ++i)
+ {
+ if(m_m & 1)
+ {
+ Em_prime += m_a[i]*detail::legendre_p_prime_imp(2*i - 1, x, policies::policy<>());
+ }
+ else
+ {
+ Em_prime += m_a[i]*detail::legendre_p_prime_imp(2*i - 2, x, policies::policy<>());
+ }
+ }
+ return Em_prime;
+ }
+
+ std::vector<Real> zeros() const
+ {
+ using boost::math::constants::half;
+
+ std::vector<Real> stieltjes_zeros;
+ std::vector<Real> legendre_zeros = legendre_p_zeros<Real>(m_m - 1);
+ int k;
+ if (m_m & 1)
+ {
+ stieltjes_zeros.resize(legendre_zeros.size() + 1, std::numeric_limits<Real>::quiet_NaN());
+ stieltjes_zeros[0] = 0;
+ k = 1;
+ }
+ else
+ {
+ stieltjes_zeros.resize(legendre_zeros.size(), std::numeric_limits<Real>::quiet_NaN());
+ k = 0;
+ }
+
+ while (k < (int)stieltjes_zeros.size())
+ {
+ Real lower_bound;
+ Real upper_bound;
+ if (m_m & 1)
+ {
+ lower_bound = legendre_zeros[k - 1];
+ if (k == (int)legendre_zeros.size())
+ {
+ upper_bound = 1;
+ }
+ else
+ {
+ upper_bound = legendre_zeros[k];
+ }
+ }
+ else
+ {
+ lower_bound = legendre_zeros[k];
+ if (k == (int)legendre_zeros.size() - 1)
+ {
+ upper_bound = 1;
+ }
+ else
+ {
+ upper_bound = legendre_zeros[k+1];
+ }
+ }
+
+ // The root bracketing is not very tight; to keep weird stuff from happening
+ // in the Newton's method, let's tighten up the tolerance using a few bisections.
+ boost::math::tools::eps_tolerance<Real> tol(6);
+ auto g = [&](Real t) { return this->operator()(t); };
+ auto p = boost::math::tools::bisect(g, lower_bound, upper_bound, tol);
+
+ Real x_nk_guess = p.first + (p.second - p.first)*half<Real>();
+ boost::uintmax_t number_of_iterations = 500;
+
+ auto f = [&] (Real x) { Real Pn = this->operator()(x);
+ Real Pn_prime = this->prime(x);
+ return std::pair<Real, Real>(Pn, Pn_prime); };
+
+ const Real x_nk = boost::math::tools::newton_raphson_iterate(f, x_nk_guess,
+ p.first, p.second,
+ 2*std::numeric_limits<Real>::digits10,
+ number_of_iterations);
+
+ BOOST_ASSERT(p.first < x_nk);
+ BOOST_ASSERT(x_nk < p.second);
+ stieltjes_zeros[k] = x_nk;
+ ++k;
+ }
+ return stieltjes_zeros;
+ }
+
+private:
+ // Coefficients of Legendre expansion
+ std::vector<Real> m_a;
+ int m_m;
+};
+
+}}
+#endif
diff --git a/boost/math/special_functions/math_fwd.hpp b/boost/math/special_functions/math_fwd.hpp
index efb8f067ac..9becef9f6d 100644
--- a/boost/math/special_functions/math_fwd.hpp
+++ b/boost/math/special_functions/math_fwd.hpp
@@ -23,6 +23,7 @@
#pragma once
#endif
+#include <vector>
#include <boost/math/special_functions/detail/round_fwd.hpp>
#include <boost/math/tools/promotion.hpp> // for argument promotion.
#include <boost/math/policies/policy.hpp>
@@ -181,10 +182,24 @@ namespace boost
template <class T>
typename tools::promote_args<T>::type
legendre_p(int l, T x);
+ template <class T>
+ typename tools::promote_args<T>::type
+ legendre_p_prime(int l, T x);
+
+
+ template <class T, class Policy>
+ inline std::vector<T> legendre_p_zeros(int l, const Policy& pol);
+
+ template <class T>
+ inline std::vector<T> legendre_p_zeros(int l);
+
#if !BOOST_WORKAROUND(BOOST_MSVC, <= 1310)
template <class T, class Policy>
typename boost::enable_if_c<policies::is_policy<Policy>::value, typename tools::promote_args<T>::type>::type
legendre_p(int l, T x, const Policy& pol);
+ template <class T, class Policy>
+ inline typename boost::enable_if_c<policies::is_policy<Policy>::value, typename tools::promote_args<T>::type>::type
+ legendre_p_prime(int l, T x, const Policy& pol);
#endif
template <class T>
typename tools::promote_args<T>::type
@@ -1147,6 +1162,10 @@ namespace boost
\
template <class T>\
inline typename boost::math::tools::promote_args<T>::type \
+ legendre_p_prime(int l, T x){ return ::boost::math::legendre_p(l, x, Policy()); }\
+\
+ template <class T>\
+ inline typename boost::math::tools::promote_args<T>::type \
legendre_q(unsigned l, T x){ return ::boost::math::legendre_q(l, x, Policy()); }\
\
using ::boost::math::legendre_next;\
@@ -1593,5 +1612,3 @@ template <class OutputIterator, class T>\
#endif // BOOST_MATH_SPECIAL_MATH_FWD_HPP
-
-
diff --git a/boost/math/special_functions/next.hpp b/boost/math/special_functions/next.hpp
index 3921b70ce0..f23ddaffe3 100644
--- a/boost/math/special_functions/next.hpp
+++ b/boost/math/special_functions/next.hpp
@@ -27,9 +27,59 @@
namespace boost{ namespace math{
+ namespace concepts {
+
+ struct real_concept;
+ struct std_real_concept;
+
+ }
+
namespace detail{
template <class T>
+struct has_hidden_guard_digits;
+template <>
+struct has_hidden_guard_digits<float> : public mpl::false_ {};
+template <>
+struct has_hidden_guard_digits<double> : public mpl::false_ {};
+template <>
+struct has_hidden_guard_digits<long double> : public mpl::false_ {};
+#ifdef BOOST_HAS_FLOAT128
+template <>
+struct has_hidden_guard_digits<__float128> : public mpl::false_ {};
+#endif
+template <>
+struct has_hidden_guard_digits<boost::math::concepts::real_concept> : public mpl::false_ {};
+template <>
+struct has_hidden_guard_digits<boost::math::concepts::std_real_concept> : public mpl::false_ {};
+
+template <class T, bool b>
+struct has_hidden_guard_digits_10 : public mpl::false_ {};
+template <class T>
+struct has_hidden_guard_digits_10<T, true> : public mpl::bool_<(std::numeric_limits<T>::digits10 != std::numeric_limits<T>::max_digits10)> {};
+
+template <class T>
+struct has_hidden_guard_digits
+ : public has_hidden_guard_digits_10<T,
+ std::numeric_limits<T>::is_specialized
+ && (std::numeric_limits<T>::radix == 10) >
+{};
+
+template <class T>
+inline const T& normalize_value(const T& val, const mpl::false_&) { return val; }
+template <class T>
+inline T normalize_value(const T& val, const mpl::true_&)
+{
+ BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
+ BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);
+
+ boost::intmax_t shift = std::numeric_limits<T>::digits - ilogb(val) - 1;
+ T result = scalbn(val, shift);
+ result = round(result);
+ return scalbn(result, -shift);
+}
+
+template <class T>
inline T get_smallest_value(mpl::true_ const&)
{
//
@@ -93,19 +143,33 @@ struct min_shift_initializer
template <class T>
const typename min_shift_initializer<T>::init min_shift_initializer<T>::initializer;
+template <class T>
+inline T calc_min_shifted(const mpl::true_&)
+{
+ BOOST_MATH_STD_USING
+ return ldexp(tools::min_value<T>(), tools::digits<T>() + 1);
+}
+template <class T>
+inline T calc_min_shifted(const mpl::false_&)
+{
+ BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
+ BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);
+
+ return scalbn(tools::min_value<T>(), std::numeric_limits<T>::digits + 1);
+}
+
template <class T>
inline T get_min_shift_value()
{
- BOOST_MATH_STD_USING
- static const T val = ldexp(tools::min_value<T>(), tools::digits<T>() + 1);
+ static const T val = calc_min_shifted<T>(mpl::bool_<!std::numeric_limits<T>::is_specialized || std::numeric_limits<T>::radix == 2>());
min_shift_initializer<T>::force_instantiate();
return val;
}
template <class T, class Policy>
-T float_next_imp(const T& val, const Policy& pol)
+T float_next_imp(const T& val, const mpl::true_&, const Policy& pol)
{
BOOST_MATH_STD_USING
int expon;
@@ -145,6 +209,54 @@ T float_next_imp(const T& val, const Policy& pol)
diff = detail::get_smallest_value<T>();
return val + diff;
} // float_next_imp
+//
+// Special version for some base other than 2:
+//
+template <class T, class Policy>
+T float_next_imp(const T& val, const mpl::false_&, const Policy& pol)
+{
+ BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
+ BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);
+
+ BOOST_MATH_STD_USING
+ boost::intmax_t expon;
+ static const char* function = "float_next<%1%>(%1%)";
+
+ int fpclass = (boost::math::fpclassify)(val);
+
+ if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE))
+ {
+ if(val < 0)
+ return -tools::max_value<T>();
+ return policies::raise_domain_error<T>(
+ function,
+ "Argument must be finite, but got %1%", val, pol);
+ }
+
+ if(val >= tools::max_value<T>())
+ return policies::raise_overflow_error<T>(function, 0, pol);
+
+ if(val == 0)
+ return detail::get_smallest_value<T>();
+
+ if((fpclass != (int)FP_SUBNORMAL) && (fpclass != (int)FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != -tools::min_value<T>()))
+ {
+ //
+ // Special case: if the value of the least significant bit is a denorm, and the result
+ // would not be a denorm, then shift the input, increment, and shift back.
+ // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
+ //
+ return scalbn(float_next(T(scalbn(val, 2 * std::numeric_limits<T>::digits)), pol), -2 * std::numeric_limits<T>::digits);
+ }
+
+ expon = 1 + ilogb(val);
+ if(-1 == scalbn(val, -expon) * std::numeric_limits<T>::radix)
+ --expon; // reduce exponent when val is a power of base, and negative.
+ T diff = scalbn(T(1), expon - std::numeric_limits<T>::digits);
+ if(diff == 0)
+ diff = detail::get_smallest_value<T>();
+ return val + diff;
+} // float_next_imp
} // namespace detail
@@ -152,7 +264,7 @@ template <class T, class Policy>
inline typename tools::promote_args<T>::type float_next(const T& val, const Policy& pol)
{
typedef typename tools::promote_args<T>::type result_type;
- return detail::float_next_imp(static_cast<result_type>(val), pol);
+ return detail::float_next_imp(detail::normalize_value(static_cast<result_type>(val), typename detail::has_hidden_guard_digits<result_type>::type()), mpl::bool_<!std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol);
}
#if 0 //def BOOST_MSVC
@@ -187,7 +299,7 @@ inline typename tools::promote_args<T>::type float_next(const T& val)
namespace detail{
template <class T, class Policy>
-T float_prior_imp(const T& val, const Policy& pol)
+T float_prior_imp(const T& val, const mpl::true_&, const Policy& pol)
{
BOOST_MATH_STD_USING
int expon;
@@ -221,13 +333,62 @@ T float_prior_imp(const T& val, const Policy& pol)
}
T remain = frexp(val, &expon);
- if(remain == 0.5)
+ if(remain == 0.5f)
--expon; // when val is a power of two we must reduce the exponent
T diff = ldexp(T(1), expon - tools::digits<T>());
if(diff == 0)
diff = detail::get_smallest_value<T>();
return val - diff;
} // float_prior_imp
+//
+// Special version for bases other than 2:
+//
+template <class T, class Policy>
+T float_prior_imp(const T& val, const mpl::false_&, const Policy& pol)
+{
+ BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
+ BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);
+
+ BOOST_MATH_STD_USING
+ boost::intmax_t expon;
+ static const char* function = "float_prior<%1%>(%1%)";
+
+ int fpclass = (boost::math::fpclassify)(val);
+
+ if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE))
+ {
+ if(val > 0)
+ return tools::max_value<T>();
+ return policies::raise_domain_error<T>(
+ function,
+ "Argument must be finite, but got %1%", val, pol);
+ }
+
+ if(val <= -tools::max_value<T>())
+ return -policies::raise_overflow_error<T>(function, 0, pol);
+
+ if(val == 0)
+ return -detail::get_smallest_value<T>();
+
+ if((fpclass != (int)FP_SUBNORMAL) && (fpclass != (int)FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != tools::min_value<T>()))
+ {
+ //
+ // Special case: if the value of the least significant bit is a denorm, and the result
+ // would not be a denorm, then shift the input, increment, and shift back.
+ // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
+ //
+ return scalbn(float_prior(T(scalbn(val, 2 * std::numeric_limits<T>::digits)), pol), -2 * std::numeric_limits<T>::digits);
+ }
+
+ expon = 1 + ilogb(val);
+ T remain = scalbn(val, -expon);
+ if(remain * std::numeric_limits<T>::radix == 1)
+ --expon; // when val is a power of two we must reduce the exponent
+ T diff = scalbn(T(1), expon - std::numeric_limits<T>::digits);
+ if(diff == 0)
+ diff = detail::get_smallest_value<T>();
+ return val - diff;
+} // float_prior_imp
} // namespace detail
@@ -235,7 +396,7 @@ template <class T, class Policy>
inline typename tools::promote_args<T>::type float_prior(const T& val, const Policy& pol)
{
typedef typename tools::promote_args<T>::type result_type;
- return detail::float_prior_imp(static_cast<result_type>(val), pol);
+ return detail::float_prior_imp(detail::normalize_value(static_cast<result_type>(val), typename detail::has_hidden_guard_digits<result_type>::type()), mpl::bool_<!std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol);
}
#if 0 //def BOOST_MSVC
@@ -283,7 +444,7 @@ inline typename tools::promote_args<T, U>::type nextafter(const T& val, const U&
namespace detail{
template <class T, class Policy>
-T float_distance_imp(const T& a, const T& b, const Policy& pol)
+T float_distance_imp(const T& a, const T& b, const mpl::true_&, const Policy& pol)
{
BOOST_MATH_STD_USING
//
@@ -331,19 +492,23 @@ T float_distance_imp(const T& a, const T& b, const Policy& pol)
frexp(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) ? tools::min_value<T>() : a, &expon);
T upper = ldexp(T(1), expon);
T result = T(0);
- expon = tools::digits<T>() - expon;
//
// If b is greater than upper, then we *must* split the calculation
// as the size of the ULP changes with each order of magnitude change:
//
if(b > upper)
{
- result = float_distance(upper, b);
+ int expon2;
+ frexp(b, &expon2);
+ T upper2 = ldexp(T(0.5), expon2);
+ result = float_distance(upper2, b);
+ result += (expon2 - expon - 1) * ldexp(T(1), tools::digits<T>() - 1);
}
//
// Use compensated double-double addition to avoid rounding
// errors in the subtraction:
//
+ expon = tools::digits<T>() - expon;
T mb, x, y, z;
if(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) || (b - a < tools::min_value<T>()))
{
@@ -380,6 +545,113 @@ T float_distance_imp(const T& a, const T& b, const Policy& pol)
BOOST_ASSERT(result == floor(result));
return result;
} // float_distance_imp
+//
+// Special versions for bases other than 2:
+//
+template <class T, class Policy>
+T float_distance_imp(const T& a, const T& b, const mpl::false_&, const Policy& pol)
+{
+ BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
+ BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);
+
+ BOOST_MATH_STD_USING
+ //
+ // Error handling:
+ //
+ static const char* function = "float_distance<%1%>(%1%, %1%)";
+ if(!(boost::math::isfinite)(a))
+ return policies::raise_domain_error<T>(
+ function,
+ "Argument a must be finite, but got %1%", a, pol);
+ if(!(boost::math::isfinite)(b))
+ return policies::raise_domain_error<T>(
+ function,
+ "Argument b must be finite, but got %1%", b, pol);
+ //
+ // Special cases:
+ //
+ if(a > b)
+ return -float_distance(b, a, pol);
+ if(a == b)
+ return T(0);
+ if(a == 0)
+ return 1 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol));
+ if(b == 0)
+ return 1 + fabs(float_distance(static_cast<T>((a < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), a, pol));
+ if(boost::math::sign(a) != boost::math::sign(b))
+ return 2 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol))
+ + fabs(float_distance(static_cast<T>((a < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), a, pol));
+ //
+ // By the time we get here, both a and b must have the same sign, we want
+ // b > a and both postive for the following logic:
+ //
+ if(a < 0)
+ return float_distance(static_cast<T>(-b), static_cast<T>(-a), pol);
+
+ BOOST_ASSERT(a >= 0);
+ BOOST_ASSERT(b >= a);
+
+ boost::intmax_t expon;
+ //
+ // Note that if a is a denorm then the usual formula fails
+ // because we actually have fewer than tools::digits<T>()
+ // significant bits in the representation:
+ //
+ expon = 1 + ilogb(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) ? tools::min_value<T>() : a);
+ T upper = scalbn(T(1), expon);
+ T result = T(0);
+ //
+ // If b is greater than upper, then we *must* split the calculation
+ // as the size of the ULP changes with each order of magnitude change:
+ //
+ if(b > upper)
+ {
+ boost::intmax_t expon2 = 1 + ilogb(b);
+ T upper2 = scalbn(T(1), expon2 - 1);
+ result = float_distance(upper2, b);
+ result += (expon2 - expon - 1) * scalbn(T(1), std::numeric_limits<T>::digits - 1);
+ }
+ //
+ // Use compensated double-double addition to avoid rounding
+ // errors in the subtraction:
+ //
+ expon = std::numeric_limits<T>::digits - expon;
+ T mb, x, y, z;
+ if(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) || (b - a < tools::min_value<T>()))
+ {
+ //
+ // Special case - either one end of the range is a denormal, or else the difference is.
+ // The regular code will fail if we're using the SSE2 registers on Intel and either
+ // the FTZ or DAZ flags are set.
+ //
+ T a2 = scalbn(a, std::numeric_limits<T>::digits);
+ T b2 = scalbn(b, std::numeric_limits<T>::digits);
+ mb = -(std::min)(T(scalbn(upper, std::numeric_limits<T>::digits)), b2);
+ x = a2 + mb;
+ z = x - a2;
+ y = (a2 - (x - z)) + (mb - z);
+
+ expon -= std::numeric_limits<T>::digits;
+ }
+ else
+ {
+ mb = -(std::min)(upper, b);
+ x = a + mb;
+ z = x - a;
+ y = (a - (x - z)) + (mb - z);
+ }
+ if(x < 0)
+ {
+ x = -x;
+ y = -y;
+ }
+ result += scalbn(x, expon) + scalbn(y, expon);
+ //
+ // Result must be an integer:
+ //
+ BOOST_ASSERT(result == floor(result));
+ return result;
+} // float_distance_imp
} // namespace detail
@@ -387,7 +659,7 @@ template <class T, class U, class Policy>
inline typename tools::promote_args<T, U>::type float_distance(const T& a, const U& b, const Policy& pol)
{
typedef typename tools::promote_args<T, U>::type result_type;
- return detail::float_distance_imp(static_cast<result_type>(a), static_cast<result_type>(b), pol);
+ return detail::float_distance_imp(detail::normalize_value(static_cast<result_type>(a), typename detail::has_hidden_guard_digits<result_type>::type()), detail::normalize_value(static_cast<result_type>(b), typename detail::has_hidden_guard_digits<result_type>::type()), mpl::bool_<!std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol);
}
template <class T, class U>
@@ -399,7 +671,7 @@ typename tools::promote_args<T, U>::type float_distance(const T& a, const U& b)
namespace detail{
template <class T, class Policy>
-T float_advance_imp(T val, int distance, const Policy& pol)
+T float_advance_imp(T val, int distance, const mpl::true_&, const Policy& pol)
{
BOOST_MATH_STD_USING
//
@@ -478,6 +750,92 @@ T float_advance_imp(T val, int distance, const Policy& pol)
diff = distance * detail::get_smallest_value<T>();
return val += diff;
} // float_advance_imp
+//
+// Special version for bases other than 2:
+//
+template <class T, class Policy>
+T float_advance_imp(T val, int distance, const mpl::false_&, const Policy& pol)
+{
+ BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
+ BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);
+
+ BOOST_MATH_STD_USING
+ //
+ // Error handling:
+ //
+ static const char* function = "float_advance<%1%>(%1%, int)";
+
+ int fpclass = (boost::math::fpclassify)(val);
+
+ if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE))
+ return policies::raise_domain_error<T>(
+ function,
+ "Argument val must be finite, but got %1%", val, pol);
+
+ if(val < 0)
+ return -float_advance(-val, -distance, pol);
+ if(distance == 0)
+ return val;
+ if(distance == 1)
+ return float_next(val, pol);
+ if(distance == -1)
+ return float_prior(val, pol);
+
+ if(fabs(val) < detail::get_min_shift_value<T>())
+ {
+ //
+ // Special case: if the value of the least significant bit is a denorm,
+ // implement in terms of float_next/float_prior.
+ // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
+ //
+ if(distance > 0)
+ {
+ do{ val = float_next(val, pol); } while(--distance);
+ }
+ else
+ {
+ do{ val = float_prior(val, pol); } while(++distance);
+ }
+ return val;
+ }
+
+ boost::intmax_t expon = 1 + ilogb(val);
+ T limit = scalbn(T(1), distance < 0 ? expon - 1 : expon);
+ if(val <= tools::min_value<T>())
+ {
+ limit = sign(T(distance)) * tools::min_value<T>();
+ }
+ T limit_distance = float_distance(val, limit);
+ while(fabs(limit_distance) < abs(distance))
+ {
+ distance -= itrunc(limit_distance);
+ val = limit;
+ if(distance < 0)
+ {
+ limit /= std::numeric_limits<T>::radix;
+ expon--;
+ }
+ else
+ {
+ limit *= std::numeric_limits<T>::radix;
+ expon++;
+ }
+ limit_distance = float_distance(val, limit);
+ if(distance && (limit_distance == 0))
+ {
+ return policies::raise_evaluation_error<T>(function, "Internal logic failed while trying to increment floating point value %1%: most likely your FPU is in non-IEEE conforming mode.", val, pol);
+ }
+ }
+ /*expon = 1 + ilogb(val);
+ if((1 == scalbn(val, 1 + expon)) && (distance < 0))
+ --expon;*/
+ T diff = 0;
+ if(val != 0)
+ diff = distance * scalbn(T(1), expon - std::numeric_limits<T>::digits);
+ if(diff == 0)
+ diff = distance * detail::get_smallest_value<T>();
+ return val += diff;
+} // float_advance_imp
} // namespace detail
@@ -485,7 +843,7 @@ template <class T, class Policy>
inline typename tools::promote_args<T>::type float_advance(T val, int distance, const Policy& pol)
{
typedef typename tools::promote_args<T>::type result_type;
- return detail::float_advance_imp(static_cast<result_type>(val), distance, pol);
+ return detail::float_advance_imp(detail::normalize_value(static_cast<result_type>(val), typename detail::has_hidden_guard_digits<result_type>::type()), distance, mpl::bool_<!std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol);
}
template <class T>
diff --git a/boost/math/special_functions/sign.hpp b/boost/math/special_functions/sign.hpp
index 3324c90a87..5cb21bac54 100644
--- a/boost/math/special_functions/sign.hpp
+++ b/boost/math/special_functions/sign.hpp
@@ -27,7 +27,7 @@ namespace detail {
template<class T>
inline int signbit_impl(T x, native_tag const&)
{
- return (std::signbit)(x);
+ return (std::signbit)(x) ? 1 : 0;
}
#endif
diff --git a/boost/math/special_functions/ulp.hpp b/boost/math/special_functions/ulp.hpp
index 3d78a1c7c2..bb0332a7e0 100644
--- a/boost/math/special_functions/ulp.hpp
+++ b/boost/math/special_functions/ulp.hpp
@@ -18,7 +18,7 @@
namespace boost{ namespace math{ namespace detail{
template <class T, class Policy>
-T ulp_imp(const T& val, const Policy& pol)
+T ulp_imp(const T& val, const mpl::true_&, const Policy& pol)
{
BOOST_MATH_STD_USING
int expon;
@@ -48,6 +48,40 @@ T ulp_imp(const T& val, const Policy& pol)
diff = detail::get_smallest_value<T>();
return diff;
}
+// non-binary version:
+template <class T, class Policy>
+T ulp_imp(const T& val, const mpl::false_&, const Policy& pol)
+{
+ BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
+ BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);
+ BOOST_MATH_STD_USING
+ int expon;
+ static const char* function = "ulp<%1%>(%1%)";
+
+ int fpclass = (boost::math::fpclassify)(val);
+
+ if(fpclass == (int)FP_NAN)
+ {
+ return policies::raise_domain_error<T>(
+ function,
+ "Argument must be finite, but got %1%", val, pol);
+ }
+ else if((fpclass == (int)FP_INFINITE) || (fabs(val) >= tools::max_value<T>()))
+ {
+ return (val < 0 ? -1 : 1) * policies::raise_overflow_error<T>(function, 0, pol);
+ }
+ else if(fpclass == FP_ZERO)
+ return detail::get_smallest_value<T>();
+ //
+ // This code is almost the same as that for float_next, except for negative integers,
+ // where we preserve the relation ulp(x) == ulp(-x) as does Java:
+ //
+ expon = 1 + ilogb(fabs(val));
+ T diff = scalbn(T(1), expon - std::numeric_limits<T>::digits);
+ if(diff == 0)
+ diff = detail::get_smallest_value<T>();
+ return diff;
+}
}
@@ -55,7 +89,7 @@ template <class T, class Policy>
inline typename tools::promote_args<T>::type ulp(const T& val, const Policy& pol)
{
typedef typename tools::promote_args<T>::type result_type;
- return detail::ulp_imp(static_cast<result_type>(val), pol);
+ return detail::ulp_imp(static_cast<result_type>(val), mpl::bool_<!std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol);
}
template <class T>
generated by cgit v1.2.3 (git 2.25.1) at 2024-05-07 05:00:57 +0000