summaryrefslogtreecommitdiff
path: root/boost/math
diff options
context:
space:
mode:
Diffstat (limited to 'boost/math')
-rw-r--r--boost/math/common_factor_ct.hpp79
-rw-r--r--boost/math/common_factor_rt.hpp419
-rw-r--r--boost/math/concepts/std_real_concept.hpp4
-rw-r--r--boost/math/constants/constants.hpp3
-rw-r--r--boost/math/interpolators/barycentric_rational.hpp70
-rw-r--r--boost/math/interpolators/cubic_b_spline.hpp78
-rw-r--r--boost/math/interpolators/detail/barycentric_rational_detail.hpp151
-rw-r--r--boost/math/interpolators/detail/cubic_b_spline_detail.hpp287
-rw-r--r--boost/math/quadrature/trapezoidal.hpp120
-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
-rw-r--r--boost/math/tools/polynomial.hpp6
-rw-r--r--boost/math/tools/polynomial_gcd.hpp209
22 files changed, 1787 insertions, 522 deletions
diff --git a/boost/math/common_factor_ct.hpp b/boost/math/common_factor_ct.hpp
index bf58b94eb2..3ca0905945 100644
--- a/boost/math/common_factor_ct.hpp
+++ b/boost/math/common_factor_ct.hpp
@@ -1,6 +1,6 @@
// Boost common_factor_ct.hpp header file ----------------------------------//
-// (C) Copyright Daryle Walker and Stephen Cleary 2001-2002.
+// (C) Copyright John Maddock 2017.
// Distributed under 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)
@@ -10,85 +10,16 @@
#ifndef BOOST_MATH_COMMON_FACTOR_CT_HPP
#define BOOST_MATH_COMMON_FACTOR_CT_HPP
-#include <boost/math_fwd.hpp> // self include
-#include <boost/config.hpp> // for BOOST_STATIC_CONSTANT, etc.
-#include <boost/mpl/integral_c.hpp>
+#include <boost/integer/common_factor_ct.hpp>
namespace boost
{
namespace math
{
-// Implementation details --------------------------------------------------//
-
-namespace detail
-{
- // Build GCD with Euclid's recursive algorithm
- template < static_gcd_type Value1, static_gcd_type Value2 >
- struct static_gcd_helper_t
- {
- private:
- BOOST_STATIC_CONSTANT( static_gcd_type, new_value1 = Value2 );
- BOOST_STATIC_CONSTANT( static_gcd_type, new_value2 = Value1 % Value2 );
-
- #ifndef __BORLANDC__
- #define BOOST_DETAIL_GCD_HELPER_VAL(Value) static_cast<static_gcd_type>(Value)
- #else
- typedef static_gcd_helper_t self_type;
- #define BOOST_DETAIL_GCD_HELPER_VAL(Value) (self_type:: Value )
- #endif
-
- typedef static_gcd_helper_t< BOOST_DETAIL_GCD_HELPER_VAL(new_value1),
- BOOST_DETAIL_GCD_HELPER_VAL(new_value2) > next_step_type;
-
- #undef BOOST_DETAIL_GCD_HELPER_VAL
-
- public:
- BOOST_STATIC_CONSTANT( static_gcd_type, value = next_step_type::value );
- };
-
- // Non-recursive case
- template < static_gcd_type Value1 >
- struct static_gcd_helper_t< Value1, 0UL >
- {
- BOOST_STATIC_CONSTANT( static_gcd_type, value = Value1 );
- };
-
- // Build the LCM from the GCD
- template < static_gcd_type Value1, static_gcd_type Value2 >
- struct static_lcm_helper_t
- {
- typedef static_gcd_helper_t<Value1, Value2> gcd_type;
-
- BOOST_STATIC_CONSTANT( static_gcd_type, value = Value1 / gcd_type::value
- * Value2 );
- };
-
- // Special case for zero-GCD values
- template < >
- struct static_lcm_helper_t< 0UL, 0UL >
- {
- BOOST_STATIC_CONSTANT( static_gcd_type, value = 0UL );
- };
-
-} // namespace detail
-
-
-// Compile-time greatest common divisor evaluator class declaration --------//
-
-template < static_gcd_type Value1, static_gcd_type Value2 >
-struct static_gcd : public mpl::integral_c<static_gcd_type, (detail::static_gcd_helper_t<Value1, Value2>::value) >
-{
-}; // boost::math::static_gcd
-
-
-// Compile-time least common multiple evaluator class declaration ----------//
-
-template < static_gcd_type Value1, static_gcd_type Value2 >
-struct static_lcm : public mpl::integral_c<static_gcd_type, (detail::static_lcm_helper_t<Value1, Value2>::value) >
-{
-}; // boost::math::static_lcm
-
+ using boost::integer::static_gcd;
+ using boost::integer::static_lcm;
+ using boost::integer::static_gcd_type;
} // namespace math
} // namespace boost
diff --git a/boost/math/common_factor_rt.hpp b/boost/math/common_factor_rt.hpp
index acde21d2c8..42d9edfc04 100644
--- a/boost/math/common_factor_rt.hpp
+++ b/boost/math/common_factor_rt.hpp
@@ -1,4 +1,4 @@
-// (C) Copyright Jeremy William Murphy 2016.
+// (C) Copyright John Maddock 2017.
// Use, modification and distribution are subject to the
// Boost Software License, Version 1.0. (See accompanying file
@@ -7,422 +7,19 @@
#ifndef BOOST_MATH_COMMON_FACTOR_RT_HPP
#define BOOST_MATH_COMMON_FACTOR_RT_HPP
-#include <boost/assert.hpp>
-#include <boost/core/enable_if.hpp>
-#include <boost/mpl/and.hpp>
-#include <boost/type_traits.hpp>
-
-#include <boost/config.hpp> // for BOOST_NESTED_TEMPLATE, etc.
-#include <boost/limits.hpp> // for std::numeric_limits
-#include <climits> // for CHAR_MIN
-#include <boost/detail/workaround.hpp>
-#include <iterator>
-#include <algorithm>
-#include <limits>
-
-#if (defined(BOOST_MSVC) || (defined(__clang__) && defined(__c2__)) || (defined(BOOST_INTEL) && defined(_MSC_VER))) && (defined(_M_IX86) || defined(_M_X64))
-#include <intrin.h>
-#endif
-
-#ifdef BOOST_MSVC
-#pragma warning(push)
-#pragma warning(disable:4127 4244) // Conditional expression is constant
-#endif
+#include <boost/integer/common_factor_rt.hpp>
namespace boost {
namespace math {
- template <class T, bool a = is_unsigned<T>::value || (std::numeric_limits<T>::is_specialized && !std::numeric_limits<T>::is_signed)>
- struct gcd_traits_abs_defaults
- {
- inline static const T& abs(const T& val) { return val; }
- };
- template <class T>
- struct gcd_traits_abs_defaults<T, false>
- {
- inline static T abs(const T& val)
- {
- using std::abs;
- return abs(val);
- }
- };
-
- template <class T>
- struct gcd_traits_defaults : public gcd_traits_abs_defaults<T>
- {
- BOOST_FORCEINLINE static unsigned make_odd(T& val)
- {
- unsigned r = 0;
- while(!(val & 1u))
- {
- val >>= 1;
- ++r;
- }
- return r;
- }
- inline static bool less(const T& a, const T& b)
- {
- return a < b;
- }
-
- enum method_type
- {
- method_euclid = 0,
- method_binary = 1,
- method_mixed = 2,
- };
-
- static const method_type method =
- boost::has_right_shift_assign<T>::value && boost::has_left_shift_assign<T>::value && boost::has_less<T>::value && boost::has_modulus<T>::value
- ? method_mixed :
- boost::has_right_shift_assign<T>::value && boost::has_left_shift_assign<T>::value && boost::has_less<T>::value
- ? method_binary : method_euclid;
- };
- //
- // Default gcd_traits just inherits from defaults:
- //
- template <class T>
- struct gcd_traits : public gcd_traits_defaults<T> {};
- //
- // Special handling for polynomials:
- //
- namespace tools {
- template <class T>
- class polynomial;
- }
-
- template <class T>
- struct gcd_traits<boost::math::tools::polynomial<T> > : public gcd_traits_defaults<T>
- {
- static const boost::math::tools::polynomial<T>& abs(const boost::math::tools::polynomial<T>& val) { return val; }
- };
- //
- // Some platforms have fast bitscan operations, that allow us to implement
- // make_odd much more efficiently:
- //
-#if (defined(BOOST_MSVC) || (defined(__clang__) && defined(__c2__)) || (defined(BOOST_INTEL) && defined(_MSC_VER))) && (defined(_M_IX86) || defined(_M_X64))
-#pragma intrinsic(_BitScanForward,)
- template <>
- struct gcd_traits<unsigned long> : public gcd_traits_defaults<unsigned long>
- {
- BOOST_FORCEINLINE static unsigned find_lsb(unsigned long val)
- {
- unsigned long result;
- _BitScanForward(&result, val);
- return result;
- }
- BOOST_FORCEINLINE static unsigned make_odd(unsigned long& val)
- {
- unsigned result = find_lsb(val);
- val >>= result;
- return result;
- }
- };
-
-#ifdef _M_X64
-#pragma intrinsic(_BitScanForward64)
- template <>
- struct gcd_traits<unsigned __int64> : public gcd_traits_defaults<unsigned __int64>
- {
- BOOST_FORCEINLINE static unsigned find_lsb(unsigned __int64 mask)
- {
- unsigned long result;
- _BitScanForward64(&result, mask);
- return result;
- }
- BOOST_FORCEINLINE static unsigned make_odd(unsigned __int64& val)
- {
- unsigned result = find_lsb(val);
- val >>= result;
- return result;
- }
- };
-#endif
- //
- // Other integer type are trivial adaptations of the above,
- // this works for signed types too, as by the time these functions
- // are called, all values are > 0.
- //
- template <> struct gcd_traits<long> : public gcd_traits_defaults<long>
- { BOOST_FORCEINLINE static unsigned make_odd(long& val){ unsigned result = gcd_traits<unsigned long>::find_lsb(val); val >>= result; return result; } };
- template <> struct gcd_traits<unsigned int> : public gcd_traits_defaults<unsigned int>
- { BOOST_FORCEINLINE static unsigned make_odd(unsigned int& val){ unsigned result = gcd_traits<unsigned long>::find_lsb(val); val >>= result; return result; } };
- template <> struct gcd_traits<int> : public gcd_traits_defaults<int>
- { BOOST_FORCEINLINE static unsigned make_odd(int& val){ unsigned result = gcd_traits<unsigned long>::find_lsb(val); val >>= result; return result; } };
- template <> struct gcd_traits<unsigned short> : public gcd_traits_defaults<unsigned short>
- { BOOST_FORCEINLINE static unsigned make_odd(unsigned short& val){ unsigned result = gcd_traits<unsigned long>::find_lsb(val); val >>= result; return result; } };
- template <> struct gcd_traits<short> : public gcd_traits_defaults<short>
- { BOOST_FORCEINLINE static unsigned make_odd(short& val){ unsigned result = gcd_traits<unsigned long>::find_lsb(val); val >>= result; return result; } };
- template <> struct gcd_traits<unsigned char> : public gcd_traits_defaults<unsigned char>
- { BOOST_FORCEINLINE static unsigned make_odd(unsigned char& val){ unsigned result = gcd_traits<unsigned long>::find_lsb(val); val >>= result; return result; } };
- template <> struct gcd_traits<signed char> : public gcd_traits_defaults<signed char>
- { BOOST_FORCEINLINE static signed make_odd(signed char& val){ signed result = gcd_traits<unsigned long>::find_lsb(val); val >>= result; return result; } };
- template <> struct gcd_traits<char> : public gcd_traits_defaults<char>
- { BOOST_FORCEINLINE static unsigned make_odd(char& val){ unsigned result = gcd_traits<unsigned long>::find_lsb(val); val >>= result; return result; } };
- template <> struct gcd_traits<wchar_t> : public gcd_traits_defaults<wchar_t>
- { BOOST_FORCEINLINE static unsigned make_odd(wchar_t& val){ unsigned result = gcd_traits<unsigned long>::find_lsb(val); val >>= result; return result; } };
-#ifdef _M_X64
- template <> struct gcd_traits<__int64> : public gcd_traits_defaults<__int64>
- { BOOST_FORCEINLINE static unsigned make_odd(__int64& val){ unsigned result = gcd_traits<unsigned __int64>::find_lsb(val); val >>= result; return result; } };
-#endif
+ using boost::integer::gcd;
+ using boost::integer::lcm;
+ using boost::integer::gcd_range;
+ using boost::integer::lcm_range;
+ using boost::integer::gcd_evaluator;
+ using boost::integer::lcm_evaluator;
-#elif defined(BOOST_GCC) || defined(__clang__) || (defined(BOOST_INTEL) && defined(__GNUC__))
-
- template <>
- struct gcd_traits<unsigned> : public gcd_traits_defaults<unsigned>
- {
- BOOST_FORCEINLINE static unsigned find_lsb(unsigned mask)
- {
- return __builtin_ctz(mask);
- }
- BOOST_FORCEINLINE static unsigned make_odd(unsigned& val)
- {
- unsigned result = find_lsb(val);
- val >>= result;
- return result;
- }
- };
- template <>
- struct gcd_traits<unsigned long> : public gcd_traits_defaults<unsigned long>
- {
- BOOST_FORCEINLINE static unsigned find_lsb(unsigned long mask)
- {
- return __builtin_ctzl(mask);
- }
- BOOST_FORCEINLINE static unsigned make_odd(unsigned long& val)
- {
- unsigned result = find_lsb(val);
- val >>= result;
- return result;
- }
- };
- template <>
- struct gcd_traits<boost::ulong_long_type> : public gcd_traits_defaults<boost::ulong_long_type>
- {
- BOOST_FORCEINLINE static unsigned find_lsb(boost::ulong_long_type mask)
- {
- return __builtin_ctzll(mask);
- }
- BOOST_FORCEINLINE static unsigned make_odd(boost::ulong_long_type& val)
- {
- unsigned result = find_lsb(val);
- val >>= result;
- return result;
- }
- };
- //
- // Other integer type are trivial adaptations of the above,
- // this works for signed types too, as by the time these functions
- // are called, all values are > 0.
- //
- template <> struct gcd_traits<boost::long_long_type> : public gcd_traits_defaults<boost::long_long_type>
- {
- BOOST_FORCEINLINE static unsigned make_odd(boost::long_long_type& val) { unsigned result = gcd_traits<boost::ulong_long_type>::find_lsb(val); val >>= result; return result; }
- };
- template <> struct gcd_traits<long> : public gcd_traits_defaults<long>
- {
- BOOST_FORCEINLINE static unsigned make_odd(long& val) { unsigned result = gcd_traits<unsigned long>::find_lsb(val); val >>= result; return result; }
- };
- template <> struct gcd_traits<int> : public gcd_traits_defaults<int>
- {
- BOOST_FORCEINLINE static unsigned make_odd(int& val) { unsigned result = gcd_traits<unsigned long>::find_lsb(val); val >>= result; return result; }
- };
- template <> struct gcd_traits<unsigned short> : public gcd_traits_defaults<unsigned short>
- {
- BOOST_FORCEINLINE static unsigned make_odd(unsigned short& val) { unsigned result = gcd_traits<unsigned>::find_lsb(val); val >>= result; return result; }
- };
- template <> struct gcd_traits<short> : public gcd_traits_defaults<short>
- {
- BOOST_FORCEINLINE static unsigned make_odd(short& val) { unsigned result = gcd_traits<unsigned>::find_lsb(val); val >>= result; return result; }
- };
- template <> struct gcd_traits<unsigned char> : public gcd_traits_defaults<unsigned char>
- {
- BOOST_FORCEINLINE static unsigned make_odd(unsigned char& val) { unsigned result = gcd_traits<unsigned>::find_lsb(val); val >>= result; return result; }
- };
- template <> struct gcd_traits<signed char> : public gcd_traits_defaults<signed char>
- {
- BOOST_FORCEINLINE static signed make_odd(signed char& val) { signed result = gcd_traits<unsigned>::find_lsb(val); val >>= result; return result; }
- };
- template <> struct gcd_traits<char> : public gcd_traits_defaults<char>
- {
- BOOST_FORCEINLINE static unsigned make_odd(char& val) { unsigned result = gcd_traits<unsigned>::find_lsb(val); val >>= result; return result; }
- };
- template <> struct gcd_traits<wchar_t> : public gcd_traits_defaults<wchar_t>
- {
- BOOST_FORCEINLINE static unsigned make_odd(wchar_t& val) { unsigned result = gcd_traits<unsigned>::find_lsb(val); val >>= result; return result; }
- };
-#endif
-
-namespace detail
-{
-
- //
- // The Mixed Binary Euclid Algorithm
- // Sidi Mohamed Sedjelmaci
- // Electronic Notes in Discrete Mathematics 35 (2009) 169-176
- //
- template <class T>
- T mixed_binary_gcd(T u, T v)
- {
- using std::swap;
- if(gcd_traits<T>::less(u, v))
- swap(u, v);
-
- unsigned shifts = 0;
-
- if(!u)
- return v;
- if(!v)
- return u;
-
- shifts = (std::min)(gcd_traits<T>::make_odd(u), gcd_traits<T>::make_odd(v));
-
- while(gcd_traits<T>::less(1, v))
- {
- u %= v;
- v -= u;
- if(!u)
- return v << shifts;
- if(!v)
- return u << shifts;
- gcd_traits<T>::make_odd(u);
- gcd_traits<T>::make_odd(v);
- if(gcd_traits<T>::less(u, v))
- swap(u, v);
- }
- return (v == 1 ? v : u) << shifts;
}
-
- /** Stein gcd (aka 'binary gcd')
- *
- * From Mathematics to Generic Programming, Alexander Stepanov, Daniel Rose
- */
- template <typename SteinDomain>
- SteinDomain Stein_gcd(SteinDomain m, SteinDomain n)
- {
- using std::swap;
- BOOST_ASSERT(m >= 0);
- BOOST_ASSERT(n >= 0);
- if (m == SteinDomain(0))
- return n;
- if (n == SteinDomain(0))
- return m;
- // m > 0 && n > 0
- int d_m = gcd_traits<SteinDomain>::make_odd(m);
- int d_n = gcd_traits<SteinDomain>::make_odd(n);
- // odd(m) && odd(n)
- while (m != n)
- {
- if (n > m)
- swap(n, m);
- m -= n;
- gcd_traits<SteinDomain>::make_odd(m);
- }
- // m == n
- m <<= (std::min)(d_m, d_n);
- return m;
- }
-
-
- /** Euclidean algorithm
- *
- * From Mathematics to Generic Programming, Alexander Stepanov, Daniel Rose
- *
- */
- template <typename EuclideanDomain>
- inline EuclideanDomain Euclid_gcd(EuclideanDomain a, EuclideanDomain b)
- {
- using std::swap;
- while (b != EuclideanDomain(0))
- {
- a %= b;
- swap(a, b);
- }
- return a;
- }
-
-
- template <typename T>
- inline BOOST_DEDUCED_TYPENAME enable_if_c<gcd_traits<T>::method == gcd_traits<T>::method_mixed, T>::type
- optimal_gcd_select(T const &a, T const &b)
- {
- return detail::mixed_binary_gcd(a, b);
- }
-
- template <typename T>
- inline BOOST_DEDUCED_TYPENAME enable_if_c<gcd_traits<T>::method == gcd_traits<T>::method_binary, T>::type
- optimal_gcd_select(T const &a, T const &b)
- {
- return detail::Stein_gcd(a, b);
- }
-
- template <typename T>
- inline BOOST_DEDUCED_TYPENAME enable_if_c<gcd_traits<T>::method == gcd_traits<T>::method_euclid, T>::type
- optimal_gcd_select(T const &a, T const &b)
- {
- return detail::Euclid_gcd(a, b);
- }
-
- template <class T>
- inline T lcm_imp(const T& a, const T& b)
- {
- T temp = boost::math::detail::optimal_gcd_select(a, b);
-#if BOOST_WORKAROUND(BOOST_GCC_VERSION, < 40500)
- return (temp != T(0)) ? T(a / temp * b) : T(0);
-#else
- return temp ? T(a / temp * b) : T(0);
-#endif
- }
-
-} // namespace detail
-
-
-template <typename Integer>
-inline Integer gcd(Integer const &a, Integer const &b)
-{
- return detail::optimal_gcd_select(static_cast<Integer>(gcd_traits<Integer>::abs(a)), static_cast<Integer>(gcd_traits<Integer>::abs(b)));
-}
-
-template <typename Integer>
-inline Integer lcm(Integer const &a, Integer const &b)
-{
- return detail::lcm_imp(static_cast<Integer>(gcd_traits<Integer>::abs(a)), static_cast<Integer>(gcd_traits<Integer>::abs(b)));
}
-/**
- * Knuth, The Art of Computer Programming: Volume 2, Third edition, 1998
- * Chapter 4.5.2, Algorithm C: Greatest common divisor of n integers.
- *
- * Knuth counts down from n to zero but we naturally go from first to last.
- * We also return the termination position because it might be useful to know.
- *
- * Partly by quirk, partly by design, this algorithm is defined for n = 1,
- * because the gcd of {x} is x. It is not defined for n = 0.
- *
- * @tparam I Input iterator.
- * @return The gcd of the range and the iterator position at termination.
- */
-template <typename I>
-std::pair<typename std::iterator_traits<I>::value_type, I>
-gcd_range(I first, I last)
-{
- BOOST_ASSERT(first != last);
- typedef typename std::iterator_traits<I>::value_type T;
-
- T d = *first++;
- while (d != T(1) && first != last)
- {
- d = gcd(d, *first);
- first++;
- }
- return std::make_pair(d, first);
-}
-
-} // namespace math
-} // namespace boost
-
-#ifdef BOOST_MSVC
-#pragma warning(pop)
-#endif
-
#endif // BOOST_MATH_COMMON_FACTOR_RT_HPP
diff --git a/boost/math/concepts/std_real_concept.hpp b/boost/math/concepts/std_real_concept.hpp
index b297501d98..b8657bc6a3 100644
--- a/boost/math/concepts/std_real_concept.hpp
+++ b/boost/math/concepts/std_real_concept.hpp
@@ -399,7 +399,3 @@ using concepts::llround;
} // namespace boost
#endif // BOOST_MATH_STD_REAL_CONCEPT_HPP
-
-
-
-
diff --git a/boost/math/constants/constants.hpp b/boost/math/constants/constants.hpp
index b5db9cedf7..8c5c4105d4 100644
--- a/boost/math/constants/constants.hpp
+++ b/boost/math/constants/constants.hpp
@@ -266,6 +266,7 @@ namespace boost{ namespace math
BOOST_DEFINE_MATH_CONSTANT(third, 3.333333333333333333333333333333333333e-01, "3.33333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333e-01")
BOOST_DEFINE_MATH_CONSTANT(twothirds, 6.666666666666666666666666666666666666e-01, "6.66666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666667e-01")
BOOST_DEFINE_MATH_CONSTANT(two_thirds, 6.666666666666666666666666666666666666e-01, "6.66666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666667e-01")
+ BOOST_DEFINE_MATH_CONSTANT(sixth, 1.66666666666666666666666666666666666666666e-01, "1.66666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666667e-01")
BOOST_DEFINE_MATH_CONSTANT(three_quarters, 7.500000000000000000000000000000000000e-01, "7.50000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e-01")
BOOST_DEFINE_MATH_CONSTANT(root_two, 1.414213562373095048801688724209698078e+00, "1.41421356237309504880168872420969807856967187537694807317667973799073247846210703885038753432764157273501384623e+00")
BOOST_DEFINE_MATH_CONSTANT(root_three, 1.732050807568877293527446341505872366e+00, "1.73205080756887729352744634150587236694280525381038062805580697945193301690880003708114618675724857567562614142e+00")
@@ -343,5 +344,3 @@ namespace boost{ namespace math
#include <boost/math/constants/calculate_constants.hpp>
#endif // BOOST_MATH_CONSTANTS_CONSTANTS_INCLUDED
-
-
diff --git a/boost/math/interpolators/barycentric_rational.hpp b/boost/math/interpolators/barycentric_rational.hpp
new file mode 100644
index 0000000000..79bab9042d
--- /dev/null
+++ b/boost/math/interpolators/barycentric_rational.hpp
@@ -0,0 +1,70 @@
+/*
+ * 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)
+ *
+ * Given N samples (t_i, y_i) which are irregularly spaced, this routine constructs an
+ * interpolant s which is constructed in O(N) time, occupies O(N) space, and can be evaluated in O(N) time.
+ * The interpolation is stable, unless one point is incredibly close to another, and the next point is incredibly far.
+ * The measure of this stability is the "local mesh ratio", which can be queried from the routine.
+ * Pictorially, the following t_i spacing is bad (has a high local mesh ratio)
+ * || | | | |
+ * and this t_i spacing is good (has a low local mesh ratio)
+ * | | | | | | | | | |
+ *
+ *
+ * If f is C^{d+2}, then the interpolant is O(h^(d+1)) accurate, where d is the interpolation order.
+ * A disadvantage of this interpolant is that it does not reproduce rational functions; for example, 1/(1+x^2) is not interpolated exactly.
+ *
+ * References:
+ * Floater, Michael S., and Kai Hormann. "Barycentric rational interpolation with no poles and high rates of approximation." Numerische Mathematik 107.2 (2007): 315-331.
+ * Press, William H., et al. "Numerical recipes third edition: the art of scientific computing." Cambridge University Press 32 (2007): 10013-2473.
+ */
+
+#ifndef BOOST_MATH_INTERPOLATORS_BARYCENTRIC_RATIONAL_HPP
+#define BOOST_MATH_INTERPOLATORS_BARYCENTRIC_RATIONAL_HPP
+
+#include <memory>
+#include <boost/math/interpolators/detail/barycentric_rational_detail.hpp>
+
+namespace boost{ namespace math{
+
+template<class Real>
+class barycentric_rational
+{
+public:
+ barycentric_rational(const Real* const x, const Real* const y, size_t n, size_t approximation_order = 3);
+
+ template <class InputIterator1, class InputIterator2>
+ barycentric_rational(InputIterator1 start_x, InputIterator1 end_x, InputIterator2 start_y, size_t approximation_order = 3, typename boost::disable_if_c<boost::is_integral<InputIterator2>::value>::type* = 0);
+
+ Real operator()(Real x) const;
+
+private:
+ std::shared_ptr<detail::barycentric_rational_imp<Real>> m_imp;
+};
+
+template <class Real>
+barycentric_rational<Real>::barycentric_rational(const Real* const x, const Real* const y, size_t n, size_t approximation_order):
+ m_imp(std::make_shared<detail::barycentric_rational_imp<Real>>(x, x + n, y, approximation_order))
+{
+ return;
+}
+
+template <class Real>
+template <class InputIterator1, class InputIterator2>
+barycentric_rational<Real>::barycentric_rational(InputIterator1 start_x, InputIterator1 end_x, InputIterator2 start_y, size_t approximation_order, typename boost::disable_if_c<boost::is_integral<InputIterator2>::value>::type*)
+ : m_imp(std::make_shared<detail::barycentric_rational_imp<Real>>(start_x, end_x, start_y, approximation_order))
+{
+}
+
+template<class Real>
+Real barycentric_rational<Real>::operator()(Real x) const
+{
+ return m_imp->operator()(x);
+}
+
+
+}}
+#endif
diff --git a/boost/math/interpolators/cubic_b_spline.hpp b/boost/math/interpolators/cubic_b_spline.hpp
new file mode 100644
index 0000000000..73ac1d0137
--- /dev/null
+++ b/boost/math/interpolators/cubic_b_spline.hpp
@@ -0,0 +1,78 @@
+// 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)
+
+// This implements the compactly supported cubic b spline algorithm described in
+// Kress, Rainer. "Numerical analysis, volume 181 of Graduate Texts in Mathematics." (1998).
+// Splines of compact support are faster to evaluate and are better conditioned than classical cubic splines.
+
+// Let f be the function we are trying to interpolate, and s be the interpolating spline.
+// The routine constructs the interpolant in O(N) time, and evaluating s at a point takes constant time.
+// The order of accuracy depends on the regularity of the f, however, assuming f is
+// four-times continuously differentiable, the error is of O(h^4).
+// In addition, we can differentiate the spline and obtain a good interpolant for f'.
+// The main restriction of this method is that the samples of f must be evenly spaced.
+// Look for barycentric rational interpolation for non-evenly sampled data.
+// Properties:
+// - s(x_j) = f(x_j)
+// - All cubic polynomials interpolated exactly
+
+#ifndef BOOST_MATH_INTERPOLATORS_CUBIC_B_SPLINE_HPP
+#define BOOST_MATH_INTERPOLATORS_CUBIC_B_SPLINE_HPP
+
+#include <boost/math/interpolators/detail/cubic_b_spline_detail.hpp>
+
+namespace boost{ namespace math{
+
+template <class Real>
+class cubic_b_spline
+{
+public:
+ // If you don't know the value of the derivative at the endpoints, leave them as nans and the routine will estimate them.
+ // f[0] = f(a), f[length -1] = b, step_size = (b - a)/(length -1).
+ template <class BidiIterator>
+ cubic_b_spline(const BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size,
+ Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
+ Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN());
+ cubic_b_spline(const Real* const f, size_t length, Real left_endpoint, Real step_size,
+ Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
+ Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN());
+
+ cubic_b_spline() = default;
+ Real operator()(Real x) const;
+
+ Real prime(Real x) const;
+
+private:
+ std::shared_ptr<detail::cubic_b_spline_imp<Real>> m_imp;
+};
+
+template<class Real>
+cubic_b_spline<Real>::cubic_b_spline(const Real* const f, size_t length, Real left_endpoint, Real step_size,
+ Real left_endpoint_derivative, Real right_endpoint_derivative) : m_imp(std::make_shared<detail::cubic_b_spline_imp<Real>>(f, f + length, left_endpoint, step_size, left_endpoint_derivative, right_endpoint_derivative))
+{
+}
+
+template <class Real>
+template <class BidiIterator>
+cubic_b_spline<Real>::cubic_b_spline(BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size,
+ Real left_endpoint_derivative, Real right_endpoint_derivative) : m_imp(std::make_shared<detail::cubic_b_spline_imp<Real>>(f, end_p, left_endpoint, step_size, left_endpoint_derivative, right_endpoint_derivative))
+{
+}
+
+template<class Real>
+Real cubic_b_spline<Real>::operator()(Real x) const
+{
+ return m_imp->operator()(x);
+}
+
+template<class Real>
+Real cubic_b_spline<Real>::prime(Real x) const
+{
+ return m_imp->prime(x);
+}
+
+}}
+#endif
diff --git a/boost/math/interpolators/detail/barycentric_rational_detail.hpp b/boost/math/interpolators/detail/barycentric_rational_detail.hpp
new file mode 100644
index 0000000000..b853901d89
--- /dev/null
+++ b/boost/math/interpolators/detail/barycentric_rational_detail.hpp
@@ -0,0 +1,151 @@
+/*
+ * 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_INTERPOLATORS_BARYCENTRIC_RATIONAL_DETAIL_HPP
+#define BOOST_MATH_INTERPOLATORS_BARYCENTRIC_RATIONAL_DETAIL_HPP
+
+#include <vector>
+#include <boost/lexical_cast.hpp>
+#include <boost/math/special_functions/fpclassify.hpp>
+#include <boost/core/demangle.hpp>
+
+namespace boost{ namespace math{ namespace detail{
+
+template<class Real>
+class barycentric_rational_imp
+{
+public:
+ template <class InputIterator1, class InputIterator2>
+ barycentric_rational_imp(InputIterator1 start_x, InputIterator1 end_x, InputIterator2 start_y, size_t approximation_order = 3);
+
+ Real operator()(Real x) const;
+
+ // The barycentric weights are not really that interesting; except to the unit tests!
+ Real weight(size_t i) const { return m_w[i]; }
+
+private:
+ // Technically, we do not need to copy over y to m_y, or x to m_x.
+ // We could simply store a pointer. However, in doing so,
+ // we'd need to make sure the user managed the lifetime of m_y,
+ // and didn't alter its data. Since we are unlikely to run out of
+ // memory for a linearly scaling algorithm, it seems best just to make a copy.
+ std::vector<Real> m_y;
+ std::vector<Real> m_x;
+ std::vector<Real> m_w;
+};
+
+template <class Real>
+template <class InputIterator1, class InputIterator2>
+barycentric_rational_imp<Real>::barycentric_rational_imp(InputIterator1 start_x, InputIterator1 end_x, InputIterator2 start_y, size_t approximation_order)
+{
+ using std::abs;
+
+ std::ptrdiff_t n = std::distance(start_x, end_x);
+
+ if (approximation_order >= (std::size_t)n)
+ {
+ throw std::domain_error("Approximation order must be < data length.");
+ }
+
+ // Big sad memcpy to make sure the object is easy to use.
+ m_x.resize(n);
+ m_y.resize(n);
+ for(unsigned i = 0; start_x != end_x; ++start_x, ++start_y, ++i)
+ {
+ // But if we're going to do a memcpy, we can do some error checking which is inexpensive relative to the copy:
+ if(boost::math::isnan(*start_x))
+ {
+ std::string msg = std::string("x[") + boost::lexical_cast<std::string>(i) + "] is a NAN";
+ throw std::domain_error(msg);
+ }
+
+ if(boost::math::isnan(*start_y))
+ {
+ std::string msg = std::string("y[") + boost::lexical_cast<std::string>(i) + "] is a NAN";
+ throw std::domain_error(msg);
+ }
+
+ m_x[i] = *start_x;
+ m_y[i] = *start_y;
+ }
+
+ m_w.resize(n, 0);
+ for(int64_t k = 0; k < n; ++k)
+ {
+ int64_t i_min = (std::max)(k - (int64_t) approximation_order, (int64_t) 0);
+ int64_t i_max = k;
+ if (k >= n - (std::ptrdiff_t)approximation_order)
+ {
+ i_max = n - approximation_order - 1;
+ }
+
+ for(int64_t i = i_min; i <= i_max; ++i)
+ {
+ Real inv_product = 1;
+ int64_t j_max = (std::min)(static_cast<int64_t>(i + approximation_order), static_cast<int64_t>(n - 1));
+ for(int64_t j = i; j <= j_max; ++j)
+ {
+ if (j == k)
+ {
+ continue;
+ }
+
+ Real diff = m_x[k] - m_x[j];
+ if (abs(diff) < std::numeric_limits<Real>::epsilon())
+ {
+ std::string msg = std::string("Spacing between x[")
+ + boost::lexical_cast<std::string>(k) + std::string("] and x[")
+ + boost::lexical_cast<std::string>(i) + std::string("] is ")
+ + boost::lexical_cast<std::string>(diff) + std::string(", which is smaller than the epsilon of ")
+ + boost::core::demangle(typeid(Real).name());
+ throw std::logic_error(msg);
+ }
+ inv_product *= diff;
+ }
+ if (i % 2 == 0)
+ {
+ m_w[k] += 1/inv_product;
+ }
+ else
+ {
+ m_w[k] -= 1/inv_product;
+ }
+ }
+ }
+}
+
+template<class Real>
+Real barycentric_rational_imp<Real>::operator()(Real x) const
+{
+ Real numerator = 0;
+ Real denominator = 0;
+ for(size_t i = 0; i < m_x.size(); ++i)
+ {
+ // Presumably we should see if the accuracy is improved by using ULP distance of say, 5 here, instead of testing for floating point equality.
+ // However, it has been shown that if x approx x_i, but x != x_i, then inaccuracy in the numerator cancels the inaccuracy in the denominator,
+ // and the result is fairly accurate. See: http://epubs.siam.org/doi/pdf/10.1137/S0036144502417715
+ if (x == m_x[i])
+ {
+ return m_y[i];
+ }
+ Real t = m_w[i]/(x - m_x[i]);
+ numerator += t*m_y[i];
+ denominator += t;
+ }
+ return numerator/denominator;
+}
+
+/*
+ * A formula for computing the derivative of the barycentric representation is given in
+ * "Some New Aspects of Rational Interpolation", by Claus Schneider and Wilhelm Werner,
+ * Mathematics of Computation, v47, number 175, 1986.
+ * http://www.ams.org/journals/mcom/1986-47-175/S0025-5718-1986-0842136-8/S0025-5718-1986-0842136-8.pdf
+ * However, this requires a lot of machinery which is not built into the library at present.
+ * So we wait until there is a requirement to interpolate the derivative.
+ */
+}}}
+#endif
diff --git a/boost/math/interpolators/detail/cubic_b_spline_detail.hpp b/boost/math/interpolators/detail/cubic_b_spline_detail.hpp
new file mode 100644
index 0000000000..f7b2d6cd29
--- /dev/null
+++ b/boost/math/interpolators/detail/cubic_b_spline_detail.hpp
@@ -0,0 +1,287 @@
+// 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 CUBIC_B_SPLINE_DETAIL_HPP
+#define CUBIC_B_SPLINE_DETAIL_HPP
+
+#include <limits>
+#include <cmath>
+#include <vector>
+#include <memory>
+#include <boost/math/constants/constants.hpp>
+#include <boost/math/special_functions/fpclassify.hpp>
+
+namespace boost{ namespace math{ namespace detail{
+
+
+template <class Real>
+class cubic_b_spline_imp
+{
+public:
+ // If you don't know the value of the derivative at the endpoints, leave them as nans and the routine will estimate them.
+ // f[0] = f(a), f[length -1] = b, step_size = (b - a)/(length -1).
+ template <class BidiIterator>
+ cubic_b_spline_imp(BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size,
+ Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
+ Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN());
+
+ Real operator()(Real x) const;
+
+ Real prime(Real x) const;
+
+private:
+ std::vector<Real> m_beta;
+ Real m_h_inv;
+ Real m_a;
+ Real m_avg;
+};
+
+
+
+template <class Real>
+Real b3_spline(Real x)
+{
+ using std::abs;
+ Real absx = abs(x);
+ if (absx < 1)
+ {
+ Real y = 2 - absx;
+ Real z = 1 - absx;
+ return boost::math::constants::sixth<Real>()*(y*y*y - 4*z*z*z);
+ }
+ if (absx < 2)
+ {
+ Real y = 2 - absx;
+ return boost::math::constants::sixth<Real>()*y*y*y;
+ }
+ return (Real) 0;
+}
+
+template<class Real>
+Real b3_spline_prime(Real x)
+{
+ if (x < 0)
+ {
+ return -b3_spline_prime(-x);
+ }
+
+ if (x < 1)
+ {
+ return x*(3*boost::math::constants::half<Real>()*x - 2);
+ }
+ if (x < 2)
+ {
+ return -boost::math::constants::half<Real>()*(2 - x)*(2 - x);
+ }
+ return (Real) 0;
+}
+
+
+template <class Real>
+template <class BidiIterator>
+cubic_b_spline_imp<Real>::cubic_b_spline_imp(BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size,
+ Real left_endpoint_derivative, Real right_endpoint_derivative) : m_a(left_endpoint), m_avg(0)
+{
+ using boost::math::constants::third;
+
+ std::size_t length = end_p - f;
+
+ if (length < 5)
+ {
+ if (boost::math::isnan(left_endpoint_derivative) || boost::math::isnan(right_endpoint_derivative))
+ {
+ throw std::logic_error("Interpolation using a cubic b spline with derivatives estimated at the endpoints requires at least 5 points.\n");
+ }
+ if (length < 3)
+ {
+ throw std::logic_error("Interpolation using a cubic b spline requires at least 3 points.\n");
+ }
+ }
+
+ if (boost::math::isnan(left_endpoint))
+ {
+ throw std::logic_error("Left endpoint is NAN; this is disallowed.\n");
+ }
+ if (left_endpoint + length*step_size >= (std::numeric_limits<Real>::max)())
+ {
+ throw std::logic_error("Right endpoint overflows the maximum representable number of the specified precision.\n");
+ }
+ if (step_size <= 0)
+ {
+ throw std::logic_error("The step size must be strictly > 0.\n");
+ }
+ // Storing the inverse of the stepsize does provide a measurable speedup.
+ // It's not huge, but nonetheless worthwhile.
+ m_h_inv = 1/step_size;
+
+ // Following Kress's notation, s'(a) = a1, s'(b) = b1
+ Real a1 = left_endpoint_derivative;
+ // See the finite-difference table on Wikipedia for reference on how
+ // to construct high-order estimates for one-sided derivatives:
+ // https://en.wikipedia.org/wiki/Finite_difference_coefficient#Forward_and_backward_finite_difference
+ // Here, we estimate then to O(h^4), as that is the maximum accuracy we could obtain from this method.
+ if (boost::math::isnan(a1))
+ {
+ // For simple functions (linear, quadratic, so on)
+ // almost all the error comes from derivative estimation.
+ // This does pairwise summation which gives us another digit of accuracy over naive summation.
+ Real t0 = 4*(f[1] + third<Real>()*f[3]);
+ Real t1 = -(25*third<Real>()*f[0] + f[4])/4 - 3*f[2];
+ a1 = m_h_inv*(t0 + t1);
+ }
+
+ Real b1 = right_endpoint_derivative;
+ if (boost::math::isnan(b1))
+ {
+ size_t n = length - 1;
+ Real t0 = 4*(f[n-3] + third<Real>()*f[n - 1]);
+ Real t1 = -(25*third<Real>()*f[n - 4] + f[n])/4 - 3*f[n - 2];
+
+ b1 = m_h_inv*(t0 + t1);
+ }
+
+ // s(x) = \sum \alpha_i B_{3}( (x- x_i - a)/h )
+ // Of course we must reindex from Kress's notation, since he uses negative indices which make C++ unhappy.
+ m_beta.resize(length + 2, std::numeric_limits<Real>::quiet_NaN());
+
+ // Since the splines have compact support, they decay to zero very fast outside the endpoints.
+ // This is often very annoying; we'd like to evaluate the interpolant a little bit outside the
+ // boundary [a,b] without massive error.
+ // A simple way to deal with this is just to subtract the DC component off the signal, so we need the average.
+ // This algorithm for computing the average is recommended in
+ // http://www.heikohoffmann.de/htmlthesis/node134.html
+ Real t = 1;
+ for (size_t i = 0; i < length; ++i)
+ {
+ if (boost::math::isnan(f[i]))
+ {
+ std::string err = "This function you are trying to interpolate is a nan at index " + std::to_string(i) + "\n";
+ throw std::logic_error(err);
+ }
+ m_avg += (f[i] - m_avg) / t;
+ t += 1;
+ }
+
+
+ // Now we must solve an almost-tridiagonal system, which requires O(N) operations.
+ // There are, in fact 5 diagonals, but they only differ from zero on the first and last row,
+ // so we can patch up the tridiagonal row reduction algorithm to deal with two special rows.
+ // See Kress, equations 8.41
+ // The the "tridiagonal" matrix is:
+ // 1 0 -1
+ // 1 4 1
+ // 1 4 1
+ // 1 4 1
+ // ....
+ // 1 4 1
+ // 1 0 -1
+ // Numerical estimate indicate that as N->Infinity, cond(A) -> 6.9, so this matrix is good.
+ std::vector<Real> rhs(length + 2, std::numeric_limits<Real>::quiet_NaN());
+ std::vector<Real> super_diagonal(length + 2, std::numeric_limits<Real>::quiet_NaN());
+
+ rhs[0] = -2*step_size*a1;
+ rhs[rhs.size() - 1] = -2*step_size*b1;
+
+ super_diagonal[0] = 0;
+
+ for(size_t i = 1; i < rhs.size() - 1; ++i)
+ {
+ rhs[i] = 6*(f[i - 1] - m_avg);
+ super_diagonal[i] = 1;
+ }
+
+
+ // One step of row reduction on the first row to patch up the 5-diagonal problem:
+ // 1 0 -1 | r0
+ // 1 4 1 | r1
+ // mapsto:
+ // 1 0 -1 | r0
+ // 0 4 2 | r1 - r0
+ // mapsto
+ // 1 0 -1 | r0
+ // 0 1 1/2| (r1 - r0)/4
+ super_diagonal[1] = 0.5;
+ rhs[1] = (rhs[1] - rhs[0])/4;
+
+ // Now do a tridiagonal row reduction the standard way, until just before the last row:
+ for (size_t i = 2; i < rhs.size() - 1; ++i)
+ {
+ Real diagonal = 4 - super_diagonal[i - 1];
+ rhs[i] = (rhs[i] - rhs[i - 1])/diagonal;
+ super_diagonal[i] /= diagonal;
+ }
+
+ // Now the last row, which is in the form
+ // 1 sd[n-3] 0 | rhs[n-3]
+ // 0 1 sd[n-2] | rhs[n-2]
+ // 1 0 -1 | rhs[n-1]
+ Real final_subdiag = -super_diagonal[rhs.size() - 3];
+ rhs[rhs.size() - 1] = (rhs[rhs.size() - 1] - rhs[rhs.size() - 3])/final_subdiag;
+ Real final_diag = -1/final_subdiag;
+ // Now we're here:
+ // 1 sd[n-3] 0 | rhs[n-3]
+ // 0 1 sd[n-2] | rhs[n-2]
+ // 0 1 final_diag | (rhs[n-1] - rhs[n-3])/diag
+
+ final_diag = final_diag - super_diagonal[rhs.size() - 2];
+ rhs[rhs.size() - 1] = rhs[rhs.size() - 1] - rhs[rhs.size() - 2];
+
+
+ // Back substitutions:
+ m_beta[rhs.size() - 1] = rhs[rhs.size() - 1]/final_diag;
+ for(size_t i = rhs.size() - 2; i > 0; --i)
+ {
+ m_beta[i] = rhs[i] - super_diagonal[i]*m_beta[i + 1];
+ }
+ m_beta[0] = m_beta[2] + rhs[0];
+}
+
+template<class Real>
+Real cubic_b_spline_imp<Real>::operator()(Real x) const
+{
+ // See Kress, 8.40: Since B3 has compact support, we don't have to sum over all terms,
+ // just the (at most 5) whose support overlaps the argument.
+ Real z = m_avg;
+ Real t = m_h_inv*(x - m_a) + 1;
+
+ using std::max;
+ using std::min;
+ using std::ceil;
+ using std::floor;
+
+ size_t k_min = (size_t) max(static_cast<long>(0), boost::math::ltrunc(ceil(t - 2)));
+ size_t k_max = (size_t) max(min(static_cast<long>(m_beta.size() - 1), boost::math::ltrunc(floor(t + 2))), (long) 0);
+ for (size_t k = k_min; k <= k_max; ++k)
+ {
+ z += m_beta[k]*b3_spline(t - k);
+ }
+
+ return z;
+}
+
+template<class Real>
+Real cubic_b_spline_imp<Real>::prime(Real x) const
+{
+ Real z = 0;
+ Real t = m_h_inv*(x - m_a) + 1;
+
+ using std::max;
+ using std::min;
+ using std::ceil;
+ using std::floor;
+
+ size_t k_min = (size_t) max(static_cast<long>(0), boost::math::ltrunc(ceil(t - 2)));
+ size_t k_max = (size_t) min(static_cast<long>(m_beta.size() - 1), boost::math::ltrunc(floor(t + 2)));
+
+ for (size_t k = k_min; k <= k_max; ++k)
+ {
+ z += m_beta[k]*b3_spline_prime(t - k);
+ }
+ return z*m_h_inv;
+}
+
+}}}
+#endif
diff --git a/boost/math/quadrature/trapezoidal.hpp b/boost/math/quadrature/trapezoidal.hpp
new file mode 100644
index 0000000000..af466a0916
--- /dev/null
+++ b/boost/math/quadrature/trapezoidal.hpp
@@ -0,0 +1,120 @@
+/*
+ * 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)
+ *
+ * Use the adaptive trapezoidal rule to estimate the integral of periodic functions over a period,
+ * or to integrate a function whose derivative vanishes at the endpoints.
+ *
+ * If your function does not satisfy these conditions, and instead is simply continuous and bounded
+ * over the whole interval, then this routine will still converge, albeit slowly. However, there
+ * are much more efficient methods in this case, including Romberg, Simpson, and double exponential quadrature.
+ */
+
+#ifndef BOOST_MATH_QUADRATURE_TRAPEZOIDAL_HPP
+#define BOOST_MATH_QUADRATURE_TRAPEZOIDAL_HPP
+
+#include <cmath>
+#include <limits>
+#include <stdexcept>
+#include <boost/math/constants/constants.hpp>
+#include <boost/math/special_functions/fpclassify.hpp>
+#include <boost/math/policies/error_handling.hpp>
+
+namespace boost{ namespace math{ namespace quadrature {
+
+template<class F, class Real, class Policy>
+Real trapezoidal(F f, Real a, Real b, Real tol, std::size_t max_refinements, Real* error_estimate, Real* L1, const Policy& pol)
+{
+ static const char* function = "boost::math::quadrature::trapezoidal<%1%>(F, %1%, %1%, %1%)";
+ using std::abs;
+ using boost::math::constants::half;
+ if(a >= b)
+ {
+ return boost::math::policies::raise_domain_error(function, "a < b for integration over the region [a, b] is required, but got a = %1%.\n", a, pol);
+ }
+ if (!(boost::math::isfinite)(a))
+ {
+ return boost::math::policies::raise_domain_error(function, "Left endpoint of integration must be finite for adaptive trapezoidal integration but got a = %1%.\n", a, pol);
+ }
+ if (!(boost::math::isfinite)(b))
+ {
+ return boost::math::policies::raise_domain_error(function, "Right endpoint of integration must be finite for adaptive trapedzoidal integration but got b = %1%.\n", b, pol);
+ }
+
+ Real ya = f(a);
+ Real yb = f(b);
+ Real h = (b - a)*half<Real>();
+ Real I0 = (ya + yb)*h;
+ Real IL0 = (abs(ya) + abs(yb))*h;
+
+ Real yh = f(a + h);
+ Real I1 = half<Real>()*I0 + yh*h;
+ Real IL1 = half<Real>()*IL0 + abs(yh)*h;
+
+ // The recursion is:
+ // I_k = 1/2 I_{k-1} + 1/2^k \sum_{j=1; j odd, j < 2^k} f(a + j(b-a)/2^k)
+ std::size_t k = 2;
+ // We want to go through at least 4 levels so we have sampled the function at least 10 times.
+ // Otherwise, we could terminate prematurely and miss essential features.
+ // This is of course possible anyway, but 10 samples seems to be a reasonable compromise.
+ Real error = abs(I0 - I1);
+ while (k < 4 || (k < max_refinements && error > tol*IL1) )
+ {
+ I0 = I1;
+ IL0 = IL1;
+
+ I1 = half<Real>()*I0;
+ IL1 = half<Real>()*IL0;
+ std::size_t p = static_cast<std::size_t>(1u) << k;
+ h *= half<Real>();
+ Real sum = 0;
+ Real absum = 0;
+
+ for(std::size_t j = 1; j < p; j += 2)
+ {
+ Real y = f(a + j*h);
+ sum += y;
+ absum += abs(y);
+ }
+
+ I1 += sum*h;
+ IL1 += absum*h;
+ ++k;
+ error = abs(I0 - I1);
+ }
+
+ if (error_estimate)
+ {
+ *error_estimate = error;
+ }
+
+ if (L1)
+ {
+ *L1 = IL1;
+ }
+
+ return I1;
+}
+#if BOOST_WORKAROUND(BOOST_MSVC, < 1800)
+// Template argument dedcution failure otherwise:
+template<class F, class Real>
+Real trapezoidal(F f, Real a, Real b, Real tol = 0, std::size_t max_refinements = 10, Real* error_estimate = 0, Real* L1 = 0)
+#elif !defined(BOOST_NO_CXX11_NULLPTR)
+template<class F, class Real>
+Real trapezoidal(F f, Real a, Real b, Real tol = boost::math::tools::root_epsilon<Real>(), std::size_t max_refinements = 10, Real* error_estimate = nullptr, Real* L1 = nullptr)
+#else
+template<class F, class Real>
+Real trapezoidal(F f, Real a, Real b, Real tol = boost::math::tools::root_epsilon<Real>(), std::size_t max_refinements = 10, Real* error_estimate = 0, Real* L1 = 0)
+#endif
+{
+#if BOOST_WORKAROUND(BOOST_MSVC, <= 1600)
+ if (tol == 0)
+ tol = boost::math::tools::root_epsilon<Real>();
+#endif
+ return trapezoidal(f, a, b, tol, max_refinements, error_estimate, L1, boost::math::policies::policy<>());
+}
+
+}}}
+#endif
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>
diff --git a/boost/math/tools/polynomial.hpp b/boost/math/tools/polynomial.hpp
index 464c334d10..56c4e6895e 100644
--- a/boost/math/tools/polynomial.hpp
+++ b/boost/math/tools/polynomial.hpp
@@ -15,7 +15,6 @@
#include <boost/assert.hpp>
#include <boost/config.hpp>
-#include <boost/config/suffix.hpp>
#include <boost/function.hpp>
#include <boost/lambda/lambda.hpp>
#include <boost/math/tools/rational.hpp>
@@ -713,6 +712,11 @@ inline std::basic_ostream<charT, traits>& operator << (std::basic_ostream<charT,
} // namespace math
} // namespace boost
+//
+// Polynomial specific overload of gcd algorithm:
+//
+#include <boost/math/tools/polynomial_gcd.hpp>
+
#endif // BOOST_MATH_TOOLS_POLYNOMIAL_HPP
diff --git a/boost/math/tools/polynomial_gcd.hpp b/boost/math/tools/polynomial_gcd.hpp
new file mode 100644
index 0000000000..fdbafda6ca
--- /dev/null
+++ b/boost/math/tools/polynomial_gcd.hpp
@@ -0,0 +1,209 @@
+// (C) Copyright Jeremy William Murphy 2016.
+
+// 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_TOOLS_POLYNOMIAL_GCD_HPP
+#define BOOST_MATH_TOOLS_POLYNOMIAL_GCD_HPP
+
+#ifdef _MSC_VER
+#pragma once
+#endif
+
+#include <boost/math/tools/polynomial.hpp>
+#include <boost/math/common_factor_rt.hpp>
+#include <boost/type_traits/is_pod.hpp>
+
+
+namespace boost{
+
+ namespace integer {
+
+ namespace gcd_detail {
+
+ template <class T>
+ struct gcd_traits;
+
+ template <class T>
+ struct gcd_traits<boost::math::tools::polynomial<T> >
+ {
+ inline static const boost::math::tools::polynomial<T>& abs(const boost::math::tools::polynomial<T>& val) { return val; }
+
+ static const method_type method = method_euclid;
+ };
+
+ }
+}
+
+
+
+namespace math{ namespace tools{
+
+/* From Knuth, 4.6.1:
+*
+* We may write any nonzero polynomial u(x) from R[x] where R is a UFD as
+*
+* u(x) = cont(u) · pp(u(x))
+*
+* where cont(u), the content of u, is an element of S, and pp(u(x)), the primitive
+* part of u(x), is a primitive polynomial over S.
+* When u(x) = 0, it is convenient to define cont(u) = pp(u(x)) = O.
+*/
+
+template <class T>
+T content(polynomial<T> const &x)
+{
+ return x ? gcd_range(x.data().begin(), x.data().end()).first : T(0);
+}
+
+// Knuth, 4.6.1
+template <class T>
+polynomial<T> primitive_part(polynomial<T> const &x, T const &cont)
+{
+ return x ? x / cont : polynomial<T>();
+}
+
+
+template <class T>
+polynomial<T> primitive_part(polynomial<T> const &x)
+{
+ return primitive_part(x, content(x));
+}
+
+
+// Trivial but useful convenience function referred to simply as l() in Knuth.
+template <class T>
+T leading_coefficient(polynomial<T> const &x)
+{
+ return x ? x.data().back() : T(0);
+}
+
+
+namespace detail
+{
+ /* Reduce u and v to their primitive parts and return the gcd of their
+ * contents. Used in a couple of gcd algorithms.
+ */
+ template <class T>
+ T reduce_to_primitive(polynomial<T> &u, polynomial<T> &v)
+ {
+ using boost::math::gcd;
+ T const u_cont = content(u), v_cont = content(v);
+ u /= u_cont;
+ v /= v_cont;
+ return gcd(u_cont, v_cont);
+ }
+}
+
+
+/**
+* Knuth, The Art of Computer Programming: Volume 2, Third edition, 1998
+* Algorithm 4.6.1C: Greatest common divisor over a unique factorization domain.
+*
+* The subresultant algorithm by George E. Collins [JACM 14 (1967), 128-142],
+* later improved by W. S. Brown and J. F. Traub [JACM 18 (1971), 505-514].
+*
+* Although step C3 keeps the coefficients to a "reasonable" size, they are
+* still potentially several binary orders of magnitude larger than the inputs.
+* Thus, this algorithm should only be used where T is a multi-precision type.
+*
+* @tparam T Polynomial coefficient type.
+* @param u First polynomial.
+* @param v Second polynomial.
+* @return Greatest common divisor of polynomials u and v.
+*/
+template <class T>
+typename enable_if_c< std::numeric_limits<T>::is_integer, polynomial<T> >::type
+subresultant_gcd(polynomial<T> u, polynomial<T> v)
+{
+ using std::swap;
+ BOOST_ASSERT(u || v);
+
+ if (!u)
+ return v;
+ if (!v)
+ return u;
+
+ typedef typename polynomial<T>::size_type N;
+
+ if (u.degree() < v.degree())
+ swap(u, v);
+
+ T const d = detail::reduce_to_primitive(u, v);
+ T g = 1, h = 1;
+ polynomial<T> r;
+ while (true)
+ {
+ BOOST_ASSERT(u.degree() >= v.degree());
+ // Pseudo-division.
+ r = u % v;
+ if (!r)
+ return d * primitive_part(v); // Attach the content.
+ if (r.degree() == 0)
+ return d * polynomial<T>(T(1)); // The content is the result.
+ N const delta = u.degree() - v.degree();
+ // Adjust remainder.
+ u = v;
+ v = r / (g * detail::integer_power(h, delta));
+ g = leading_coefficient(u);
+ T const tmp = detail::integer_power(g, delta);
+ if (delta <= N(1))
+ h = tmp * detail::integer_power(h, N(1) - delta);
+ else
+ h = tmp / detail::integer_power(h, delta - N(1));
+ }
+}
+
+
+/**
+ * @brief GCD for polynomials with unbounded multi-precision integral coefficients.
+ *
+ * The multi-precision constraint is enforced via numeric_limits.
+ *
+ * Note that intermediate terms in the evaluation can grow arbitrarily large, hence the need for
+ * unbounded integers, otherwise numeric loverflow would break the algorithm.
+ *
+ * @tparam T A multi-precision integral type.
+ */
+template <typename T>
+typename enable_if_c<std::numeric_limits<T>::is_integer && !std::numeric_limits<T>::is_bounded, polynomial<T> >::type
+gcd(polynomial<T> const &u, polynomial<T> const &v)
+{
+ return subresultant_gcd(u, v);
+}
+// GCD over bounded integers is not currently allowed:
+template <typename T>
+typename enable_if_c<std::numeric_limits<T>::is_integer && std::numeric_limits<T>::is_bounded, polynomial<T> >::type
+gcd(polynomial<T> const &u, polynomial<T> const &v)
+{
+ BOOST_STATIC_ASSERT_MSG(sizeof(v) == 0, "GCD on polynomials of bounded integers is disallowed due to the excessive growth in the size of intermediate terms.");
+ return subresultant_gcd(u, v);
+}
+// GCD over polynomials of floats can go via the Euclid algorithm:
+template <typename T>
+typename enable_if_c<!std::numeric_limits<T>::is_integer && (std::numeric_limits<T>::min_exponent != std::numeric_limits<T>::max_exponent) && !std::numeric_limits<T>::is_exact, polynomial<T> >::type
+gcd(polynomial<T> const &u, polynomial<T> const &v)
+{
+ return boost::integer::gcd_detail::Euclid_gcd(u, v);
+}
+
+}
+//
+// Using declaration so we overload the default implementation in this namespace:
+//
+using boost::math::tools::gcd;
+
+}
+
+namespace integer
+{
+ //
+ // Using declaration so we overload the default implementation in this namespace:
+ //
+ using boost::math::tools::gcd;
+}
+
+} // namespace boost::math::tools
+
+#endif