diff options
Diffstat (limited to 'boost/math')
21 files changed, 856 insertions, 58 deletions
diff --git a/boost/math/common_factor.hpp b/boost/math/common_factor.hpp index 21a71e4948..0c7cf432a8 100644 --- a/boost/math/common_factor.hpp +++ b/boost/math/common_factor.hpp @@ -5,7 +5,7 @@ // accompanying file LICENSE_1_0.txt or copy at // http://www.boost.org/LICENSE_1_0.txt) -// See http://www.boost.org for updates, documentation, and revision history. +// See http://www.boost.org for updates, documentation, and revision history. #ifndef BOOST_MATH_COMMON_FACTOR_HPP #define BOOST_MATH_COMMON_FACTOR_HPP @@ -13,4 +13,7 @@ #include <boost/math/common_factor_ct.hpp> #include <boost/math/common_factor_rt.hpp> +BOOST_HEADER_DEPRECATED("<boost/integer/common_factor.hpp>"); + + #endif // BOOST_MATH_COMMON_FACTOR_HPP diff --git a/boost/math/common_factor_ct.hpp b/boost/math/common_factor_ct.hpp index 3ca0905945..2d0870e12a 100644 --- a/boost/math/common_factor_ct.hpp +++ b/boost/math/common_factor_ct.hpp @@ -5,12 +5,15 @@ // accompanying file LICENSE_1_0.txt or copy at // http://www.boost.org/LICENSE_1_0.txt) -// See http://www.boost.org for updates, documentation, and revision history. +// See http://www.boost.org for updates, documentation, and revision history. #ifndef BOOST_MATH_COMMON_FACTOR_CT_HPP #define BOOST_MATH_COMMON_FACTOR_CT_HPP #include <boost/integer/common_factor_ct.hpp> +#include <boost/config/header_deprecated.hpp> + +BOOST_HEADER_DEPRECATED("<boost/integer/common_factor_ct.hpp>"); namespace boost { diff --git a/boost/math/common_factor_rt.hpp b/boost/math/common_factor_rt.hpp index 42d9edfc04..c1f60aca15 100644 --- a/boost/math/common_factor_rt.hpp +++ b/boost/math/common_factor_rt.hpp @@ -8,17 +8,18 @@ #define BOOST_MATH_COMMON_FACTOR_RT_HPP #include <boost/integer/common_factor_rt.hpp> +#include <boost/config/header_deprecated.hpp> + +BOOST_HEADER_DEPRECATED("<boost/integer/common_factor_rt.hpp>"); namespace boost { namespace math { - 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; - } } diff --git a/boost/math/concepts/std_real_concept.hpp b/boost/math/concepts/std_real_concept.hpp index d855d2d55f..b116282888 100644 --- a/boost/math/concepts/std_real_concept.hpp +++ b/boost/math/concepts/std_real_concept.hpp @@ -226,6 +226,8 @@ inline boost::math::concepts::std_real_concept sqrt(boost::math::concepts::std_r { return std::sqrt(a.value()); } inline boost::math::concepts::std_real_concept tanh(boost::math::concepts::std_real_concept a) { return std::tanh(a.value()); } +inline boost::math::concepts::std_real_concept (nextafter)(boost::math::concepts::std_real_concept a, boost::math::concepts::std_real_concept b) +{ return (boost::math::nextafter)(a, b); } // // C++11 ism's // Note that these must not actually call the std:: versions as that precludes using this @@ -237,6 +239,10 @@ inline boost::math::concepts::std_real_concept acosh(boost::math::concepts::std_ { return boost::math::acosh(a.value(), boost::math::policies::make_policy(boost::math::policies::overflow_error<boost::math::policies::ignore_error>())); } inline boost::math::concepts::std_real_concept atanh(boost::math::concepts::std_real_concept a) { return boost::math::atanh(a.value(), boost::math::policies::make_policy(boost::math::policies::overflow_error<boost::math::policies::ignore_error>())); } +inline bool (isfinite)(boost::math::concepts::std_real_concept a) +{ + return (boost::math::isfinite)(a.value()); +} } // namespace std @@ -396,6 +402,13 @@ inline BOOST_MATH_CONSTEXPR int digits<concepts::std_real_concept>(BOOST_MATH_EX return digits<concepts::std_real_concept_base_type>(); } +template <> +inline double real_cast<double, concepts::std_real_concept>(concepts::std_real_concept r) +{ + return static_cast<double>(r.value()); +} + + } // namespace tools #if BOOST_WORKAROUND(BOOST_MSVC, <= 1310) diff --git a/boost/math/quadrature/detail/exp_sinh_detail.hpp b/boost/math/quadrature/detail/exp_sinh_detail.hpp index 3ce7a2d217..4d852393aa 100644 --- a/boost/math/quadrature/detail/exp_sinh_detail.hpp +++ b/boost/math/quadrature/detail/exp_sinh_detail.hpp @@ -228,7 +228,7 @@ Real exp_sinh_detail<Real, Policy>::integrate(const F& f, Real* error, Real* L1, //std::cout << "Estimate: " << I1 << " Error estimate at level " << i << " = " << err << std::endl; if (!isfinite(I1)) { - return policies::raise_evaluation_error(function, "The exp_sinh quadrature evaluated your function at a singular point and resulted in %1%. Please ensure your function evaluates to a finite number of its entire domain.", I1, Policy()); + return policies::raise_evaluation_error(function, "The exp_sinh quadrature evaluated your function at a singular point and returned %1%. Please ensure your function evaluates to a finite number over its entire domain.", I1, Policy()); } if (err <= tolerance*L1_I1) { diff --git a/boost/math/quadrature/naive_monte_carlo.hpp b/boost/math/quadrature/naive_monte_carlo.hpp new file mode 100644 index 0000000000..5fc2455024 --- /dev/null +++ b/boost/math/quadrature/naive_monte_carlo.hpp @@ -0,0 +1,426 @@ +/* + * Copyright Nick Thompson, 2018 + * 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_QUADRATURE_NAIVE_MONTE_CARLO_HPP +#define BOOST_MATH_QUADRATURE_NAIVE_MONTE_CARLO_HPP +#include <sstream> +#include <algorithm> +#include <vector> +#include <boost/atomic.hpp> +#include <functional> +#include <future> +#include <thread> +#include <initializer_list> +#include <utility> +#include <random> +#include <chrono> +#include <map> +#include <boost/math/policies/error_handling.hpp> + +namespace boost { namespace math { namespace quadrature { + +namespace detail { + enum class limit_classification {FINITE, + LOWER_BOUND_INFINITE, + UPPER_BOUND_INFINITE, + DOUBLE_INFINITE}; +} + +template<class Real, class F, class RandomNumberGenerator = std::mt19937_64, class Policy = boost::math::policies::policy<>> +class naive_monte_carlo +{ +public: + naive_monte_carlo(const F& integrand, + std::vector<std::pair<Real, Real>> const & bounds, + Real error_goal, + bool singular = true, + size_t threads = std::thread::hardware_concurrency(), + size_t seed = 0): m_num_threads{threads}, m_seed{seed} + { + using std::numeric_limits; + using std::sqrt; + size_t n = bounds.size(); + m_lbs.resize(n); + m_dxs.resize(n); + m_limit_types.resize(n); + m_volume = 1; + static const char* function = "boost::math::quadrature::naive_monte_carlo<%1%>"; + for (size_t i = 0; i < n; ++i) + { + if (bounds[i].second <= bounds[i].first) + { + boost::math::policies::raise_domain_error(function, "The upper bound is <= the lower bound.\n", bounds[i].second, Policy()); + return; + } + if (bounds[i].first == -numeric_limits<Real>::infinity()) + { + if (bounds[i].second == numeric_limits<Real>::infinity()) + { + m_limit_types[i] = detail::limit_classification::DOUBLE_INFINITE; + } + else + { + m_limit_types[i] = detail::limit_classification::LOWER_BOUND_INFINITE; + // Ok ok this is bad to use the second bound as the lower limit and then reflect. + m_lbs[i] = bounds[i].second; + m_dxs[i] = numeric_limits<Real>::quiet_NaN(); + } + } + else if (bounds[i].second == numeric_limits<Real>::infinity()) + { + m_limit_types[i] = detail::limit_classification::UPPER_BOUND_INFINITE; + if (singular) + { + // I've found that it's easier to sample on a closed set and perturb the boundary + // than to try to sample very close to the boundary. + m_lbs[i] = std::nextafter(bounds[i].first, (std::numeric_limits<Real>::max)()); + } + else + { + m_lbs[i] = bounds[i].first; + } + m_dxs[i] = numeric_limits<Real>::quiet_NaN(); + } + else + { + m_limit_types[i] = detail::limit_classification::FINITE; + if (singular) + { + if (bounds[i].first == 0) + { + m_lbs[i] = std::numeric_limits<Real>::epsilon(); + } + else + { + m_lbs[i] = std::nextafter(bounds[i].first, (std::numeric_limits<Real>::max)()); + } + + m_dxs[i] = std::nextafter(bounds[i].second, std::numeric_limits<Real>::lowest()) - m_lbs[i]; + } + else + { + m_lbs[i] = bounds[i].first; + m_dxs[i] = bounds[i].second - bounds[i].first; + } + m_volume *= m_dxs[i]; + } + } + + m_integrand = [this, &integrand](std::vector<Real> & x)->Real + { + Real coeff = m_volume; + for (size_t i = 0; i < x.size(); ++i) + { + // Variable transformation are listed at: + // https://en.wikipedia.org/wiki/Numerical_integration + // However, we've made some changes to these so that we can evaluate on a compact domain. + if (m_limit_types[i] == detail::limit_classification::FINITE) + { + x[i] = m_lbs[i] + x[i]*m_dxs[i]; + } + else if (m_limit_types[i] == detail::limit_classification::UPPER_BOUND_INFINITE) + { + Real t = x[i]; + Real z = 1/(1 + numeric_limits<Real>::epsilon() - t); + coeff *= (z*z)*(1 + numeric_limits<Real>::epsilon()); + x[i] = m_lbs[i] + t*z; + } + else if (m_limit_types[i] == detail::limit_classification::LOWER_BOUND_INFINITE) + { + Real t = x[i]; + Real z = 1/(t+sqrt((numeric_limits<Real>::min)())); + coeff *= (z*z); + x[i] = m_lbs[i] + (t-1)*z; + } + else + { + Real t1 = 1/(1+numeric_limits<Real>::epsilon() - x[i]); + Real t2 = 1/(x[i]+numeric_limits<Real>::epsilon()); + x[i] = (2*x[i]-1)*t1*t2/4; + coeff *= (t1*t1+t2*t2)/4; + } + } + return coeff*integrand(x); + }; + + // If we don't do a single function call in the constructor, + // we can't do a restart. + std::vector<Real> x(m_lbs.size()); + + // If the seed is zero, that tells us to choose a random seed for the user: + if (seed == 0) + { + std::random_device rd; + seed = rd(); + } + + RandomNumberGenerator gen(seed); + Real inv_denom = 1/static_cast<Real>(((gen.max)()-(gen.min)())); + + m_num_threads = (std::max)(m_num_threads, (size_t) 1); + m_thread_calls.reset(new boost::atomic<size_t>[threads]); + m_thread_Ss.reset(new boost::atomic<Real>[threads]); + m_thread_averages.reset(new boost::atomic<Real>[threads]); + + Real avg = 0; + for (size_t i = 0; i < m_num_threads; ++i) + { + for (size_t j = 0; j < m_lbs.size(); ++j) + { + x[j] = (gen()-(gen.min)())*inv_denom; + } + Real y = m_integrand(x); + m_thread_averages[i] = y; // relaxed store + m_thread_calls[i] = 1; + m_thread_Ss[i] = 0; + avg += y; + } + avg /= m_num_threads; + m_avg = avg; // relaxed store + + m_error_goal = error_goal; // relaxed store + m_start = std::chrono::system_clock::now(); + m_done = false; // relaxed store + m_total_calls = m_num_threads; // relaxed store + m_variance = (numeric_limits<Real>::max)(); + } + + std::future<Real> integrate() + { + // Set done to false in case we wish to restart: + m_done.store(false); // relaxed store, no worker threads yet + return std::async(std::launch::async, + &naive_monte_carlo::m_integrate, this); + } + + void cancel() + { + // If seed = 0 (meaning have the routine pick the seed), this leaves the seed the same. + // If seed != 0, then the seed is changed, so a restart doesn't do the exact same thing. + m_seed = m_seed*m_seed; + m_done = true; // relaxed store, worker threads will get the message eventually + } + + Real variance() const + { + return m_variance.load(); + } + + Real current_error_estimate() const + { + using std::sqrt; + // + // There is a bug here: m_variance and m_total_calls get updated asynchronously + // and may be out of synch when we compute the error estimate, not sure if it matters though... + // + return sqrt(m_variance.load()/m_total_calls.load()); + } + + std::chrono::duration<Real> estimated_time_to_completion() const + { + auto now = std::chrono::system_clock::now(); + std::chrono::duration<Real> elapsed_seconds = now - m_start; + Real r = this->current_error_estimate()/m_error_goal.load(); // relaxed load + if (r*r <= 1) { + return 0*elapsed_seconds; + } + return (r*r - 1)*elapsed_seconds; + } + + void update_target_error(Real new_target_error) + { + m_error_goal = new_target_error; // relaxed store + } + + Real progress() const + { + Real r = m_error_goal.load()/this->current_error_estimate(); // relaxed load + if (r*r >= 1) + { + return 1; + } + return r*r; + } + + Real current_estimate() const + { + return m_avg.load(); + } + + size_t calls() const + { + return m_total_calls.load(); // relaxed load + } + +private: + + Real m_integrate() + { + m_start = std::chrono::system_clock::now(); + std::vector<std::thread> threads(m_num_threads); + size_t seed; + // If the user tells us to pick a seed, pick a seed: + if (m_seed == 0) + { + std::random_device rd; + seed = rd(); + } + else // use the seed we are given: + { + seed = m_seed; + } + RandomNumberGenerator gen(seed); + for (size_t i = 0; i < threads.size(); ++i) + { + threads[i] = std::thread(&naive_monte_carlo::m_thread_monte, this, i, gen()); + } + do { + std::this_thread::sleep_for(std::chrono::milliseconds(100)); + size_t total_calls = 0; + for (size_t i = 0; i < m_num_threads; ++i) + { + size_t t_calls = m_thread_calls[i].load(boost::memory_order::consume); + total_calls += t_calls; + } + Real variance = 0; + Real avg = 0; + for (size_t i = 0; i < m_num_threads; ++i) + { + size_t t_calls = m_thread_calls[i].load(boost::memory_order::consume); + // Will this overflow? Not hard to remove . . . + avg += m_thread_averages[i].load(boost::memory_order::relaxed)*( (Real) t_calls/ (Real) total_calls); + variance += m_thread_Ss[i].load(boost::memory_order::relaxed); + } + m_avg.store(avg, boost::memory_order::release); + m_variance.store(variance/(total_calls - 1), boost::memory_order::release); + m_total_calls = total_calls; // relaxed store, it's just for user feedback + // Allow cancellation: + if (m_done) // relaxed load + { + break; + } + } while (this->current_error_estimate() > m_error_goal.load(boost::memory_order::consume)); + // Error bound met; signal the threads: + m_done = true; // relaxed store, threads will get the message in the end + std::for_each(threads.begin(), threads.end(), + std::mem_fn(&std::thread::join)); + if (m_exception) + { + std::rethrow_exception(m_exception); + } + // Incorporate their work into the final estimate: + size_t total_calls = 0; + for (size_t i = 0; i < m_num_threads; ++i) + { + size_t t_calls = m_thread_calls[i].load(boost::memory_order::consume); + total_calls += t_calls; + } + Real variance = 0; + Real avg = 0; + + for (size_t i = 0; i < m_num_threads; ++i) + { + size_t t_calls = m_thread_calls[i].load(boost::memory_order::consume); + // Averages weighted by the number of calls the thread made: + avg += m_thread_averages[i].load(boost::memory_order::relaxed)*( (Real) t_calls/ (Real) total_calls); + variance += m_thread_Ss[i].load(boost::memory_order::relaxed); + } + m_avg.store(avg, boost::memory_order::release); + m_variance.store(variance/(total_calls - 1), boost::memory_order::release); + m_total_calls = total_calls; // relaxed store, this is just user feedback + + return m_avg.load(boost::memory_order::consume); + } + + void m_thread_monte(size_t thread_index, size_t seed) + { + using std::numeric_limits; + try + { + std::vector<Real> x(m_lbs.size()); + RandomNumberGenerator gen(seed); + Real inv_denom = (Real) 1/(Real)( (gen.max)() - (gen.min)() ); + Real M1 = m_thread_averages[thread_index].load(boost::memory_order::consume); + Real S = m_thread_Ss[thread_index].load(boost::memory_order::consume); + // Kahan summation is required or the value of the integrand will go on a random walk during long computations. + // See the implementation discussion. + // The idea is that the unstabilized additions have error sigma(f)/sqrt(N) + epsilon*N, which diverges faster than it converges! + // Kahan summation turns this to sigma(f)/sqrt(N) + epsilon^2*N, and the random walk occurs on a timescale of 10^14 years (on current hardware) + Real compensator = 0; + size_t k = m_thread_calls[thread_index].load(boost::memory_order::consume); + while (!m_done) // relaxed load + { + int j = 0; + // If we don't have a certain number of calls before an update, we can easily terminate prematurely + // because the variance estimate is way too low. This magic number is a reasonable compromise, as 1/sqrt(2048) = 0.02, + // so it should recover 2 digits if the integrand isn't poorly behaved, and if it is, it should discover that before premature termination. + // Of course if the user has 64 threads, then this number is probably excessive. + int magic_calls_before_update = 2048; + while (j++ < magic_calls_before_update) + { + for (size_t i = 0; i < m_lbs.size(); ++i) + { + x[i] = (gen() - (gen.min)())*inv_denom; + } + Real f = m_integrand(x); + using std::isfinite; + if (!isfinite(f)) + { + // The call to m_integrand transform x, so this error message states the correct node. + std::stringstream os; + os << "Your integrand was evaluated at {"; + for (size_t i = 0; i < x.size() -1; ++i) + { + os << x[i] << ", "; + } + os << x[x.size() -1] << "}, and returned " << f << std::endl; + static const char* function = "boost::math::quadrature::naive_monte_carlo<%1%>"; + boost::math::policies::raise_domain_error(function, os.str().c_str(), /*this is a dummy arg to make it compile*/ 7.2, Policy()); + } + ++k; + Real term = (f - M1)/k; + Real y1 = term - compensator; + Real M2 = M1 + y1; + compensator = (M2 - M1) - y1; + S += (f - M1)*(f - M2); + M1 = M2; + } + m_thread_averages[thread_index].store(M1, boost::memory_order::release); + m_thread_Ss[thread_index].store(S, boost::memory_order::release); + m_thread_calls[thread_index].store(k, boost::memory_order::release); + } + } + catch (...) + { + // Signal the other threads that the computation is ruined: + m_done = true; // relaxed store + m_exception = std::current_exception(); + } + } + + std::function<Real(std::vector<Real> &)> m_integrand; + size_t m_num_threads; + size_t m_seed; + boost::atomic<Real> m_error_goal; + boost::atomic<bool> m_done; + std::vector<Real> m_lbs; + std::vector<Real> m_dxs; + std::vector<detail::limit_classification> m_limit_types; + Real m_volume; + boost::atomic<size_t> m_total_calls; + // I wanted these to be vectors rather than maps, + // but you can't resize a vector of atomics. + std::unique_ptr<boost::atomic<size_t>[]> m_thread_calls; + boost::atomic<Real> m_variance; + std::unique_ptr<boost::atomic<Real>[]> m_thread_Ss; + boost::atomic<Real> m_avg; + std::unique_ptr<boost::atomic<Real>[]> m_thread_averages; + std::chrono::time_point<std::chrono::system_clock> m_start; + std::exception_ptr m_exception; +}; + +}}} +#endif diff --git a/boost/math/quaternion.hpp b/boost/math/quaternion.hpp index 7b7a1aa17d..a69147c539 100644 --- a/boost/math/quaternion.hpp +++ b/boost/math/quaternion.hpp @@ -11,6 +11,7 @@ #define BOOST_QUATERNION_HPP #include <boost/config.hpp> // for BOOST_NO_STD_LOCALE +#include <boost/math_fwd.hpp> #include <boost/detail/workaround.hpp> #include <boost/type_traits/is_convertible.hpp> #include <boost/utility/enable_if.hpp> diff --git a/boost/math/special_functions/chebyshev_transform.hpp b/boost/math/special_functions/chebyshev_transform.hpp index d4ce106d72..e50c40bd92 100644 --- a/boost/math/special_functions/chebyshev_transform.hpp +++ b/boost/math/special_functions/chebyshev_transform.hpp @@ -227,16 +227,6 @@ public: return dzdx*(z*d1 - d2 + b1); } - void print_coefficients() const - { - std::cout << "{"; - for(auto const & coeff : m_coeffs) { - std::cout << coeff << ", "; - } - std::cout << "}\n"; - } - - private: std::vector<Real> m_coeffs; Real m_a; diff --git a/boost/math/special_functions/detail/bernoulli_details.hpp b/boost/math/special_functions/detail/bernoulli_details.hpp index 41a59e53c6..75fadbf34a 100644 --- a/boost/math/special_functions/detail/bernoulli_details.hpp +++ b/boost/math/special_functions/detail/bernoulli_details.hpp @@ -192,9 +192,18 @@ struct fixed_vector : private std::allocator<T> } ~fixed_vector() { +#ifdef BOOST_NO_CXX11_ALLOCATOR for(unsigned i = 0; i < m_used; ++i) this->destroy(&m_data[i]); this->deallocate(m_data, m_capacity); +#else + typedef std::allocator<T> allocator_type; + typedef std::allocator_traits<allocator_type> allocator_traits; + allocator_type& alloc = *this; + for(unsigned i = 0; i < m_used; ++i) + allocator_traits::destroy(alloc, &m_data[i]); + allocator_traits::deallocate(alloc, m_data, m_capacity); +#endif } T& operator[](unsigned n) { BOOST_ASSERT(n < m_used); return m_data[n]; } const T& operator[](unsigned n)const { BOOST_ASSERT(n < m_used); return m_data[n]; } diff --git a/boost/math/special_functions/detail/fp_traits.hpp b/boost/math/special_functions/detail/fp_traits.hpp index 09dc5169aa..c957022223 100644 --- a/boost/math/special_functions/detail/fp_traits.hpp +++ b/boost/math/special_functions/detail/fp_traits.hpp @@ -556,7 +556,8 @@ struct select_native<long double> && !defined(__FAST_MATH__)\ && !defined(BOOST_MATH_DISABLE_STD_FPCLASSIFY)\ && !defined(BOOST_INTEL)\ - && !defined(sun) + && !defined(sun)\ + && !defined(__VXWORKS__) # define BOOST_MATH_USE_STD_FPCLASSIFY #endif diff --git a/boost/math/special_functions/detail/unchecked_bernoulli.hpp b/boost/math/special_functions/detail/unchecked_bernoulli.hpp index 03c376678d..5a16d9df17 100644 --- a/boost/math/special_functions/detail/unchecked_bernoulli.hpp +++ b/boost/math/special_functions/detail/unchecked_bernoulli.hpp @@ -19,6 +19,13 @@ #include <boost/mpl/int.hpp> #include <boost/type_traits/is_convertible.hpp> +#ifdef BOOST_MATH_HAVE_CONSTEXPR_TABLES +#include <array> +#else +#include <boost/array.hpp> +#endif + + namespace boost { namespace math { namespace detail { @@ -83,9 +90,13 @@ struct max_bernoulli_b2n : public detail::max_bernoulli_index<detail::bernoulli_ namespace detail{ template <class T> -inline T unchecked_bernoulli_imp(std::size_t n, const mpl::int_<0>& ) +inline BOOST_MATH_CONSTEXPR_TABLE_FUNCTION T unchecked_bernoulli_imp(std::size_t n, const mpl::int_<0>& ) { +#ifdef BOOST_MATH_HAVE_CONSTEXPR_TABLES + constexpr std::array<boost::int64_t, 1 + max_bernoulli_b2n<T>::value> numerators = +#else static const boost::array<boost::int64_t, 1 + max_bernoulli_b2n<T>::value> numerators = +#endif {{ boost::int64_t( +1LL), boost::int64_t( +1LL), @@ -107,7 +118,11 @@ inline T unchecked_bernoulli_imp(std::size_t n, const mpl::int_<0>& ) boost::int64_t(+2577687858367LL) }}; +#ifdef BOOST_MATH_HAVE_CONSTEXPR_TABLES + constexpr std::array<boost::int64_t, 1 + max_bernoulli_b2n<T>::value> denominators = +#else static const boost::array<boost::int64_t, 1 + max_bernoulli_b2n<T>::value> denominators = +#endif {{ boost::int64_t( 1LL), boost::int64_t( 6LL), @@ -132,9 +147,13 @@ inline T unchecked_bernoulli_imp(std::size_t n, const mpl::int_<0>& ) } template <class T> -inline T unchecked_bernoulli_imp(std::size_t n, const mpl::int_<1>& ) +inline BOOST_MATH_CONSTEXPR_TABLE_FUNCTION T unchecked_bernoulli_imp(std::size_t n, const mpl::int_<1>& ) { +#ifdef BOOST_MATH_HAVE_CONSTEXPR_TABLES + constexpr std::array<float, 1 + max_bernoulli_b2n<T>::value> bernoulli_data = +#else static const boost::array<float, 1 + max_bernoulli_b2n<T>::value> bernoulli_data = +#endif {{ +1.00000000000000000000000000000000000000000F, +0.166666666666666666666666666666666666666667F, @@ -176,9 +195,13 @@ inline T unchecked_bernoulli_imp(std::size_t n, const mpl::int_<1>& ) template <class T> -inline T unchecked_bernoulli_imp(std::size_t n, const mpl::int_<2>& ) +inline BOOST_MATH_CONSTEXPR_TABLE_FUNCTION T unchecked_bernoulli_imp(std::size_t n, const mpl::int_<2>& ) { +#ifdef BOOST_MATH_HAVE_CONSTEXPR_TABLES + constexpr std::array<double, 1 + max_bernoulli_b2n<T>::value> bernoulli_data = +#else static const boost::array<double, 1 + max_bernoulli_b2n<T>::value> bernoulli_data = +#endif {{ +1.00000000000000000000000000000000000000000, +0.166666666666666666666666666666666666666667, @@ -316,9 +339,13 @@ inline T unchecked_bernoulli_imp(std::size_t n, const mpl::int_<2>& ) } template <class T> -inline T unchecked_bernoulli_imp(std::size_t n, const mpl::int_<3>& ) +inline BOOST_MATH_CONSTEXPR_TABLE_FUNCTION T unchecked_bernoulli_imp(std::size_t n, const mpl::int_<3>& ) { +#ifdef BOOST_MATH_HAVE_CONSTEXPR_TABLES + constexpr std::array<long double, 1 + max_bernoulli_b2n<T>::value> bernoulli_data = +#else static const boost::array<long double, 1 + max_bernoulli_b2n<T>::value> bernoulli_data = +#endif {{ +1.00000000000000000000000000000000000000000L, +0.166666666666666666666666666666666666666667L, @@ -688,7 +715,7 @@ inline T unchecked_bernoulli_imp(std::size_t n, const mpl::int_<4>& ) } // namespace detail template<class T> -inline T unchecked_bernoulli_b2n(const std::size_t n) +inline BOOST_MATH_CONSTEXPR_TABLE_FUNCTION T unchecked_bernoulli_b2n(const std::size_t n) { typedef mpl::int_<detail::bernoulli_imp_variant<T>::value> tag_type; diff --git a/boost/math/special_functions/detail/unchecked_factorial.hpp b/boost/math/special_functions/detail/unchecked_factorial.hpp index ca9a752291..17366742c4 100644 --- a/boost/math/special_functions/detail/unchecked_factorial.hpp +++ b/boost/math/special_functions/detail/unchecked_factorial.hpp @@ -10,7 +10,6 @@ #pragma once #endif -#include <boost/array.hpp> #ifdef BOOST_MSVC #pragma warning(push) // Temporary until lexical cast fixed. #pragma warning(disable: 4127 4701) @@ -21,9 +20,15 @@ #ifdef BOOST_MSVC #pragma warning(pop) #endif -#include <boost/config/no_tr1/cmath.hpp> +#include <cmath> #include <boost/math/special_functions/math_fwd.hpp> +#ifdef BOOST_MATH_HAVE_CONSTEXPR_TABLES +#include <array> +#else +#include <boost/array.hpp> +#endif + namespace boost { namespace math { // Forward declarations: @@ -32,9 +37,13 @@ struct max_factorial; // Definitions: template <> -inline float unchecked_factorial<float>(unsigned i BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE_SPEC(float)) +inline BOOST_MATH_CONSTEXPR_TABLE_FUNCTION float unchecked_factorial<float>(unsigned i BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE_SPEC(float)) { +#ifdef BOOST_MATH_HAVE_CONSTEXPR_TABLES + constexpr std::array<float, 35> factorials = { { +#else static const boost::array<float, 35> factorials = {{ +#endif 1.0F, 1.0F, 2.0F, @@ -83,9 +92,13 @@ struct max_factorial<float> template <> -inline long double unchecked_factorial<long double>(unsigned i BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE_SPEC(long double)) +inline BOOST_MATH_CONSTEXPR_TABLE_FUNCTION long double unchecked_factorial<long double>(unsigned i BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE_SPEC(long double)) { +#ifdef BOOST_MATH_HAVE_CONSTEXPR_TABLES + constexpr std::array<long double, 171> factorials = { { +#else static const boost::array<long double, 171> factorials = {{ +#endif 1L, 1L, 2L, @@ -271,9 +284,13 @@ struct max_factorial<long double> #ifdef BOOST_MATH_USE_FLOAT128 template <> -inline BOOST_MATH_FLOAT128_TYPE unchecked_factorial<BOOST_MATH_FLOAT128_TYPE>(unsigned i) +inline BOOST_MATH_CONSTEXPR_TABLE_FUNCTION BOOST_MATH_FLOAT128_TYPE unchecked_factorial<BOOST_MATH_FLOAT128_TYPE>(unsigned i) { +#ifdef BOOST_MATH_HAVE_CONSTEXPR_TABLES + constexpr std::array<BOOST_MATH_FLOAT128_TYPE, 171> factorials = { { +#else static const boost::array<BOOST_MATH_FLOAT128_TYPE, 171> factorials = { { +#endif 1, 1, 2, @@ -459,7 +476,7 @@ struct max_factorial<BOOST_MATH_FLOAT128_TYPE> #endif template <> -inline double unchecked_factorial<double>(unsigned i BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE_SPEC(double)) +inline BOOST_MATH_CONSTEXPR_TABLE_FUNCTION double unchecked_factorial<double>(unsigned i BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE_SPEC(double)) { return static_cast<double>(boost::math::unchecked_factorial<long double>(i)); } diff --git a/boost/math/special_functions/ellint_d.hpp b/boost/math/special_functions/ellint_d.hpp index bc5a4b2a56..fa5c53db18 100644 --- a/boost/math/special_functions/ellint_d.hpp +++ b/boost/math/special_functions/ellint_d.hpp @@ -82,8 +82,12 @@ T ellint_d_imp(T phi, T k, const Policy& pol) s = -1; rphi = constants::half_pi<T>() - rphi; } + BOOST_MATH_INSTRUMENT_VARIABLE(rphi); + BOOST_MATH_INSTRUMENT_VARIABLE(m); T sinp = sin(rphi); T cosp = cos(rphi); + BOOST_MATH_INSTRUMENT_VARIABLE(sinp); + BOOST_MATH_INSTRUMENT_VARIABLE(cosp); T c = 1 / (sinp * sinp); T cm1 = cosp * cosp / (sinp * sinp); // c - 1 T k2 = k * k; @@ -99,6 +103,7 @@ T ellint_d_imp(T phi, T k, const Policy& pol) { // http://dlmf.nist.gov/19.25#E10 result = s * ellint_rd_imp(cm1, T(c - k2), c, pol) / 3; + BOOST_MATH_INSTRUMENT_VARIABLE(result); } if(m != 0) result += m * ellint_d_imp(k, pol); diff --git a/boost/math/special_functions/ellint_rd.hpp b/boost/math/special_functions/ellint_rd.hpp index 03b73b159f..c08430d545 100644 --- a/boost/math/special_functions/ellint_rd.hpp +++ b/boost/math/special_functions/ellint_rd.hpp @@ -133,6 +133,7 @@ T ellint_rd_imp(T x, T y, T z, const Policy& pol) T A0 = An; // This has an extra 1.2 fudge factor which is really only needed when x, y and z are close in magnitude: T Q = pow(tools::epsilon<T>() / 4, -T(1) / 8) * (std::max)((std::max)(An - x, An - y), An - z) * 1.2f; + BOOST_MATH_INSTRUMENT_VARIABLE(Q); T lambda, rx, ry, rz; unsigned k = 0; T fn = 1; @@ -151,6 +152,9 @@ T ellint_rd_imp(T x, T y, T z, const Policy& pol) zn = (zn + lambda) / 4; fn /= 4; Q /= 4; + BOOST_MATH_INSTRUMENT_VARIABLE(k); + BOOST_MATH_INSTRUMENT_VARIABLE(RD_sum); + BOOST_MATH_INSTRUMENT_VARIABLE(Q); if(Q < An) break; } @@ -168,6 +172,7 @@ T ellint_rd_imp(T x, T y, T z, const Policy& pol) T result = fn * pow(An, T(-3) / 2) * (1 - 3 * E2 / 14 + E3 / 6 + 9 * E2 * E2 / 88 - 3 * E4 / 22 - 9 * E2 * E3 / 52 + 3 * E5 / 26 - E2 * E2 * E2 / 16 + 3 * E3 * E3 / 40 + 3 * E2 * E4 / 20 + 45 * E2 * E2 * E3 / 272 - 9 * (E3 * E4 + E2 * E5) / 68); + BOOST_MATH_INSTRUMENT_VARIABLE(result); result += 3 * RD_sum; return result; diff --git a/boost/math/special_functions/lanczos.hpp b/boost/math/special_functions/lanczos.hpp index c9f4f90dcf..c1ff86930b 100644 --- a/boost/math/special_functions/lanczos.hpp +++ b/boost/math/special_functions/lanczos.hpp @@ -1284,7 +1284,7 @@ struct lanczos } // namespace boost #if !defined(_CRAYC) && !defined(__CUDACC__) && (!defined(__GNUC__) || (__GNUC__ > 3) || ((__GNUC__ == 3) && (__GNUC_MINOR__ > 3))) -#if (defined(_M_IX86_FP) && (_M_IX86_FP >= 2)) || defined(__SSE2__) || defined(_M_AMD64) || defined(_M_X64) +#if ((defined(_M_IX86_FP) && (_M_IX86_FP >= 2)) || defined(__SSE2__) || defined(_M_AMD64) || defined(_M_X64)) && !defined(_MANAGED) #include <boost/math/special_functions/detail/lanczos_sse2.hpp> #endif #endif diff --git a/boost/math/special_functions/math_fwd.hpp b/boost/math/special_functions/math_fwd.hpp index 03c55bd1a4..4f44f56113 100644 --- a/boost/math/special_functions/math_fwd.hpp +++ b/boost/math/special_functions/math_fwd.hpp @@ -1057,7 +1057,7 @@ namespace boost typename tools::promote_args<T, U>::type epsilon_difference(const T&, const U&); template<class T> - T unchecked_bernoulli_b2n(const std::size_t n); + BOOST_MATH_CONSTEXPR_TABLE_FUNCTION T unchecked_bernoulli_b2n(const std::size_t n); template <class T, class Policy> T bernoulli_b2n(const int i, const Policy &pol); template <class T> diff --git a/boost/math/special_functions/next.hpp b/boost/math/special_functions/next.hpp index 606b356542..a63983e1c3 100644 --- a/boost/math/special_functions/next.hpp +++ b/boost/math/special_functions/next.hpp @@ -30,7 +30,7 @@ namespace boost{ namespace math{ namespace concepts { class real_concept; - struct std_real_concept; + class std_real_concept; } diff --git a/boost/math/special_functions/prime.hpp b/boost/math/special_functions/prime.hpp index 94c28f9842..858d96d10b 100644 --- a/boost/math/special_functions/prime.hpp +++ b/boost/math/special_functions/prime.hpp @@ -8,15 +8,19 @@ #ifndef BOOST_MATH_SF_PRIME_HPP #define BOOST_MATH_SF_PRIME_HPP -#include <boost/array.hpp> #include <boost/cstdint.hpp> #include <boost/math/policies/error_handling.hpp> #include <boost/math/special_functions/math_fwd.hpp> +#ifdef BOOST_MATH_HAVE_CONSTEXPR_TABLES +#include <array> +#else +#include <boost/array.hpp> +#endif namespace boost{ namespace math{ template <class Policy> - boost::uint32_t prime(unsigned n, const Policy& pol) + BOOST_MATH_CONSTEXPR_TABLE_FUNCTION boost::uint32_t prime(unsigned n, const Policy& pol) { // // This is basically three big tables which together @@ -26,10 +30,17 @@ namespace boost{ namespace math{ // That gives us the first 10000 primes with the largest // being 104729: // +#ifdef BOOST_MATH_HAVE_CONSTEXPR_TABLES + constexpr unsigned b1 = 53; + constexpr unsigned b2 = 6541; + constexpr unsigned b3 = 10000; + constexpr std::array<unsigned char, 54> a1 = {{ +#else static const unsigned b1 = 53; static const unsigned b2 = 6541; static const unsigned b3 = 10000; static const boost::array<unsigned char, 54> a1 = {{ +#endif 2u, 3u, 5u, 7u, 11u, 13u, 17u, 19u, 23u, 29u, 31u, 37u, 41u, 43u, 47u, 53u, 59u, 61u, 67u, 71u, 73u, 79u, 83u, 89u, 97u, 101u, 103u, 107u, 109u, 113u, @@ -37,7 +48,11 @@ namespace boost{ namespace math{ 167u, 173u, 179u, 181u, 191u, 193u, 197u, 199u, 211u, 223u, 227u, 229u, 233u, 239u, 241u, 251u }}; +#ifdef BOOST_MATH_HAVE_CONSTEXPR_TABLES + constexpr std::array<boost::uint16_t, 6488> a2 = {{ +#else static const boost::array<boost::uint16_t, 6488> a2 = {{ +#endif 257u, 263u, 269u, 271u, 277u, 281u, 283u, 293u, 307u, 311u, 313u, 317u, 331u, 337u, 347u, 349u, 353u, 359u, 367u, 373u, 379u, 383u, 389u, 397u, 401u, 409u, @@ -760,7 +775,11 @@ namespace boost{ namespace math{ 65323u, 65327u, 65353u, 65357u, 65371u, 65381u, 65393u, 65407u, 65413u, 65419u, 65423u, 65437u, 65447u, 65449u, 65479u, 65497u, 65519u, 65521u }}; +#ifdef BOOST_MATH_HAVE_CONSTEXPR_TABLES + constexpr std::array<boost::uint16_t, 3458> a3 = {{ +#else static const boost::array<boost::uint16_t, 3458> a3 = {{ +#endif 2u, 4u, 8u, 16u, 22u, 28u, 44u, 46u, 52u, 64u, 74u, 82u, 94u, 98u, 112u, 116u, 122u, 142u, 152u, 164u, 166u, 172u, 178u, @@ -1208,7 +1227,7 @@ namespace boost{ namespace math{ return static_cast<boost::uint32_t>(a3[n - b2 - 1]) + 0xFFFFu; } - inline boost::uint32_t prime(unsigned n) + inline BOOST_MATH_CONSTEXPR_TABLE_FUNCTION boost::uint32_t prime(unsigned n) { return boost::math::prime(n, boost::math::policies::policy<>()); } diff --git a/boost/math/tools/config.hpp b/boost/math/tools/config.hpp index 8131facb98..016211cfa8 100644 --- a/boost/math/tools/config.hpp +++ b/boost/math/tools/config.hpp @@ -31,7 +31,7 @@ #if (defined(__CYGWIN__) || defined(__FreeBSD__) || defined(__NetBSD__) \ || (defined(__hppa) && !defined(__OpenBSD__)) || (defined(__NO_LONG_DOUBLE_MATH) && (DBL_MANT_DIG != LDBL_MANT_DIG))) \ && !defined(BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS) -//# define BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS +# define BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS #endif #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582)) // @@ -451,6 +451,17 @@ namespace boost{ namespace math{ # define BOOST_MATH_THREAD_LOCAL #endif +// +// Can we have constexpr tables? +// +#if (!defined(BOOST_NO_CXX11_HDR_ARRAY) && !defined(BOOST_NO_CXX14_CONSTEXPR)) || BOOST_WORKAROUND(BOOST_MSVC, >= 1910) +#define BOOST_MATH_HAVE_CONSTEXPR_TABLES +#define BOOST_MATH_CONSTEXPR_TABLE_FUNCTION constexpr +#else +#define BOOST_MATH_CONSTEXPR_TABLE_FUNCTION +#endif + + #endif // BOOST_MATH_TOOLS_CONFIG_HPP diff --git a/boost/math/tools/numerical_differentiation.hpp b/boost/math/tools/numerical_differentiation.hpp new file mode 100644 index 0000000000..34fef0db87 --- /dev/null +++ b/boost/math/tools/numerical_differentiation.hpp @@ -0,0 +1,267 @@ +// (C) Copyright Nick Thompson 2018. +// 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_NUMERICAL_DIFFERENTIATION_HPP +#define BOOST_MATH_TOOLS_NUMERICAL_DIFFERENTIATION_HPP + +/* + * Performs numerical differentiation by finite-differences. + * + * All numerical differentiation using finite-differences are ill-conditioned, and these routines are no exception. + * A simple argument demonstrates that the error is unbounded as h->0. + * Take the one sides finite difference formula f'(x) = (f(x+h)-f(x))/h. + * The evaluation of f induces an error as well as the error from the finite-difference approximation, giving + * |f'(x) - (f(x+h) -f(x))/h| < h|f''(x)|/2 + (|f(x)|+|f(x+h)|)eps/h =: g(h), where eps is the unit roundoff for the type. + * It is reasonable to choose h in a way that minimizes the maximum error bound g(h). + * The value of h that minimizes g is h = sqrt(2eps(|f(x)| + |f(x+h)|)/|f''(x)|), and for this value of h the error bound is + * sqrt(2eps(|f(x+h) +f(x)||f''(x)|)). + * In fact it is not necessary to compute the ratio (|f(x+h)| + |f(x)|)/|f''(x)|; the error bound of ~\sqrt{\epsilon} still holds if we set it to one. + * + * + * For more details on this method of analysis, see + * + * http://www.uio.no/studier/emner/matnat/math/MAT-INF1100/h08/kompendiet/diffint.pdf + * http://web.archive.org/web/20150420195907/http://www.uio.no/studier/emner/matnat/math/MAT-INF1100/h08/kompendiet/diffint.pdf + * + * + * It can be shown on general grounds that when choosing the optimal h, the maximum error in f'(x) is ~(|f(x)|eps)^k/k+1|f^(k-1)(x)|^1/k+1. + * From this we can see that full precision can be recovered in the limit k->infinity. + * + * References: + * + * 1) Fornberg, Bengt. "Generation of finite difference formulas on arbitrarily spaced grids." Mathematics of computation 51.184 (1988): 699-706. + * + * + * The second algorithm, the complex step derivative, is not ill-conditioned. + * However, it requires that your function can be evaluated at complex arguments. + * The idea is that f(x+ih) = f(x) +ihf'(x) - h^2f''(x) + ... so f'(x) \approx Im[f(x+ih)]/h. + * No subtractive cancellation occurs. The error is ~ eps|f'(x)| + eps^2|f'''(x)|/6; hard to beat that. + * + * References: + * + * 1) Squire, William, and George Trapp. "Using complex variables to estimate derivatives of real functions." Siam Review 40.1 (1998): 110-112. + */ + +#include <complex> +#include <boost/math/special_functions/next.hpp> + +// Make sure everyone is informed that C++17 is required: +namespace boost{ namespace math{ namespace tools { + +namespace detail { + template<class Real> + Real make_xph_representable(Real x, Real h) + { + using std::numeric_limits; + // Redefine h so that x + h is representable. Not using this trick leads to large error. + // The compiler flag -ffast-math evaporates these operations . . . + Real temp = x + h; + h = temp - x; + // Handle the case x + h == x: + if (h == 0) + { + h = boost::math::nextafter(x, (numeric_limits<Real>::max)()) - x; + } + return h; + } +} + +template<class F, class Real> +Real complex_step_derivative(const F f, Real x) +{ + // Is it really this easy? Yes. + // Note that some authors recommend taking the stepsize h to be smaller than epsilon(), some recommending use of the min(). + // This idea was tested over a few billion test cases and found the make the error *much* worse. + // Even 2eps and eps/2 made the error worse, which was surprising. + using std::complex; + using std::numeric_limits; + constexpr const Real step = (numeric_limits<Real>::epsilon)(); + constexpr const Real inv_step = 1/(numeric_limits<Real>::epsilon)(); + return f(complex<Real>(x, step)).imag()*inv_step; +} + +namespace detail { + + template <unsigned> + struct fd_tag {}; + + template<class F, class Real> + Real finite_difference_derivative(const F f, Real x, Real* error, const fd_tag<1>&) + { + using std::sqrt; + using std::pow; + using std::abs; + using std::numeric_limits; + + const Real eps = (numeric_limits<Real>::epsilon)(); + // Error bound ~eps^1/2 + // Note that this estimate of h differs from the best estimate by a factor of sqrt((|f(x)| + |f(x+h)|)/|f''(x)|). + // Since this factor is invariant under the scaling f -> kf, then we are somewhat justified in approximating it by 1. + // This approximation will get better as we move to higher orders of accuracy. + Real h = 2 * sqrt(eps); + h = detail::make_xph_representable(x, h); + + Real yh = f(x + h); + Real y0 = f(x); + Real diff = yh - y0; + if (error) + { + Real ym = f(x - h); + Real ypph = abs(yh - 2 * y0 + ym) / h; + // h*|f''(x)|*0.5 + (|f(x+h)+|f(x)|)*eps/h + *error = ypph / 2 + (abs(yh) + abs(y0))*eps / h; + } + return diff / h; + } + + template<class F, class Real> + Real finite_difference_derivative(const F f, Real x, Real* error, const fd_tag<2>&) + { + using std::sqrt; + using std::pow; + using std::abs; + using std::numeric_limits; + + const Real eps = (numeric_limits<Real>::epsilon)(); + // Error bound ~eps^2/3 + // See the previous discussion to understand determination of h and the error bound. + // Series[(f[x+h] - f[x-h])/(2*h), {h, 0, 4}] + Real h = pow(3 * eps, static_cast<Real>(1) / static_cast<Real>(3)); + h = detail::make_xph_representable(x, h); + + Real yh = f(x + h); + Real ymh = f(x - h); + Real diff = yh - ymh; + if (error) + { + Real yth = f(x + 2 * h); + Real ymth = f(x - 2 * h); + *error = eps * (abs(yh) + abs(ymh)) / (2 * h) + abs((yth - ymth) / 2 - diff) / (6 * h); + } + + return diff / (2 * h); + } + + template<class F, class Real> + Real finite_difference_derivative(const F f, Real x, Real* error, const fd_tag<4>&) + { + using std::sqrt; + using std::pow; + using std::abs; + using std::numeric_limits; + + const Real eps = (numeric_limits<Real>::epsilon)(); + // Error bound ~eps^4/5 + Real h = pow(11.25*eps, (Real)1 / (Real)5); + h = detail::make_xph_representable(x, h); + Real ymth = f(x - 2 * h); + Real yth = f(x + 2 * h); + Real yh = f(x + h); + Real ymh = f(x - h); + Real y2 = ymth - yth; + Real y1 = yh - ymh; + if (error) + { + // Mathematica code to extract the remainder: + // Series[(f[x-2*h]+ 8*f[x+h] - 8*f[x-h] - f[x+2*h])/(12*h), {h, 0, 7}] + Real y_three_h = f(x + 3 * h); + Real y_m_three_h = f(x - 3 * h); + // Error from fifth derivative: + *error = abs((y_three_h - y_m_three_h) / 2 + 2 * (ymth - yth) + 5 * (yh - ymh) / 2) / (30 * h); + // Error from function evaluation: + *error += eps * (abs(yth) + abs(ymth) + 8 * (abs(ymh) + abs(yh))) / (12 * h); + } + return (y2 + 8 * y1) / (12 * h); + } + + template<class F, class Real> + Real finite_difference_derivative(const F f, Real x, Real* error, const fd_tag<6>&) + { + using std::sqrt; + using std::pow; + using std::abs; + using std::numeric_limits; + + const Real eps = (numeric_limits<Real>::epsilon)(); + // Error bound ~eps^6/7 + // Error: h^6f^(7)(x)/140 + 5|f(x)|eps/h + Real h = pow(eps / 168, (Real)1 / (Real)7); + h = detail::make_xph_representable(x, h); + + Real yh = f(x + h); + Real ymh = f(x - h); + Real y1 = yh - ymh; + Real y2 = f(x - 2 * h) - f(x + 2 * h); + Real y3 = f(x + 3 * h) - f(x - 3 * h); + + if (error) + { + // Mathematica code to generate fd scheme for 7th derivative: + // Sum[(-1)^i*Binomial[7, i]*(f[x+(3-i)*h] + f[x+(4-i)*h])/2, {i, 0, 7}] + // Mathematica to demonstrate that this is a finite difference formula for 7th derivative: + // Series[(f[x+4*h]-f[x-4*h] + 6*(f[x-3*h] - f[x+3*h]) + 14*(f[x-h] - f[x+h] + f[x+2*h] - f[x-2*h]))/2, {h, 0, 15}] + Real y7 = (f(x + 4 * h) - f(x - 4 * h) - 6 * y3 - 14 * y1 - 14 * y2) / 2; + *error = abs(y7) / (140 * h) + 5 * (abs(yh) + abs(ymh))*eps / h; + } + return (y3 + 9 * y2 + 45 * y1) / (60 * h); + } + + template<class F, class Real> + Real finite_difference_derivative(const F f, Real x, Real* error, const fd_tag<8>&) + { + using std::sqrt; + using std::pow; + using std::abs; + using std::numeric_limits; + + const Real eps = (numeric_limits<Real>::epsilon)(); + // Error bound ~eps^8/9. + // In double precision, we only expect to lose two digits of precision while using this formula, at the cost of 8 function evaluations. + // Error: h^8|f^(9)(x)|/630 + 7|f(x)|eps/h assuming 7 unstabilized additions. + // Mathematica code to get the error: + // Series[(f[x+h]-f[x-h])*(4/5) + (1/5)*(f[x-2*h] - f[x+2*h]) + (4/105)*(f[x+3*h] - f[x-3*h]) + (1/280)*(f[x-4*h] - f[x+4*h]), {h, 0, 9}] + // If we used Kahan summation, we could get the max error down to h^8|f^(9)(x)|/630 + |f(x)|eps/h. + Real h = pow(551.25*eps, (Real)1 / (Real)9); + h = detail::make_xph_representable(x, h); + + Real yh = f(x + h); + Real ymh = f(x - h); + Real y1 = yh - ymh; + Real y2 = f(x - 2 * h) - f(x + 2 * h); + Real y3 = f(x + 3 * h) - f(x - 3 * h); + Real y4 = f(x - 4 * h) - f(x + 4 * h); + + Real tmp1 = 3 * y4 / 8 + 4 * y3; + Real tmp2 = 21 * y2 + 84 * y1; + + if (error) + { + // Mathematica code to generate fd scheme for 7th derivative: + // Sum[(-1)^i*Binomial[9, i]*(f[x+(4-i)*h] + f[x+(5-i)*h])/2, {i, 0, 9}] + // Mathematica to demonstrate that this is a finite difference formula for 7th derivative: + // Series[(f[x+5*h]-f[x- 5*h])/2 + 4*(f[x-4*h] - f[x+4*h]) + 27*(f[x+3*h] - f[x-3*h])/2 + 24*(f[x-2*h] - f[x+2*h]) + 21*(f[x+h] - f[x-h]), {h, 0, 15}] + Real f9 = (f(x + 5 * h) - f(x - 5 * h)) / 2 + 4 * y4 + 27 * y3 / 2 + 24 * y2 + 21 * y1; + *error = abs(f9) / (630 * h) + 7 * (abs(yh) + abs(ymh))*eps / h; + } + return (tmp1 + tmp2) / (105 * h); + } + + template<class F, class Real, class tag> + Real finite_difference_derivative(const F f, Real x, Real* error, const tag&) + { + // Always fails, but condition is template-arg-dependent so only evaluated if we get instantiated. + BOOST_STATIC_ASSERT_MSG(sizeof(Real) == 0, "Finite difference not implemented for this order: try 1, 2, 4, 6 or 8"); + } + +} + +template<class F, class Real, size_t order=6> +inline Real finite_difference_derivative(const F f, Real x, Real* error = nullptr) +{ + return detail::finite_difference_derivative(f, x, error, detail::fd_tag<order>()); +} + +}}} // namespaces +#endif diff --git a/boost/math/tools/polynomial_gcd.hpp b/boost/math/tools/polynomial_gcd.hpp index f8e8334869..fd0a8c2740 100644 --- a/boost/math/tools/polynomial_gcd.hpp +++ b/boost/math/tools/polynomial_gcd.hpp @@ -12,12 +12,12 @@ #endif #include <boost/math/tools/polynomial.hpp> -#include <boost/math/common_factor_rt.hpp> +#include <boost/integer/common_factor_rt.hpp> #include <boost/type_traits/is_pod.hpp> -namespace boost{ - +namespace boost{ + namespace integer { namespace gcd_detail { @@ -35,26 +35,26 @@ namespace boost{ } } - - - + + + 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. +* 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); + return x ? boost::integer::gcd_range(x.data().begin(), x.data().end()).first : T(0); } // Knuth, 4.6.1 @@ -82,13 +82,13 @@ T leading_coefficient(polynomial<T> const &x) namespace detail { - /* Reduce u and v to their primitive parts and return the gcd of their + /* 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; + using boost::integer::gcd; T const u_cont = content(u), v_cont = content(v); u /= u_cont; v /= v_cont; @@ -100,14 +100,14 @@ namespace detail /** * 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], +* +* 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. @@ -119,17 +119,17 @@ 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; @@ -154,16 +154,16 @@ subresultant_gcd(polynomial<T> u, polynomial<T> v) 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> |