diff options
author | DongHun Kwak <dh0128.kwak@samsung.com> | 2017-09-13 11:24:46 +0900 |
---|---|---|
committer | DongHun Kwak <dh0128.kwak@samsung.com> | 2017-09-13 11:25:39 +0900 |
commit | 4fadd968fa12130524c8380f33fcfe25d4de79e5 (patch) | |
tree | fd26a490cd15388d42fc6652b3c5c13012e7f93e /boost/math | |
parent | b5c87084afaef42b2d058f68091be31988a6a874 (diff) | |
download | boost-4fadd968fa12130524c8380f33fcfe25d4de79e5.tar.gz boost-4fadd968fa12130524c8380f33fcfe25d4de79e5.tar.bz2 boost-4fadd968fa12130524c8380f33fcfe25d4de79e5.zip |
Imported Upstream version 1.65.0upstream/1.65.0
Change-Id: Icf8400b375482cb11bcf77440a6934ba360d6ba4
Signed-off-by: DongHun Kwak <dh0128.kwak@samsung.com>
Diffstat (limited to 'boost/math')
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 |