bgirgis/cpp-thirdparty
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

update boost

Browse Source
This commit is contained in:
Bassem Girgis
2023-11-24 12:56:13 -06:00
parent cfc99971af
commit 19d727037a
9260 changed files with 849256 additions and 299957 deletions
Download Patch File Download Diff File Expand all files Collapse all files

131
linx64/include/boost/math/quadrature/detail/exp_sinh_detail.hpp
View File

File diff suppressed because one or more lines are too long

676
linx64/include/boost/math/quadrature/detail/ooura_fourier_integrals_detail.hpp Normal file
View File

@@ -0,0 +1,676 @@
// Copyright Nick Thompson, 2019
// 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_DETAIL_OOURA_FOURIER_INTEGRALS_DETAIL_HPP
#define BOOST_MATH_QUADRATURE_DETAIL_OOURA_FOURIER_INTEGRALS_DETAIL_HPP
#include <utility> // for std::pair.
#include <vector>
#include <iostream>
#include <boost/math/special_functions/expm1.hpp>
#include <boost/math/special_functions/sin_pi.hpp>
#include <boost/math/special_functions/cos_pi.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/tools/config.hpp>
#ifdef BOOST_HAS_THREADS
#include <mutex>
#include <atomic>
#endif
namespace boost { namespace math { namespace quadrature { namespace detail {
// Ooura and Mori, A robust double exponential formula for Fourier-type integrals,
// eta is the argument to the exponential in equation 3.3:
template<class Real>
std::pair<Real, Real> ooura_eta(Real x, Real alpha) {
using std::expm1;
using std::exp;
using std::abs;
Real expx = exp(x);
Real eta_prime = 2 + alpha/expx + expx/4;
Real eta;
// This is the fast branch:
if (abs(x) > 0.125) {
eta = 2*x - alpha*(1/expx - 1) + (expx - 1)/4;
}
else {// this is the slow branch using expm1 for small x:
eta = 2*x - alpha*expm1(-x) + expm1(x)/4;
}
return {eta, eta_prime};
}
// Ooura and Mori, A robust double exponential formula for Fourier-type integrals,
// equation 3.6:
template<class Real>
Real calculate_ooura_alpha(Real h)
{
using boost::math::constants::pi;
using std::log1p;
using std::sqrt;
Real x = sqrt(16 + 4*log1p(pi<Real>()/h)/h);
return 1/x;
}
template<class Real>
std::pair<Real, Real> ooura_sin_node_and_weight(long n, Real h, Real alpha)
{
using std::expm1;
using std::exp;
using std::abs;
using boost::math::constants::pi;
using std::isnan;
if (n == 0) {
// Equation 44 of https://arxiv.org/pdf/0911.4796.pdf
// Fourier Transform of the Stretched Exponential Function: Analytic Error Bounds,
// Double Exponential Transform, and Open-Source Implementation,
// Joachim Wuttke,
// The C library libkww provides functions to compute the Kohlrausch-Williams-Watts function,
// the Laplace-Fourier transform of the stretched (or compressed) exponential function exp(-t^beta)
// for exponent beta between 0.1 and 1.9 with sixteen decimal digits accuracy.
Real eta_prime_0 = Real(2) + alpha + Real(1)/Real(4);
Real node = pi<Real>()/(eta_prime_0*h);
Real weight = pi<Real>()*boost::math::sin_pi(1/(eta_prime_0*h));
Real eta_dbl_prime = -alpha + Real(1)/Real(4);
Real phi_prime_0 = (1 - eta_dbl_prime/(eta_prime_0*eta_prime_0))/2;
weight *= phi_prime_0;
return {node, weight};
}
Real x = n*h;
auto p = ooura_eta(x, alpha);
auto eta = p.first;
auto eta_prime = p.second;
Real expm1_meta = expm1(-eta);
Real exp_meta = exp(-eta);
Real node = -n*pi<Real>()/expm1_meta;
// I have verified that this is not a significant source of inaccuracy in the weight computation:
Real phi_prime = -(expm1_meta + x*exp_meta*eta_prime)/(expm1_meta*expm1_meta);
// The main source of inaccuracy is in computation of sin_pi.
// But I've agonized over this, and I think it's as good as it can get:
Real s = pi<Real>();
Real arg;
if(eta > 1) {
arg = n/( 1/exp_meta - 1 );
s *= boost::math::sin_pi(arg);
if (n&1) {
s *= -1;
}
}
else if (eta < -1) {
arg = n/(1-exp_meta);
s *= boost::math::sin_pi(arg);
}
else {
arg = -n*exp_meta/expm1_meta;
s *= boost::math::sin_pi(arg);
if (n&1) {
s *= -1;
}
}
Real weight = s*phi_prime;
return {node, weight};
}
#ifdef BOOST_MATH_INSTRUMENT_OOURA
template<class Real>
void print_ooura_estimate(size_t i, Real I0, Real I1, Real omega) {
using std::abs;
std::cout << std::defaultfloat
<< std::setprecision(std::numeric_limits<Real>::digits10)
<< std::fixed;
std::cout << "h = " << Real(1)/Real(1<<i) << ", I_h = " << I0/omega
<< " = " << std::hexfloat << I0/omega << ", absolute error estimate = "
<< std::defaultfloat << std::scientific << abs(I0-I1) << std::endl;
}
#endif
template<class Real>
std::pair<Real, Real> ooura_cos_node_and_weight(long n, Real h, Real alpha)
{
using std::expm1;
using std::exp;
using std::abs;
using boost::math::constants::pi;
Real x = h*(n-Real(1)/Real(2));
auto p = ooura_eta(x, alpha);
auto eta = p.first;
auto eta_prime = p.second;
Real expm1_meta = expm1(-eta);
Real exp_meta = exp(-eta);
Real node = pi<Real>()*(Real(1)/Real(2)-n)/expm1_meta;
Real phi_prime = -(expm1_meta + x*exp_meta*eta_prime)/(expm1_meta*expm1_meta);
// Takuya Ooura and Masatake Mori,
// Journal of Computational and Applied Mathematics, 112 (1999) 229-241.
// A robust double exponential formula for Fourier-type integrals.
// Equation 4.6
Real s = pi<Real>();
Real arg;
if (eta < -1) {
arg = -(n-Real(1)/Real(2))/expm1_meta;
s *= boost::math::cos_pi(arg);
}
else {
arg = -(n-Real(1)/Real(2))*exp_meta/expm1_meta;
s *= boost::math::sin_pi(arg);
if (n&1) {
s *= -1;
}
}
Real weight = s*phi_prime;
return {node, weight};
}
template<class Real>
class ooura_fourier_sin_detail {
public:
ooura_fourier_sin_detail(const Real relative_error_goal, size_t levels) {
#ifdef BOOST_MATH_INSTRUMENT_OOURA
std::cout << "ooura_fourier_sin with relative error goal " << relative_error_goal
<< " & " << levels << " levels." << std::endl;
#endif // BOOST_MATH_INSTRUMENT_OOURA
if (relative_error_goal < std::numeric_limits<Real>::epsilon() * 2) {
throw std::domain_error("The relative error goal cannot be smaller than the unit roundoff.");
}
using std::abs;
requested_levels_ = levels;
starting_level_ = 0;
rel_err_goal_ = relative_error_goal;
big_nodes_.reserve(levels);
bweights_.reserve(levels);
little_nodes_.reserve(levels);
lweights_.reserve(levels);
for (size_t i = 0; i < levels; ++i) {
if (std::is_same<Real, float>::value) {
add_level<double>(i);
}
else if (std::is_same<Real, double>::value) {
add_level<long double>(i);
}
else {
add_level<Real>(i);
}
}
}
std::vector<std::vector<Real>> const & big_nodes() const {
return big_nodes_;
}
std::vector<std::vector<Real>> const & weights_for_big_nodes() const {
return bweights_;
}
std::vector<std::vector<Real>> const & little_nodes() const {
return little_nodes_;
}
std::vector<std::vector<Real>> const & weights_for_little_nodes() const {
return lweights_;
}
template<class F>
std::pair<Real,Real> integrate(F const & f, Real omega) {
using std::abs;
using std::max;
using boost::math::constants::pi;
if (omega == 0) {
return {Real(0), Real(0)};
}
if (omega < 0) {
auto p = this->integrate(f, -omega);
return {-p.first, p.second};
}
Real I1 = std::numeric_limits<Real>::quiet_NaN();
Real relative_error_estimate = std::numeric_limits<Real>::quiet_NaN();
// As we compute integrals, we learn about their structure.
// Assuming we compute f(t)sin(wt) for many different omega, this gives some
// a posteriori ability to choose a refinement level that is roughly appropriate.
size_t i = starting_level_;
do {
Real I0 = estimate_integral(f, omega, i);
#ifdef BOOST_MATH_INSTRUMENT_OOURA
print_ooura_estimate(i, I0, I1, omega);
#endif
Real absolute_error_estimate = abs(I0-I1);
Real scale = (max)(abs(I0), abs(I1));
if (!isnan(I1) && absolute_error_estimate <= rel_err_goal_*scale) {
starting_level_ = (max)(long(i) - 1, long(0));
return {I0/omega, absolute_error_estimate/scale};
}
I1 = I0;
} while(++i < big_nodes_.size());
// We've used up all our precomputed levels.
// Now we need to add more.
// It might seems reasonable to just keep adding levels indefinitely, if that's what the user wants.
// But in fact the nodes and weights just merge into each other and the error gets worse after a certain number.
// This value for max_additional_levels was chosen by observation of a slowly converging oscillatory integral:
// f(x) := cos(7cos(x))sin(x)/x
size_t max_additional_levels = 4;
while (big_nodes_.size() < requested_levels_ + max_additional_levels) {
size_t ii = big_nodes_.size();
if (std::is_same<Real, float>::value) {
add_level<double>(ii);
}
else if (std::is_same<Real, double>::value) {
add_level<long double>(ii);
}
else {
add_level<Real>(ii);
}
Real I0 = estimate_integral(f, omega, ii);
Real absolute_error_estimate = abs(I0-I1);
Real scale = (max)(abs(I0), abs(I1));
#ifdef BOOST_MATH_INSTRUMENT_OOURA
print_ooura_estimate(ii, I0, I1, omega);
#endif
if (absolute_error_estimate <= rel_err_goal_*scale) {
starting_level_ = (max)(long(ii) - 1, long(0));
return {I0/omega, absolute_error_estimate/scale};
}
I1 = I0;
++ii;
}
starting_level_ = static_cast<long>(big_nodes_.size() - 2);
return {I1/omega, relative_error_estimate};
}
private:
template<class PreciseReal>
void add_level(size_t i) {
using std::abs;
size_t current_num_levels = big_nodes_.size();
Real unit_roundoff = std::numeric_limits<Real>::epsilon()/2;
// h0 = 1. Then all further levels have h_i = 1/2^i.
// Since the nodes don't nest, we could conceivably divide h by (say) 1.5, or 3.
// It's not clear how much benefit (or loss) would be obtained from this.
PreciseReal h = PreciseReal(1)/PreciseReal(1<<i);
std::vector<Real> bnode_row;
std::vector<Real> bweight_row;
// This is a pretty good estimate for how many elements will be placed in the vector:
bnode_row.reserve((static_cast<size_t>(1)<<i)*sizeof(Real));
bweight_row.reserve((static_cast<size_t>(1)<<i)*sizeof(Real));
std::vector<Real> lnode_row;
std::vector<Real> lweight_row;
lnode_row.reserve((static_cast<size_t>(1)<<i)*sizeof(Real));
lweight_row.reserve((static_cast<size_t>(1)<<i)*sizeof(Real));
Real max_weight = 1;
auto alpha = calculate_ooura_alpha(h);
long n = 0;
Real w;
do {
auto precise_nw = ooura_sin_node_and_weight(n, h, alpha);
Real node = static_cast<Real>(precise_nw.first);
Real weight = static_cast<Real>(precise_nw.second);
w = weight;
if (bnode_row.size() == bnode_row.capacity()) {
bnode_row.reserve(2*bnode_row.size());
bweight_row.reserve(2*bnode_row.size());
}
bnode_row.push_back(node);
bweight_row.push_back(weight);
if (abs(weight) > max_weight) {
max_weight = abs(weight);
}
++n;
// f(t)->0 as t->infty, which is why the weights are computed up to the unit roundoff.
} while(abs(w) > unit_roundoff*max_weight);
// This class tends to consume a lot of memory; shrink the vectors back down to size:
bnode_row.shrink_to_fit();
bweight_row.shrink_to_fit();
// Why we are splitting the nodes into regimes where t_n >> 1 and t_n << 1?
// It will create the opportunity to sensibly truncate the quadrature sum to significant terms.
n = -1;
do {
auto precise_nw = ooura_sin_node_and_weight(n, h, alpha);
Real node = static_cast<Real>(precise_nw.first);
if (node <= 0) {
break;
}
Real weight = static_cast<Real>(precise_nw.second);
w = weight;
using std::isnan;
if (isnan(node)) {
// This occurs at n = -11 in quad precision:
break;
}
if (lnode_row.size() > 0) {
if (lnode_row[lnode_row.size()-1] == node) {
// The nodes have fused into each other:
break;
}
}
if (lnode_row.size() == lnode_row.capacity()) {
lnode_row.reserve(2*lnode_row.size());
lweight_row.reserve(2*lnode_row.size());
}
lnode_row.push_back(node);
lweight_row.push_back(weight);
if (abs(weight) > max_weight) {
max_weight = abs(weight);
}
--n;
// f(t)->infty is possible as t->0, hence compute up to the min.
} while(abs(w) > (std::numeric_limits<Real>::min)()*max_weight);
lnode_row.shrink_to_fit();
lweight_row.shrink_to_fit();
#ifdef BOOST_HAS_THREADS
// std::scoped_lock once C++17 is more common?
std::lock_guard<std::mutex> lock(node_weight_mutex_);
#endif
// Another thread might have already finished this calculation and appended it to the nodes/weights:
if (current_num_levels == big_nodes_.size()) {
big_nodes_.push_back(bnode_row);
bweights_.push_back(bweight_row);
little_nodes_.push_back(lnode_row);
lweights_.push_back(lweight_row);
}
}
template<class F>
Real estimate_integral(F const & f, Real omega, size_t i) {
// Because so few function evaluations are required to get high accuracy on the integrals in the tests,
// Kahan summation doesn't really help.
//auto cond = boost::math::tools::summation_condition_number<Real, true>(0);
Real I0 = 0;
auto const & b_nodes = big_nodes_[i];
auto const & b_weights = bweights_[i];
// Will benchmark if this is helpful:
Real inv_omega = 1/omega;
for(size_t j = 0 ; j < b_nodes.size(); ++j) {
I0 += f(b_nodes[j]*inv_omega)*b_weights[j];
}
auto const & l_nodes = little_nodes_[i];
auto const & l_weights = lweights_[i];
// If f decays rapidly as |t|->infty, not all of these calls are necessary.
for (size_t j = 0; j < l_nodes.size(); ++j) {
I0 += f(l_nodes[j]*inv_omega)*l_weights[j];
}
return I0;
}
#ifdef BOOST_HAS_THREADS
std::mutex node_weight_mutex_;
#endif
// Nodes for n >= 0, giving t_n = pi*phi(nh)/h. Generally t_n >> 1.
std::vector<std::vector<Real>> big_nodes_;
// The term bweights_ will indicate that these are weights corresponding
// to the big nodes:
std::vector<std::vector<Real>> bweights_;
// Nodes for n < 0: Generally t_n << 1, and an invariant is that t_n > 0.
std::vector<std::vector<Real>> little_nodes_;
std::vector<std::vector<Real>> lweights_;
Real rel_err_goal_;
#ifdef BOOST_HAS_THREADS
std::atomic<long> starting_level_{};
#else
long starting_level_;
#endif
size_t requested_levels_;
};
template<class Real>
class ooura_fourier_cos_detail {
public:
ooura_fourier_cos_detail(const Real relative_error_goal, size_t levels) {
#ifdef BOOST_MATH_INSTRUMENT_OOURA
std::cout << "ooura_fourier_cos with relative error goal " << relative_error_goal
<< " & " << levels << " levels." << std::endl;
std::cout << "epsilon for type = " << std::numeric_limits<Real>::epsilon() << std::endl;
#endif // BOOST_MATH_INSTRUMENT_OOURA
if (relative_error_goal < std::numeric_limits<Real>::epsilon() * 2) {
throw std::domain_error("The relative error goal cannot be smaller than the unit roundoff!");
}
using std::abs;
requested_levels_ = levels;
starting_level_ = 0;
rel_err_goal_ = relative_error_goal;
big_nodes_.reserve(levels);
bweights_.reserve(levels);
little_nodes_.reserve(levels);
lweights_.reserve(levels);
for (size_t i = 0; i < levels; ++i) {
if (std::is_same<Real, float>::value) {
add_level<double>(i);
}
else if (std::is_same<Real, double>::value) {
add_level<long double>(i);
}
else {
add_level<Real>(i);
}
}
}
template<class F>
std::pair<Real,Real> integrate(F const & f, Real omega) {
using std::abs;
using std::max;
using boost::math::constants::pi;
if (omega == 0) {
throw std::domain_error("At omega = 0, the integral is not oscillatory. The user must choose an appropriate method for this case.\n");
}
if (omega < 0) {
return this->integrate(f, -omega);
}
Real I1 = std::numeric_limits<Real>::quiet_NaN();
Real absolute_error_estimate = std::numeric_limits<Real>::quiet_NaN();
Real scale = std::numeric_limits<Real>::quiet_NaN();
size_t i = starting_level_;
do {
Real I0 = estimate_integral(f, omega, i);
#ifdef BOOST_MATH_INSTRUMENT_OOURA
print_ooura_estimate(i, I0, I1, omega);
#endif
absolute_error_estimate = abs(I0-I1);
scale = (max)(abs(I0), abs(I1));
if (!isnan(I1) && absolute_error_estimate <= rel_err_goal_*scale) {
starting_level_ = (max)(long(i) - 1, long(0));
return {I0/omega, absolute_error_estimate/scale};
}
I1 = I0;
} while(++i < big_nodes_.size());
size_t max_additional_levels = 4;
while (big_nodes_.size() < requested_levels_ + max_additional_levels) {
size_t ii = big_nodes_.size();
if (std::is_same<Real, float>::value) {
add_level<double>(ii);
}
else if (std::is_same<Real, double>::value) {
add_level<long double>(ii);
}
else {
add_level<Real>(ii);
}
Real I0 = estimate_integral(f, omega, ii);
#ifdef BOOST_MATH_INSTRUMENT_OOURA
print_ooura_estimate(ii, I0, I1, omega);
#endif
absolute_error_estimate = abs(I0-I1);
scale = (max)(abs(I0), abs(I1));
if (absolute_error_estimate <= rel_err_goal_*scale) {
starting_level_ = (max)(long(ii) - 1, long(0));
return {I0/omega, absolute_error_estimate/scale};
}
I1 = I0;
++ii;
}
starting_level_ = static_cast<long>(big_nodes_.size() - 2);
return {I1/omega, absolute_error_estimate/scale};
}
private:
template<class PreciseReal>
void add_level(size_t i) {
using std::abs;
size_t current_num_levels = big_nodes_.size();
Real unit_roundoff = std::numeric_limits<Real>::epsilon()/2;
PreciseReal h = PreciseReal(1)/PreciseReal(1<<i);
std::vector<Real> bnode_row;
std::vector<Real> bweight_row;
bnode_row.reserve((static_cast<size_t>(1)<<i)*sizeof(Real));
bweight_row.reserve((static_cast<size_t>(1)<<i)*sizeof(Real));
std::vector<Real> lnode_row;
std::vector<Real> lweight_row;
lnode_row.reserve((static_cast<size_t>(1)<<i)*sizeof(Real));
lweight_row.reserve((static_cast<size_t>(1)<<i)*sizeof(Real));
Real max_weight = 1;
auto alpha = calculate_ooura_alpha(h);
long n = 0;
Real w;
do {
auto precise_nw = ooura_cos_node_and_weight(n, h, alpha);
Real node = static_cast<Real>(precise_nw.first);
Real weight = static_cast<Real>(precise_nw.second);
w = weight;
if (bnode_row.size() == bnode_row.capacity()) {
bnode_row.reserve(2*bnode_row.size());
bweight_row.reserve(2*bnode_row.size());
}
bnode_row.push_back(node);
bweight_row.push_back(weight);
if (abs(weight) > max_weight) {
max_weight = abs(weight);
}
++n;
// f(t)->0 as t->infty, which is why the weights are computed up to the unit roundoff.
} while(abs(w) > unit_roundoff*max_weight);
bnode_row.shrink_to_fit();
bweight_row.shrink_to_fit();
n = -1;
do {
auto precise_nw = ooura_cos_node_and_weight(n, h, alpha);
Real node = static_cast<Real>(precise_nw.first);
// The function cannot be singular at zero,
// so zero is not a unreasonable node,
// unlike in the case of the Fourier Sine.
// Hence only break if the node is negative.
if (node < 0) {
break;
}
Real weight = static_cast<Real>(precise_nw.second);
w = weight;
if (lnode_row.size() > 0) {
if (lnode_row.back() == node) {
// The nodes have fused into each other:
break;
}
}
if (lnode_row.size() == lnode_row.capacity()) {
lnode_row.reserve(2*lnode_row.size());
lweight_row.reserve(2*lnode_row.size());
}
lnode_row.push_back(node);
lweight_row.push_back(weight);
if (abs(weight) > max_weight) {
max_weight = abs(weight);
}
--n;
} while(abs(w) > (std::numeric_limits<Real>::min)()*max_weight);
lnode_row.shrink_to_fit();
lweight_row.shrink_to_fit();
#ifdef BOOST_HAS_THREADS
std::lock_guard<std::mutex> lock(node_weight_mutex_);
#endif
// Another thread might have already finished this calculation and appended it to the nodes/weights:
if (current_num_levels == big_nodes_.size()) {
big_nodes_.push_back(bnode_row);
bweights_.push_back(bweight_row);
little_nodes_.push_back(lnode_row);
lweights_.push_back(lweight_row);
}
}
template<class F>
Real estimate_integral(F const & f, Real omega, size_t i) {
Real I0 = 0;
auto const & b_nodes = big_nodes_[i];
auto const & b_weights = bweights_[i];
Real inv_omega = 1/omega;
for(size_t j = 0 ; j < b_nodes.size(); ++j) {
I0 += f(b_nodes[j]*inv_omega)*b_weights[j];
}
auto const & l_nodes = little_nodes_[i];
auto const & l_weights = lweights_[i];
for (size_t j = 0; j < l_nodes.size(); ++j) {
I0 += f(l_nodes[j]*inv_omega)*l_weights[j];
}
return I0;
}
#ifdef BOOST_HAS_THREADS
std::mutex node_weight_mutex_;
#endif
std::vector<std::vector<Real>> big_nodes_;
std::vector<std::vector<Real>> bweights_;
std::vector<std::vector<Real>> little_nodes_;
std::vector<std::vector<Real>> lweights_;
Real rel_err_goal_;
#ifdef BOOST_HAS_THREADS
std::atomic<long> starting_level_{};
#else
long starting_level_;
#endif
size_t requested_levels_;
};
}}}}
#endif

82
linx64/include/boost/math/quadrature/detail/sinh_sinh_detail.hpp
View File

File diff suppressed because one or more lines are too long

228
linx64/include/boost/math/quadrature/detail/tanh_sinh_constants.hpp
View File

@@ -1,228 +0,0 @@
// 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_QUADRATURE_DETAIL_TANH_SINH_CONSTANTS_HPP
#define BOOST_MATH_QUADRATURE_DETAIL_TANH_SINH_CONSTANTS_HPP
#include <boost/type_traits/is_constructible.hpp>
#include <boost/lexical_cast.hpp>
namespace boost {
namespace math {
namespace quadrature {
namespace detail {
template <class Real>
inline Real abscissa_at_t(const Real& t)
{
using std::tanh;
using std::sinh;
return tanh(constants::half_pi<Real>()*sinh(t));
}
template <class Real>
inline Real weight_at_t(const Real& t)
{
using std::cosh;
using std::sinh;
Real cs = cosh(constants::half_pi<Real>() * sinh(t));
return constants::half_pi<Real>() * cosh(t) / (cs * cs);
}
template <class Real>
inline Real abscissa_complement_at_t(const Real& t)
{
using std::cosh;
using std::exp;
using std::sinh;
Real u2 = constants::half_pi<Real>() * sinh(t);
return 1 / (exp(u2) *cosh(u2));
}
template <class Real>
inline Real t_from_abscissa_complement(const Real& x)
{
using std::log;
using std::sqrt;
Real l = log(sqrt((2 - x) / x));
return log((sqrt(4 * l * l + constants::pi<Real>() * constants::pi<Real>()) + 2 * l) / constants::pi<Real>());
};
template<class Real, int precision>
class tanh_sinh_detail_constants
{
static const int initializer_selector =
!std::numeric_limits<Real>::is_specialized || (std::numeric_limits<Real>::radix != 2) ?
0 :
(std::numeric_limits<Real>::digits < 30) && (std::numeric_limits<Real>::max_exponent <= 128) ?
1 : 0;
public:
tanh_sinh_detail_constants(std::size_t max_refinements, std::size_t initial_commit = 4) : m_max_refinements(max_refinements), m_committed_refinements(initial_commit)
{
typedef mpl::int_<initializer_selector> tag_type;
init(tag_type());
}
void init(const mpl::int_<0>&)
{
using std::tanh;
using std::sinh;
using std::asinh;
using std::atanh;
using boost::math::constants::half_pi;
using boost::math::constants::pi;
using boost::math::constants::two_div_pi;
m_inital_row_length = 7;
std::size_t first_complement = 0;
m_t_max = m_inital_row_length;
m_t_crossover = t_from_abscissa_complement(Real(0.5f));
m_abscissas.assign(m_max_refinements + 1, std::vector<Real>());
m_weights.assign(m_max_refinements + 1, std::vector<Real>());
m_first_complements.assign(m_max_refinements + 1, 0);
//
// First row is special:
//
Real h = m_t_max / m_inital_row_length;
std::vector<Real> temp(m_inital_row_length + 1, Real(0));
for (std::size_t i = 0; i < m_inital_row_length; ++i)
{
Real t = h * i;
if ((first_complement < 0) && (t < m_t_crossover))
first_complement = i;
temp[i] = t < m_t_crossover ? abscissa_at_t(t) : -abscissa_complement_at_t(t);
}
temp[m_inital_row_length] = abscissa_complement_at_t(m_t_max);
m_abscissas[0].swap(temp);
m_first_complements[0] = first_complement;
temp.assign(m_inital_row_length + 1, Real(0));
for (std::size_t i = 0; i < m_inital_row_length; ++i)
temp[i] = weight_at_t(Real(h * i));
temp[m_inital_row_length] = weight_at_t(m_t_max);
m_weights[0].swap(temp);
for (std::size_t row = 1; row <= m_committed_refinements; ++row)
{
h /= 2;
first_complement = 0;
for (Real pos = h; pos < m_t_max; pos += 2 * h)
{
if (pos < m_t_crossover)
++first_complement;
temp.push_back(pos < m_t_crossover ? abscissa_at_t(pos) : -abscissa_complement_at_t(pos));
}
m_abscissas[row].swap(temp);
m_first_complements[row] = first_complement;
for (Real pos = h; pos < m_t_max; pos += 2 * h)
temp.push_back(weight_at_t(pos));
m_weights[row].swap(temp);
}
}
void init(const mpl::int_<1>&)
{
m_inital_row_length = 4;
m_abscissas.reserve(m_max_refinements + 1);
m_weights.reserve(m_max_refinements + 1);
m_first_complements.reserve(m_max_refinements + 1);
m_abscissas = {
{ 0.0f, -0.04863203593f, -2.252280754e-05f, -4.294161056e-14f, -1.167648898e-37f, },
{ -0.3257285078f, -0.002485143543f, -1.112433512e-08f, -5.378491591e-23f, },
{ 0.3772097382f, -0.1404309413f, -0.01295943949f, -0.0003117359716f, -7.952628853e-07f, -4.714355182e-11f, -5.415222824e-18f, -2.040300394e-29f, },
{ 0.1943570033f, -0.4608532946f, -0.219392561f, -0.08512073674f, -0.0260331318f, -0.005944493369f, -0.0009348035442f, -9.061530486e-05f, -4.683958779e-06f, -1.072183876e-07f, -8.572949078e-10f, -1.767834693e-12f, -6.374878495e-16f, -2.442937279e-20f, -5.251546473e-26f, -2.789467162e-33f, },
{ 0.09792388529f, 0.2878799327f, 0.4612535439f, -0.3897263425f, -0.2689819652f, -0.1766829945f, -0.1101085972f, -0.06483914248f, -0.03588783578f, -0.01854517332f, -0.008873007558f, -0.003891334562f, -0.001545791232f, -0.0005485655647f, -0.0001711779271f, -4.612899437e-05f, -1.051798518e-05f, -1.982859405e-06f, -3.011058474e-07f, -3.576091908e-08f, -3.212800902e-09f, -2.102671378e-10f, -9.606066479e-12f, -2.919026641e-13f, -5.586116639e-15f, -6.328207135e-17f, -3.956461339e-19f, -1.260975845e-21f, -1.872492175e-24f, -1.170030294e-27f, -2.740986178e-31f, -2.112282721e-35f, },
{ 0.04905596731f, 0.1464179843f, 0.2415663195f, 0.3331422646f, 0.4199521113f, -0.4989866106f, -0.4244155094f, -0.356823241f, -0.2964499949f, -0.2433060914f, -0.1972012587f, -0.1577807536f, -0.1245646024f, -0.09698671849f, -0.07443136593f, -0.05626521395f, -0.04186397729f, -0.0306332671f, -0.02202376481f, -0.01554116883f, -0.01075156891f, -0.007283002803f, -0.004823973845f, -0.003119681872f, -0.001966663685f, -0.001206465701f, -0.0007188880782f, -0.0004152496485f, -0.0002320284004f, -0.0001251349512f, -6.498007492e-05f, -3.240693206e-05f, -1.548009773e-05f, -7.062123337e-06f, -3.06755081e-06f, -1.264528134e-06f, -4.929942806e-07f, -1.811062872e-07f, -6.244592163e-08f, -2.01254968e-08f, -6.035865798e-09f, -1.676638052e-09f, -4.292122274e-10f, -1.007222767e-10f, -2.154466259e-11f, -4.175393124e-12f, -7.284737264e-13f, -1.136386985e-13f, -1.57355351e-14f, -1.91921891e-15f, -2.044959541e-16f, -1.886958307e-17f, -1.493871649e-18f, -1.00469381e-19f, -5.679929178e-21f, -2.669103953e-22f, -1.030178138e-23f, -3.224484876e-25f, -8.0747829e-27f, -1.594654602e-28f, -2.445723975e-30f, -2.865919772e-32f, -2.521682525e-34f, -1.63551444e-36f, },
{ 0.02453976357f, 0.07352512299f, 0.1222291222f, 0.1704679724f, 0.2180634735f, 0.2648450766f, 0.3106517806f, 0.3553338252f, 0.3987541505f, 0.440789599f, 0.4813318461f, -0.4797119493f, -0.4424187717f, -0.4068496464f, -0.3730497919f, -0.3410490083f, -0.3108622749f, -0.2824905325f, -0.2559216165f, -0.2311313132f, -0.2080845076f, -0.1867363915f, -0.1670337061f, -0.148915992f, -0.1323168242f, -0.1171650118f, -0.1033857457f, -0.09090168184f, -0.07963394697f, -0.069503062f, -0.06042977607f, -0.05233580938f, -0.04514450419f, -0.03878138485f, -0.03317462969f, -0.02825545843f, -0.02395843974f, -0.0202217242f, -0.01698720852f, -0.01420063697f, -0.0118116462f, -0.009773759532f, -0.008044336997f, -0.006584486831f, -0.005358944287f, -0.004335923183f, -0.00348694536f, -0.002786652957f, -0.002212608041f, -0.001745083828f, -0.001366851359f, -0.001062965166f, -0.0008205510651f, -0.0006285988591f, -0.0004777623488f, -0.0003601686544f, -0.0002692384802f, -0.0001995185689f, -0.0001465272269f, -0.0001066134524f, -7.682987071e-05f, -5.481938554e-05f, -3.871519214e-05f, -2.705357477e-05f, -1.869872988e-05f, -1.2778718e-05f, -8.631551655e-06f, -5.760372383e-06f, -3.796652834e-06f, -2.470376195e-06f, -1.586189035e-06f, -1.00458931e-06f, -6.272926646e-07f, -3.860114498e-07f, -2.339766676e-07f, -1.396287854e-07f, -8.199520529e-08f, -4.735733554e-08f, -2.688676406e-08f, -1.499692369e-08f, -8.213543901e-09f, -4.414366384e-09f, -2.326763262e-09f, -1.2020165e-09f, -6.082231242e-10f, -3.012456307e-10f, -1.459438845e-10f, -6.911160499e-11f, -3.196678326e-11f, -1.443120992e-11f, -6.353676128e-12f, -2.725950519e-12f, -1.138734572e-12f, -4.627734622e-13f, -1.827990225e-13f, -7.012046738e-14f, -2.609605366e-14f, -9.413361335e-15f, -3.287925695e-15f, -1.11086294e-15f, -3.626602264e-16f, -1.142787653e-16f, -3.471902269e-17f, -1.015781907e-17f, -2.858532825e-18f, -7.727834787e-19f, -2.004423992e-19f, -4.981568376e-20f, -1.184668492e-20f, -2.691984904e-21f, -5.836667442e-22f, -1.205661589e-22f, -2.369109351e-23f, -4.421339462e-24f, -7.823835004e-25f, -1.310533144e-25f, -2.074345013e-26f, -3.096968326e-27f, -4.353200528e-28f, -5.749955668e-29f, -7.122715927e-30f, -8.25783542e-31f, -8.941541007e-32f, -9.02279665e-33f, -8.46603012e-34f, -7.369272807e-35f, -5.936632356e-36f, -4.415277839e-37f, },
{ 0.01227135512f, 0.03680228095f, 0.06129788941f, 0.08573475488f, 0.1100896299f, 0.1343395153f, 0.1584617283f, 0.1824339697f, 0.2062343883f, 0.2298416433f, 0.2532349634f, 0.2763942036f, 0.2992998981f, 0.3219333097f, 0.3442764756f, 0.3663122492f, 0.3880243378f, 0.4093973357f, 0.4304167537f, 0.4510690435f, 0.4713416183f, 0.4912228687f, -0.489297826f, -0.4702300899f, -0.4515825479f, -0.4333628262f, -0.4155775577f, -0.39823239f, -0.3813319982f, -0.3648801001f, -0.3488794756f, -0.3333319877f, -0.318238608f, -0.3035994429f, -0.2894137636f, -0.275680037f, -0.2623959593f, -0.2495584902f, -0.2371638893f, -0.2252077531f, -0.2136850525f, -0.2025901721f, -0.1919169483f, -0.1816587092f, -0.1718083133f, -0.1623581887f, -0.1533003716f, -0.1446265448f, -0.1363280753f, -0.1283960505f, -0.120821315f, -0.113594505f, -0.1067060822f, -0.1001463668f, -0.09390556883f, -0.08797381831f, -0.08234119416f, -0.07699775172f, -0.07193354881f, -0.06713867034f, -0.0626032516f, -0.0583175f, -0.0542717154f, -0.050456309f, -0.0468618209f, -0.0434789361f, -0.0402984993f, -0.03731152826f, -0.0345092259f, -0.03188299107f, -0.02942442821f, -0.02712535566f, -0.02497781299f, -0.02297406711f, -0.02110661737f, -0.01936819969f, -0.01775178968f, -0.01625060483f, -0.01485810598f, -0.01356799776f, -0.01237422849f, -0.01127098916f, -0.01025271192f, -0.009314067826f, -0.008449964034f, -0.007655540506f, -0.006926166179f, -0.006257434709f, -0.005645159804f, -0.005085370197f, -0.004574304294f, -0.004108404533f, -0.003684311494f, -0.003298857801f, -0.002949061834f, -0.002632121306f, -0.002345406714f, -0.002086454715f, -0.001852961444f, -0.001642775797f, -0.001453892723f, -0.001284446528f, -0.001132704236f, -0.0009970590063f, -0.0008760236476f, -0.0007682242321f, -0.0006723938372f, -0.000587366425f, -0.0005120708762f, -0.0004455251889f, -0.0003868308567f, -0.0003351674323f, -0.0002897872882f, -0.0002500105784f, -0.0002152204083f, -0.0001848582177f, -0.0001584193765f, -0.0001354489984f, -0.0001155379702f, -9.8319198e-05f, -8.346406605e-05f, -7.067910812e-05f, -5.970288552e-05f, -5.030306831e-05f, -4.227371392e-05f, -3.543273737e-05f, -2.961956641e-05f, -2.469297451e-05f, -2.052908425e-05f, -1.701953324e-05f, -1.406979453e-05f, -1.159764317e-05f, -9.531760786e-06f, -7.810469566e-06f, -6.380587461e-06f, -5.196396351e-06f, -4.218715072e-06f, -3.414069429e-06f, -2.753951574e-06f, -2.214161365e-06f, -1.774222689e-06f, -1.416868016e-06f, -1.127584838e-06f, -8.942180042e-07f, -7.066223315e-07f, -5.563602711e-07f, -4.364397681e-07f, -3.410878407e-07f, -2.65555767e-07f, -2.059521239e-07f, -1.591002641e-07f, -1.224171449e-07f, -9.381073151e-08f, -7.159348586e-08f, -5.440972705e-08f, -4.117489851e-08f, -3.102500976e-08f, -2.327473287e-08f, -1.738282654e-08f, -1.292373523e-08f, -9.564367214e-09f, -7.045195508e-09f, -5.164949276e-09f, -3.768272997e-09f, -2.735826254e-09f, -1.976380568e-09f, -1.420541964e-09f, -1.015790099e-09f, -7.225780068e-10f, -5.112816947e-10f, -3.598270551e-10f, -2.518536224e-10f, -1.753014775e-10f, -1.213298105e-10f, -8.349395075e-11f, -5.712266027e-11f, -3.8849686e-11f, -2.626342738e-11f, -1.764649978e-11f, -1.178329832e-11f, -7.818680628e-12f, -5.154836404e-12f, -3.376501461e-12f, -2.19707447e-12f, -1.420047367e-12f, -9.115800793e-13f, -5.811306251e-13f, -3.67867911e-13f, -2.312068904e-13f, -1.442617249e-13f, -8.934966665e-14f, -5.492567772e-14f, -3.35078831e-14f, -2.02840764e-14f, -1.218279513e-14f, -7.258853736e-15f, -4.290062787e-15f, -2.514652539e-15f, -1.461687113e-15f, -8.424330958e-16f, -4.81351759e-16f, -2.726317943e-16f, -1.530442297e-16f, -8.513785725e-17f, -4.692795039e-17f, -2.562597595e-17f, -1.386133467e-17f, -7.425750192e-18f, -3.939309402e-18f, -2.069079146e-18f, -1.075828622e-18f, -5.536671017e-19f, -2.819834004e-19f, -1.421006375e-19f, -7.084243878e-20f, -3.493350557e-20f, -1.703597751e-20f, -8.214706044e-21f, -3.915973438e-21f, -1.845149816e-21f, -8.591874581e-22f, -3.953006958e-22f, -1.7966698e-22f, -8.065408821e-23f, -3.575337212e-23f, -1.564782132e-23f, -6.760041764e-24f, -2.882136323e-24f, -1.212436233e-24f, -5.031421831e-25f, -2.059286312e-25f, -8.310787798e-26f, -3.306518068e-26f, -1.296597346e-26f, -5.010091032e-27f, -1.907179385e-27f, -7.150567334e-28f, -2.639906037e-28f, -9.594664542e-29f, -3.432078099e-29f, -1.207985066e-29f, -4.182464721e-30f, -1.424153707e-30f, -4.767833382e-31f, -1.568949178e-31f, -5.073438855e-32f, -1.611691916e-32f, -5.028356694e-33f, -1.540320437e-33f, -4.631392443e-34f, -1.366470225e-34f, -3.955008527e-35f, -1.122591825e-35f, -3.123848743e-36f, -8.519540247e-37f, -2.276473481e-37f, },
};
m_weights = {
{ 1.570796327f, 0.2300223945f, 0.0002662005138f, 1.358178427e-12f, 1.001741678e-35f, },
{ 0.9659765794f, 0.01834316699f, 2.143120456e-07f, 2.800315102e-21f, },
{ 1.389614759f, 0.5310782754f, 0.07638574357f, 0.002902517748f, 1.198370136e-05f, 1.163116581e-09f, 2.197079236e-16f, 1.363510331e-27f, },
{ 1.523283719f, 1.193463026f, 0.7374378484f, 0.3604614185f, 0.1374221077f, 0.03917500549f, 0.007742601026f, 0.0009499468043f, 6.248255924e-05f, 1.826332059e-06f, 1.868728227e-08f, 4.937853878e-11f, 2.28349267e-14f, 1.122753143e-18f, 3.09765397e-24f, 2.112123344e-31f, },
{ 1.558773356f, 1.466014427f, 1.29747575f, 1.081634985f, 0.8501728565f, 0.6304051352f, 0.4408332363f, 0.2902406793f, 0.1793244121f, 0.1034321542f, 0.05528968374f, 0.02713351001f, 0.0120835436f, 0.004816298144f, 0.001690873998f, 0.0005133938241f, 0.0001320523413f, 2.811016433e-05f, 4.823718203e-06f, 6.477756604e-07f, 6.583518513e-08f, 4.876006097e-09f, 2.521634792e-10f, 8.675931415e-12f, 1.880207173e-13f, 2.412423038e-15f, 1.708453277e-17f, 6.168256849e-20f, 1.037679724e-22f, 7.345984103e-26f, 1.949783362e-29f, 1.702438776e-33f, },
{ 1.567781431f, 1.543881116f, 1.497226223f, 1.430008355f, 1.345278885f, 1.246701207f, 1.138272243f, 1.024044933f, 0.9078793792f, 0.7932427008f, 0.6830685163f, 0.5796781031f, 0.4847580912f, 0.3993847415f, 0.3240825396f, 0.2589046395f, 0.2035239989f, 0.1573262035f, 0.1194974113f, 0.08910313924f, 0.06515553343f, 0.04666820805f, 0.03269873273f, 0.02237947106f, 0.0149378351f, 0.009707223739f, 0.006130037632f, 0.003754250977f, 0.002225082706f, 0.001273327945f, 0.0007018595157f, 0.0003716669362f, 0.0001885644298f, 9.139081749e-05f, 4.218318384e-05f, 1.84818136e-05f, 7.659575853e-06f, 2.991661588e-06f, 1.096883513e-06f, 3.759541186e-07f, 1.199244278e-07f, 3.543477717e-08f, 9.649888896e-09f, 2.409177326e-09f, 5.48283578e-10f, 1.130605535e-10f, 2.09893354e-11f, 3.484193767e-12f, 5.134127525e-13f, 6.663992283e-14f, 7.556721776e-15f, 7.420993231e-16f, 6.252804845e-17f, 4.475759507e-18f, 2.693120661e-19f, 1.346994157e-20f, 5.533583499e-22f, 1.843546975e-23f, 4.913936871e-25f, 1.032939131e-26f, 1.686277004e-28f, 2.103305749e-30f, 1.96992098e-32f, 1.359989462e-34f, },
{ 1.570042029f, 1.564021404f, 1.55205317f, 1.534281738f, 1.510919723f, 1.482243298f, 1.448586255f, 1.410332971f, 1.367910512f, 1.321780117f, 1.272428346f, 1.22035811f, 1.16607987f, 1.110103194f, 1.05292888f, 0.995041804f, 0.9369046127f, 0.8789523456f, 0.8215880353f, 0.7651792989f, 0.7100559012f, 0.6565082461f, 0.6047867306f, 0.555101878f, 0.5076251588f, 0.4624903981f, 0.4197956684f, 0.3796055694f, 0.3419537959f, 0.3068459094f, 0.2742622297f, 0.2441607779f, 0.2164802091f, 0.1911426841f, 0.1680566379f, 0.1471194133f, 0.1282197336f, 0.111239999f, 0.09605839187f, 0.08255078811f, 0.07059246991f, 0.06005964236f, 0.05083075757f, 0.04278765216f, 0.0358165056f, 0.02980862812f, 0.02466108731f, 0.02027718382f, 0.01656678625f, 0.01344653661f, 0.01083993717f, 0.00867733075f, 0.006895785969f, 0.005438899798f, 0.004256529599f, 0.003304466994f, 0.002544065768f, 0.001941835776f, 0.00146901436f, 0.001101126113f, 0.0008175410133f, 0.0006010398799f, 0.0004373949562f, 0.0003149720919f, 0.0002243596521f, 0.000158027884f, 0.0001100211285f, 7.568399659e-05f, 5.142149745e-05f, 3.449212476e-05f, 2.283211811e-05f, 1.490851403e-05f, 9.598194128e-06f, 6.089910032e-06f, 3.806198326e-06f, 2.342166721e-06f, 1.418306716e-06f, 8.447375638e-07f, 4.94582887e-07f, 2.844992366e-07f, 1.606939458e-07f, 8.907139514e-08f, 4.84209502e-08f, 2.579956823e-08f, 1.346464552e-08f, 6.878461096e-09f, 3.437185674e-09f, 1.678889768e-09f, 8.009978448e-10f, 3.729950184e-10f, 1.693945779e-10f, 7.496739757e-11f, 3.230446433e-11f, 1.354251291e-11f, 5.518236947e-12f, 2.18359221e-12f, 8.383128961e-13f, 3.119497729e-13f, 1.124020896e-13f, 3.917679451e-14f, 1.319434223e-14f, 4.289196222e-15f, 1.344322288e-15f, 4.057557702e-16f, 1.177981213e-16f, 3.285386163e-17f, 8.791316559e-18f, 2.25407483e-18f, 5.530176913e-19f, 1.296452714e-19f, 2.899964556e-20f, 6.180143249e-21f, 1.252867643e-21f, 2.412250547e-22f, 4.4039067e-23f, 7.610577808e-24f, 1.242805165e-24f, 1.91431069e-25f, 2.776125103e-26f, 3.783124073e-27f, 4.834910155e-28f, 5.783178697e-29f, 6.460575703e-30f, 6.72603739e-31f, 6.511153451e-32f, 5.847409075e-33f, 4.860046055e-34f, 3.72923953e-35f, },
{ 1.570607717f, 1.569099695f, 1.566088239f, 1.561582493f, 1.555596115f, 1.548147191f, 1.539258145f, 1.528955608f, 1.517270275f, 1.504236738f, 1.489893298f, 1.474281762f, 1.457447221f, 1.439437815f, 1.420304486f, 1.400100716f, 1.378882264f, 1.35670689f, 1.333634075f, 1.309724744f, 1.285040985f, 1.259645765f, 1.233602657f, 1.206975567f, 1.179828472f, 1.152225159f, 1.124228984f, 1.09590263f, 1.067307886f, 1.038505436f, 1.00955466f, 0.9805134517f, 0.951438051f, 0.9223828892f, 0.8934004523f, 0.8645411596f, 0.8358532563f, 0.807382723f, 0.7791731997f, 0.7512659245f, 0.723699687f, 0.6965107951f, 0.6697330554f, 0.6433977657f, 0.6175337199f, 0.5921672237f, 0.5673221206f, 0.5430198278f, 0.5192793805f, 0.4961174844f, 0.4735485755f, 0.4515848861f, 0.4302365164f, 0.4095115109f, 0.3894159397f, 0.3699539819f, 0.3511280132f, 0.3329386948f, 0.3153850641f, 0.2984646265f, 0.2821734476f, 0.2665062456f, 0.2514564831f, 0.2370164583f, 0.2231773949f, 0.2099295305f, 0.1972622032f, 0.1851639366f, 0.1736225217f, 0.1626250975f, 0.1521582278f, 0.1422079761f, 0.1327599774f, 0.1237995069f, 0.1153115463f, 0.1072808458f, 0.09969198461f, 0.09252942711f, 0.08577757654f, 0.0794208254f, 0.07344360286f, 0.06783041903f, 0.06256590638f, 0.05763485811f, 0.05302226366f, 0.04871334138f, 0.04469356846f, 0.04094870813f, 0.0374648342f, 0.03422835312f, 0.03122602351f, 0.02844497325f, 0.02587271434f, 0.02349715546f, 0.02130661237f, 0.01928981624f, 0.01743592007f, 0.01573450311f, 0.01417557353f, 0.01274956936f, 0.01144735783f, 0.01026023317f, 0.009179912924f, 0.008198533005f, 0.007308641451f, 0.006503191044f, 0.005775530877f, 0.005119396961f, 0.004528901979f, 0.003998524263f, 0.00352309611f, 0.003097791523f, 0.002718113458f, 0.002379880688f, 0.002079214354f, 0.001812524299f, 0.001576495262f, 0.00136807301f, 0.001184450486f, 0.001023054043f, 0.0008815298242f, 0.0007577303578f, 0.0006497014187f, 0.0005556692074f, 0.000474027894f, 0.0004033275645f, 0.0003422626065f, 0.0002896605611f, 0.0002444714673f, 0.0002057577147f, 0.0001726844199f, 0.0001445103343f, 0.0001205792873f, 0.0001003121646f, 8.319941724e-05f, 6.879409311e-05f, 5.670537985e-05f, 4.659264463e-05f, 3.815995412e-05f, 3.115105568e-05f, 2.534479897e-05f, 2.055097594e-05f, 1.660655598e-05f, 1.337229228e-05f, 1.072967541e-05f, 8.578209354e-06f, 6.832986277e-06f, 5.422535892e-06f, 4.286926494e-06f, 3.376095235e-06f, 2.648386225e-06f, 2.069276126e-06f, 1.610268009e-06f, 1.247935512e-06f, 9.631005212e-07f, 7.401289349e-07f, 5.66330284e-07f, 4.314482559e-07f, 3.272303733e-07f, 2.470662451e-07f, 1.856849137e-07f, 1.3890287e-07f, 1.034152804e-07f, 7.662387397e-08f, 5.649576387e-08f, 4.144823356e-08f, 3.025519646e-08f, 2.197164892e-08f, 1.587297809e-08f, 1.140646555e-08f, 8.152746483e-09f, 5.795349573e-09f, 4.096757914e-09f, 2.879701346e-09f, 2.012621022e-09f, 1.398441431e-09f, 9.659485186e-10f, 6.632086347e-10f, 4.52575761e-10f, 3.069270208e-10f, 2.068420354e-10f, 1.385028753e-10f, 9.214056423e-11f, 6.089338706e-11f, 3.997338952e-11f, 2.60619605e-11f, 1.687451934e-11f, 1.084916183e-11f, 6.925528015e-12f, 4.38886519e-12f, 2.760858767e-12f, 1.723764404e-12f, 1.068075044e-12f, 6.56694435e-13f, 4.00598538e-13f, 2.424296605e-13f, 1.455249916e-13f, 8.663812725e-14f, 5.114974901e-14f, 2.99421776e-14f, 1.737681695e-14f, 9.99642401e-15f, 5.699626666e-15f, 3.220432513e-15f, 1.802958964e-15f, 9.999957344e-16f, 5.493978397e-16f, 2.989420886e-16f, 1.610765424e-16f, 8.593209748e-17f, 4.538246827e-17f, 2.372253167e-17f, 1.227167167e-17f, 6.281229049e-18f, 3.180614714e-18f, 1.593049257e-18f, 7.890855159e-19f, 3.864733103e-19f, 1.87127733e-19f, 8.955739455e-20f, 4.235742852e-20f, 1.979436202e-20f, 9.138078558e-21f, 4.166641158e-21f, 1.876075055e-21f, 8.339901949e-22f, 3.659575236e-22f, 1.584785218e-22f, 6.771575694e-23f, 2.854281708e-23f, 1.186583858e-23f, 4.864069936e-24f, 1.965643419e-24f, 7.829165625e-25f, 3.072789229e-25f, 1.188107615e-25f, 4.524619749e-26f, 1.696710187e-26f, 6.263641003e-27f, 2.275790793e-27f, 8.136077716e-28f, 2.861306549e-28f, 9.896184197e-29f, 3.365200893e-29f, 1.124807055e-29f, 3.694460433e-30f, 1.192093301e-30f, 3.777757876e-31f, 1.175436379e-31f, 3.589879078e-32f, 1.075842686e-32f, 3.162835126e-33f, 9.118674189e-34f, 2.577393168e-34f, 7.139829504e-35f, 1.937828921e-35f, },
};
m_first_complements = {
1, 0, 1, 1, 3, 5, 11, 22,
};
while (m_abscissas.size() <= m_committed_refinements)
extend_refinements();
m_committed_refinements = m_abscissas.size() - 1;
if (m_max_refinements >= m_abscissas.size())
{
m_abscissas.resize(m_max_refinements + 1);
m_weights.resize(m_max_refinements + 1);
m_first_complements.resize(m_max_refinements + 1);
}
else
{
m_max_refinements = m_abscissas.size() - 1;
}
m_t_max = m_inital_row_length;
m_t_crossover = t_from_abscissa_complement(Real(0.5f));
}
void extend_refinements()const
{
++m_committed_refinements;
std::size_t row = m_committed_refinements;
Real h = ldexp(Real(1), -static_cast<int>(row));
std::size_t first_complement = 0;
for (Real pos = h; pos < m_t_max; pos += 2 * h)
{
if (pos < m_t_crossover)
++first_complement;
m_abscissas[row].push_back(pos < m_t_crossover ? abscissa_at_t(pos) : -abscissa_complement_at_t(pos));
}
m_first_complements[row] = first_complement;
for (Real pos = h; pos < m_t_max; pos += 2 * h)
m_weights[row].push_back(weight_at_t(pos));
}
const std::vector<Real>& get_abscissa_row(std::size_t n)const
{
if (m_committed_refinements < n)
extend_refinements();
BOOST_ASSERT(m_committed_refinements >= n);
return m_abscissas[n];
}
const std::vector<Real>& get_weight_row(std::size_t n)const
{
if (m_committed_refinements < n)
extend_refinements();
BOOST_ASSERT(m_committed_refinements >= n);
return m_weights[n];
}
std::size_t get_first_complement_index(std::size_t n)const
{
if (m_committed_refinements < n)
extend_refinements();
BOOST_ASSERT(m_committed_refinements >= n);
return m_first_complements[n];
}
protected:
mutable std::vector<std::vector<Real>> m_abscissas;
mutable std::vector<std::vector<Real>> m_weights;
mutable std::vector<std::size_t> m_first_complements;
std::size_t m_max_refinements, m_inital_row_length;
mutable std::size_t m_committed_refinements;
Real m_t_max, m_t_crossover;
};
}
}
}
}
#endif

312
linx64/include/boost/math/quadrature/detail/tanh_sinh_detail.hpp
View File

@@ -9,11 +9,15 @@
#include <cmath>
#include <vector>
#include <boost/math/tools/atomic.hpp>
#include <boost/detail/lightweight_mutex.hpp>
#include <typeinfo>
#include <boost/math/tools/atomic.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/special_functions/next.hpp>
#include <boost/math/tools/config.hpp>
#ifdef BOOST_HAS_THREADS
#include <mutex>
#endif
namespace boost{ namespace math{ namespace quadrature { namespace detail{
@@ -40,7 +44,7 @@ class tanh_sinh_detail
public:
tanh_sinh_detail(size_t max_refinements, const Real& min_complement) : m_max_refinements(max_refinements)
{
typedef mpl::int_<initializer_selector> tag_type;
typedef std::integral_constant<int, initializer_selector> tag_type;
init(min_complement, tag_type());
}
@@ -50,61 +54,61 @@ public:
private:
const std::vector<Real>& get_abscissa_row(std::size_t n)const
{
#ifndef BOOST_MATH_NO_ATOMIC_INT
#if !defined(BOOST_MATH_NO_ATOMIC_INT) && defined(BOOST_HAS_THREADS)
if (m_committed_refinements.load() < n)
extend_refinements();
BOOST_ASSERT(m_committed_refinements.load() >= n);
BOOST_MATH_ASSERT(m_committed_refinements.load() >= n);
#else
if (m_committed_refinements < n)
extend_refinements();
BOOST_ASSERT(m_committed_refinements >= n);
BOOST_MATH_ASSERT(m_committed_refinements >= n);
#endif
return m_abscissas[n];
}
const std::vector<Real>& get_weight_row(std::size_t n)const
{
#ifndef BOOST_MATH_NO_ATOMIC_INT
#if !defined(BOOST_MATH_NO_ATOMIC_INT) && defined(BOOST_HAS_THREADS)
if (m_committed_refinements.load() < n)
extend_refinements();
BOOST_ASSERT(m_committed_refinements.load() >= n);
BOOST_MATH_ASSERT(m_committed_refinements.load() >= n);
#else
if (m_committed_refinements < n)
extend_refinements();
BOOST_ASSERT(m_committed_refinements >= n);
BOOST_MATH_ASSERT(m_committed_refinements >= n);
#endif
return m_weights[n];
}
std::size_t get_first_complement_index(std::size_t n)const
{
#ifndef BOOST_MATH_NO_ATOMIC_INT
#if !defined(BOOST_MATH_NO_ATOMIC_INT) && defined(BOOST_HAS_THREADS)
if (m_committed_refinements.load() < n)
extend_refinements();
BOOST_ASSERT(m_committed_refinements.load() >= n);
BOOST_MATH_ASSERT(m_committed_refinements.load() >= n);
#else
if (m_committed_refinements < n)
extend_refinements();
BOOST_ASSERT(m_committed_refinements >= n);
BOOST_MATH_ASSERT(m_committed_refinements >= n);
#endif
return m_first_complements[n];
}
void init(const Real& min_complement, const mpl::int_<0>&);
void init(const Real& min_complement, const mpl::int_<1>&);
void init(const Real& min_complement, const mpl::int_<2>&);
void init(const Real& min_complement, const mpl::int_<3>&);
void init(const Real& min_complement, const std::integral_constant<int, 0>&);
void init(const Real& min_complement, const std::integral_constant<int, 1>&);
void init(const Real& min_complement, const std::integral_constant<int, 2>&);
void init(const Real& min_complement, const std::integral_constant<int, 3>&);
#ifdef BOOST_HAS_FLOAT128
void init(const Real& min_complement, const mpl::int_<4>&);
void init(const Real& min_complement, const std::integral_constant<int, 4>&);
#endif
void prune_to_min_complement(const Real& m);
void extend_refinements()const
{
#ifndef BOOST_MATH_NO_ATOMIC_INT
boost::detail::lightweight_mutex::scoped_lock guard(m_mutex);
#if !defined(BOOST_MATH_NO_ATOMIC_INT) && defined(BOOST_HAS_THREADS)
std::lock_guard<std::mutex> guard(m_mutex);
#endif
//
// Check some other thread hasn't got here after we read the atomic variable, but before we got here:
//
#ifndef BOOST_MATH_NO_ATOMIC_INT
#if !defined(BOOST_MATH_NO_ATOMIC_INT) && defined(BOOST_HAS_THREADS)
if (m_committed_refinements.load() >= m_max_refinements)
return;
#else
@@ -115,12 +119,12 @@ private:
using std::ldexp;
using std::ceil;
++m_committed_refinements;
#ifndef BOOST_MATH_NO_ATOMIC_INT
#if !defined(BOOST_MATH_NO_ATOMIC_INT) && defined(BOOST_HAS_THREADS)
std::size_t row = m_committed_refinements.load();
#else
std::size_t row = m_committed_refinements;
#endif
Real h = ldexp(Real(1), -static_cast<int>(row));
Real h = ldexp(static_cast<Real>(1), -static_cast<int>(row));
std::size_t first_complement = 0;
std::size_t n = boost::math::itrunc(ceil((m_t_max - h) / (2 * h)));
m_abscissas[row].reserve(n);
@@ -140,39 +144,43 @@ private:
{
using std::tanh;
using std::sinh;
return tanh(constants::half_pi<Real>()*sinh(t));
using boost::math::constants::half_pi;
return tanh(half_pi<Real>()*sinh(t));
}
static inline Real weight_at_t(const Real& t)
{
using std::cosh;
using std::sinh;
Real cs = cosh(constants::half_pi<Real>() * sinh(t));
return constants::half_pi<Real>() * cosh(t) / (cs * cs);
using boost::math::constants::half_pi;
Real cs = cosh(half_pi<Real>() * sinh(t));
return half_pi<Real>() * cosh(t) / (cs * cs);
}
static inline Real abscissa_complement_at_t(const Real& t)
{
using std::cosh;
using std::exp;
using std::sinh;
Real u2 = constants::half_pi<Real>() * sinh(t);
using boost::math::constants::half_pi;
Real u2 = half_pi<Real>() * sinh(t);
return 1 / (exp(u2) *cosh(u2));
}
static inline Real t_from_abscissa_complement(const Real& x)
{
using std::log;
using std::sqrt;
Real l = log(sqrt((2 - x) / x));
return log((sqrt(4 * l * l + constants::pi<Real>() * constants::pi<Real>()) + 2 * l) / constants::pi<Real>());
using boost::math::constants::pi;
Real l = log(2-x) - log(x);
return log((sqrt(l * l + pi<Real>() * pi<Real>()) + l) / pi<Real>());
};
mutable std::vector<std::vector<Real>> m_abscissas;
mutable std::vector<std::vector<Real>> m_weights;
mutable std::vector<std::size_t> m_first_complements;
std::size_t m_max_refinements, m_inital_row_length;
#ifndef BOOST_MATH_NO_ATOMIC_INT
mutable boost::math::detail::atomic_unsigned_type m_committed_refinements;
mutable boost::detail::lightweight_mutex m_mutex;
std::size_t m_max_refinements, m_inital_row_length{};
#if !defined(BOOST_MATH_NO_ATOMIC_INT) && defined(BOOST_HAS_THREADS)
mutable boost::math::detail::atomic_unsigned_type m_committed_refinements{};
mutable std::mutex m_mutex;
#else
mutable unsigned m_committed_refinements;
#endif
@@ -191,6 +199,14 @@ decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) tanh_sin
using std::sqrt;
using boost::math::constants::half;
using boost::math::constants::half_pi;
//
// The type of the result:
typedef decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) result_type;
Real h = m_t_max / m_inital_row_length;
result_type I0 = half_pi<Real>() * f(0, 1);
Real L1_I0 = abs(I0);
//
// We maintain 4 integer values:
// max_left_position is the logical index of the abscissa value closest to the
@@ -201,7 +217,7 @@ decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) tanh_sin
// indexes in each row, this may actually be one position to the left of max_left_position
// in the case that is even. Then, if we only evaluate f(-x_i) for abscissa values
// i <= max_left_index we will never evaluate f(-x_i) at the left endpoint.
// max_right_position and max_right_index are defained similarly for the right boundary
// max_right_position and max_right_index are defined similarly for the right boundary
// and are used to guard evaluation of f(x_i).
//
// max_left_position and max_right_position start off as the last element in row zero:
@@ -213,27 +229,69 @@ decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) tanh_sin
// of the abscissa value is greater than the smallest permitted (as specified
// by the function caller):
//
while (max_left_position && fabs(m_abscissas[0][max_left_position]) < left_min_complement)
while ((max_left_position > 1) && fabs(m_abscissas[0][max_left_position]) < left_min_complement)
--max_left_position;
while (max_right_position && fabs(m_abscissas[0][max_right_position]) < right_min_complement)
while ((max_right_position > 1) && fabs(m_abscissas[0][max_right_position]) < right_min_complement)
--max_right_position;
//
// Check for non-finite values at the end points:
//
result_type yp{ f(-1 - m_abscissas[0][max_left_position], m_abscissas[0][max_left_position]) };
result_type ym{ f(1 + m_abscissas[0][max_right_position], -m_abscissas[0][max_right_position]) };
result_type tail_tolerance{ (std::max)(boost::math::tools::epsilon<Real>(), Real(tolerance * tolerance)) };
while (max_left_position)
{
if ((boost::math::isfinite)(yp))
break;
--max_left_position;
yp = f(-1 - m_abscissas[0][max_left_position], m_abscissas[0][max_left_position]);
}
//
// Also remove points which are insignificant or zero:
//
while (max_left_position > 1)
{
if (abs(yp * m_weights[0][max_left_position]) > abs(L1_I0 * tail_tolerance))
break;
--max_left_position;
yp = f(-1 - m_abscissas[0][max_left_position], m_abscissas[0][max_left_position]);
}
//
// Over again for the right hand side:
//
while (max_right_position)
{
if ((boost::math::isfinite)(ym))
break;
--max_right_position;
ym = f(1 + m_abscissas[0][max_right_position], -m_abscissas[0][max_right_position]);
}
while (max_right_position > 1)
{
if (abs(ym * m_weights[0][max_right_position]) > abs(L1_I0 * tail_tolerance))
break;
--max_right_position;
ym = f(1 + m_abscissas[0][max_right_position], -m_abscissas[0][max_right_position]);
}
if ((max_left_position == 0) || (max_right_position == 0))
{
return policies::raise_evaluation_error(function, "The tanh_sinh quadrature found your function to be non-finite everywhere! Please check your function for singularities.", ym, Policy());
}
I0 += yp * m_weights[0][max_left_position] + ym * m_weights[0][max_right_position];
L1_I0 += abs(yp * m_weights[0][max_left_position]) + abs(ym * m_weights[0][max_right_position]);
//
// Assumption: left_min_complement/right_min_complement are sufficiently small that we only
// ever decrement through the stored values that are complements (the negative ones), and
// never ever hit the true abscissa values (positive stored values).
//
BOOST_ASSERT(m_abscissas[0][max_left_position] < 0);
BOOST_ASSERT(m_abscissas[0][max_right_position] < 0);
//
// The type of the result:
typedef decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) result_type;
BOOST_MATH_ASSERT(m_abscissas[0][max_left_position] < 0);
BOOST_MATH_ASSERT(m_abscissas[0][max_right_position] < 0);
Real h = m_t_max / m_inital_row_length;
result_type I0 = half_pi<Real>()*f(0, 1);
Real L1_I0 = abs(I0);
for(size_t i = 1; i < m_abscissas[0].size(); ++i)
{
if ((i > max_right_position) && (i > max_left_position))
if ((i >= max_right_position) && (i >= max_left_position))
break;
Real x = m_abscissas[0][i];
Real xc = x;
@@ -245,9 +303,8 @@ decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) tanh_sin
}
else
xc = x - 1;
result_type yp, ym;
yp = i <= max_right_position ? f(x, -xc) : 0;
ym = i <= max_left_position ? f(-x, xc) : 0;
yp = i < max_right_position ? f(x, -xc) : 0;
ym = i < max_left_position ? f(-x, xc) : 0;
I0 += (yp + ym)*w;
L1_I0 += (abs(yp) + abs(ym))*w;
}
@@ -264,7 +321,7 @@ decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) tanh_sin
Real err = 0;
//
// thrash_count is a heuristic - it counts how many time the error has actually increased
// rather than descreased, if this gets too high we abort...
// rather than decreased, if this gets too high we abort...
//
unsigned thrash_count = 0;
@@ -278,6 +335,7 @@ decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) tanh_sin
h *= half<Real>();
result_type sum = 0;
Real absum = 0;
Real endpoint_error = 0;
auto const& abscissa_row = this->get_abscissa_row(k);
auto const& weight_row = this->get_weight_row(k);
std::size_t first_complement_index = this->get_first_complement_index(k);
@@ -285,7 +343,7 @@ decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) tanh_sin
// At the start of each new row we need to update the max left/right indexes
// at which we can evaluate f(x_i). The new logical position is simply twice
// the old value. The new max index is one position to the left of the new
// logical value (remmember each row contains only odd numbered positions).
// logical value (remember each row contains only odd numbered positions).
// Then we have to make a single check, to see if one position to the right
// is also in bounds (this is the new abscissa value in this row which is
// known to be in between a value known to be in bounds, and one known to be
@@ -293,8 +351,10 @@ decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) tanh_sin
// Thus, we filter which abscissa values generate a call to f(x_i), with a single
// floating point comparison per loop. Everything else is integer logic.
//
BOOST_MATH_ASSERT(max_left_position);
max_left_index = max_left_position - 1;
max_left_position *= 2;
BOOST_MATH_ASSERT(max_right_position);
max_right_index = max_right_position - 1;
max_right_position *= 2;
if ((abscissa_row.size() > max_left_index + 1) && (fabs(abscissa_row[max_left_index + 1]) > left_min_complement))
@@ -307,12 +367,58 @@ decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) tanh_sin
++max_right_position;
++max_right_index;
}
//
// We also check that our endpoints don't hit singularities:
//
do
{
yp = f(-1 - abscissa_row[max_left_index], abscissa_row[max_left_index]);
if ((boost::math::isfinite)(yp))
break;
if(max_left_position <= 2)
{
return policies::raise_evaluation_error(function, "The tanh_sinh quadrature found your function to be non-finite everywhere! Please check your function for singularities.", ym, Policy());
}
max_left_position -= 2;
--max_left_index;
} while (abscissa_row[max_left_index] < 0);
bool truncate_left(false), truncate_right(false);
if (abs(L1_I1 * tail_tolerance) > abs(yp * weight_row[max_left_index]))
truncate_left = true;
do
{
ym = f(1 + abscissa_row[max_right_index], -abscissa_row[max_right_index]);
if ((boost::math::isfinite)(ym))
break;
if (max_right_position <= 2)
{
return policies::raise_evaluation_error(function, "The tanh_sinh quadrature found your function to be non-finite everywhere! Please check your function for singularities.", ym, Policy());
}
--max_right_index;
max_right_position -= 2;
} while (abscissa_row[max_right_index] < 0);
if (abs(L1_I1 * tail_tolerance) > abs(ym * weight_row[max_right_index]))
truncate_right = true;
sum += yp * weight_row[max_left_index] + ym * weight_row[max_right_index];
absum += abs(yp * weight_row[max_left_index]) + abs(ym * weight_row[max_right_index]);
//
// We estimate the error due to truncation as the value contributed by the two most extreme points.
// In most cases this is tiny and can be ignored, if it is significant then either the area of the
// integral is so far our in the tails that our exponent range can't reach it (example x^-p at double
// precision and p ~ 1), or our function is truncated near epsilon, and we have had to narrow our endpoints.
// In this latter case we may over-estimate the error, but this is the best we can do.
// In any event, we do not add endpoint_error to the error estimate until we terminate the main loop,
// otherwise it can make things appear to be non-converged, when in reality, they are as converged as they
// will ever be.
//
endpoint_error = absum;
for(size_t j = 0; j < weight_row.size(); ++j)
{
// If both left and right abscissa values are out of bounds at this step
// we can just stop this loop right now:
if ((j > max_left_index) && (j > max_right_index))
if ((j >= max_left_index) && (j >= max_right_index))
break;
Real x = abscissa_row[j];
Real xc = x;
@@ -320,17 +426,17 @@ decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) tanh_sin
if (j >= first_complement_index)
{
// We have stored x - 1:
BOOST_ASSERT(x < 0);
BOOST_MATH_ASSERT(x < 0);
x = 1 + xc;
}
else
{
BOOST_ASSERT(x >= 0);
BOOST_MATH_ASSERT(x >= 0);
xc = x - 1;
}
result_type yp = j > max_right_index ? 0 : f(x, -xc);
result_type ym = j > max_left_index ? 0 : f(-x, xc);
yp = j >= max_right_index ? 0 : f(x, -xc);
ym = j >= max_left_index ? 0 : f(-x, xc);
result_type term = (yp + ym)*w;
sum += term;
@@ -343,10 +449,10 @@ decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) tanh_sin
I1 += sum*h;
L1_I1 += absum*h;
++k;
Real last_err = err;
err = abs(I0 - I1);
// std::cout << "Estimate: " << I1 << " Error estimate at level " << k << " = " << err << std::endl;
if (!(boost::math::isfinite)(I1))
{
@@ -355,15 +461,38 @@ decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) tanh_sin
//
// If the error is increasing, and we're past level 4, something bad is very likely happening:
//
if ((err > last_err) && (k > 4) && (++thrash_count > 1))
if ((err * 1.5 > last_err) && (k > 4))
{
// We could raise an evaluation_error, but since we likely have some sort of result, just return the last one
// (ie before the error started going up)
I1 = I0;
L1_I1 = L1_I0;
--k;
err = last_err;
break;
bool terminate = false;
if ((++thrash_count > 1) && (last_err < 1e-3))
// Probably just thrashing, abort:
terminate = true;
else if(thrash_count > 2)
// OK, terrible error, but giving up anyway!
terminate = true;
else if (last_err < boost::math::tools::root_epsilon<Real>())
// Trying to squeeze precision that probably isn't there, abort:
terminate = true;
else
{
// Take a look at the end points, if there's significant new area being
// discovered, then we're not able to get close enough to the endpoints
// to ever find the integral:
if (abs(endpoint_error / sum) > err)
terminate = true;
}
if (terminate)
{
// We could raise an evaluation_error, but since we likely have some sort of result, just return the last one
// (ie before the error started going up)
I1 = I0;
L1_I1 = L1_I0;
--k;
err = last_err + endpoint_error;
break;
}
// Fall through and keep going, assume we've discovered a new feature of f(x)....
}
//
// Termination condition:
@@ -374,11 +503,41 @@ decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) tanh_sin
// parameters. We could keep hunting until we find something, but that would handicap
// integrals which really are zero.... so a compromise then!
//
if (err <= abs(tolerance*L1_I1))
if ((err <= abs(tolerance*L1_I1)) && (k >= 4))
{
//
// A quick sanity check: have we at some point narrowed our boundaries as a result
// of non-finite values? If so let's check that the area isn't on an increasing
// trajectory at our new end point, and increase our error estimate by the last
// good value as an estimate for what we may have discarded.
//
if (max_left_index && (max_left_index < abscissa_row.size() - 1) && (abs(abscissa_row[max_left_index + 1]) > left_min_complement))
{
yp = f(-1 - abscissa_row[max_left_index], abscissa_row[max_left_index]) * weight_row[max_left_index];
ym = f(-1 - abscissa_row[max_left_index - 1], abscissa_row[max_left_index - 1]) * weight_row[max_left_index - 1];
if (abs(yp) > abs(ym))
{
return policies::raise_evaluation_error(function, "The tanh_sinh quadrature evaluated your function at a singular point and got %1%. Integration bounds were automatically narrowed, but the integral was found to be increasing at the new endpoint. Please check your function, and consider providing a 2-argument functor.", I1, Policy());
}
}
if (max_right_index && (max_right_index < abscissa_row.size() - 1) && (abs(abscissa_row[max_right_index + 1]) > right_min_complement))
{
yp = f(1 + abscissa_row[max_right_index], -abscissa_row[max_right_index]) * weight_row[max_right_index];
ym = f(1 + abscissa_row[max_right_index - 1], -abscissa_row[max_right_index - 1]) * weight_row[max_right_index - 1];
if (abs(yp) > abs(ym))
{
return policies::raise_evaluation_error(function, "The tanh_sinh quadrature evaluated your function at a singular point and got %1%. Integration bounds were automatically narrowed, but the integral was found to be increasing at the new endpoint. Please check your function, and consider providing a 2-argument functor.", I1, Policy());
}
}
err += endpoint_error;
break;
}
if ((truncate_left) && (max_left_position > 1))
--max_left_position;
if ((truncate_right) && (max_right_position > 1))
--max_right_position;
}
if (error)
{
@@ -399,7 +558,7 @@ decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) tanh_sin
}
template<class Real, class Policy>
void tanh_sinh_detail<Real, Policy>::init(const Real& min_complement, const mpl::int_<0>&)
void tanh_sinh_detail<Real, Policy>::init(const Real& min_complement, const std::integral_constant<int, 0>&)
{
using std::tanh;
using std::sinh;
@@ -409,6 +568,7 @@ void tanh_sinh_detail<Real, Policy>::init(const Real& min_complement, const mpl:
using boost::math::constants::half_pi;
using boost::math::constants::pi;
using boost::math::constants::two_div_pi;
using boost::math::lltrunc;
m_committed_refinements = 4;
//
@@ -417,7 +577,7 @@ void tanh_sinh_detail<Real, Policy>::init(const Real& min_complement, const mpl:
// to ensure full precision as otherwise we chop off quite a chunk of the
// range in *subsequent* rows.
//
m_inital_row_length = itrunc(ceil(t_from_abscissa_complement(min_complement)));
m_inital_row_length = lltrunc(ceil(t_from_abscissa_complement(min_complement)));
std::size_t first_complement = 0;
m_t_max = m_inital_row_length;
m_t_crossover = t_from_abscissa_complement(Real(0.5f));
@@ -447,7 +607,7 @@ void tanh_sinh_detail<Real, Policy>::init(const Real& min_complement, const mpl:
temp[m_inital_row_length] = weight_at_t(m_t_max);
m_weights[0].swap(temp);
#ifndef BOOST_MATH_NO_ATOMIC_INT
#if !defined(BOOST_MATH_NO_ATOMIC_INT) && defined(BOOST_HAS_THREADS)
for (std::size_t row = 1; row <= m_committed_refinements.load(); ++row)
#else
for (std::size_t row = 1; row <= m_committed_refinements; ++row)
@@ -479,7 +639,7 @@ void tanh_sinh_detail<Real, Policy>::init(const Real& min_complement, const mpl:
#endif
template<class Real, class Policy>
void tanh_sinh_detail<Real, Policy>::init(const Real& min_complement, const mpl::int_<1>&)
void tanh_sinh_detail<Real, Policy>::init(const Real& min_complement, const std::integral_constant<int, 1>&)
{
m_inital_row_length = 4;
m_abscissas.reserve(m_max_refinements + 1);
@@ -508,7 +668,7 @@ void tanh_sinh_detail<Real, Policy>::init(const Real& min_complement, const mpl:
m_first_complements = {
1, 0, 1, 1, 3, 5, 11, 22,
};
#ifndef BOOST_MATH_NO_ATOMIC_INT
#if !defined(BOOST_MATH_NO_ATOMIC_INT) && defined(BOOST_HAS_THREADS)
m_committed_refinements = static_cast<boost::math::detail::atomic_unsigned_integer_type>(m_abscissas.size() - 1);
#else
m_committed_refinements = m_abscissas.size() - 1;
@@ -532,9 +692,9 @@ void tanh_sinh_detail<Real, Policy>::init(const Real& min_complement, const mpl:
}
template<class Real, class Policy>
void tanh_sinh_detail<Real, Policy>::init(const Real& min_complement, const mpl::int_<2>&)
void tanh_sinh_detail<Real, Policy>::init(const Real& min_complement, const std::integral_constant<int, 2>&)
{
m_inital_row_length = 5;
m_inital_row_length = 6;
m_abscissas.reserve(m_max_refinements + 1);
m_weights.reserve(m_max_refinements + 1);
m_first_complements.reserve(m_max_refinements + 1);
@@ -561,7 +721,7 @@ void tanh_sinh_detail<Real, Policy>::init(const Real& min_complement, const mpl:
m_first_complements = {
1, 0, 1, 1, 3, 5, 11, 22,
};
#ifndef BOOST_MATH_NO_ATOMIC_INT
#if !defined(BOOST_MATH_NO_ATOMIC_INT) && defined(BOOST_HAS_THREADS)
m_committed_refinements = static_cast<boost::math::detail::atomic_unsigned_integer_type>(m_abscissas.size() - 1);
#else
m_committed_refinements = m_abscissas.size() - 1;
@@ -584,7 +744,7 @@ void tanh_sinh_detail<Real, Policy>::init(const Real& min_complement, const mpl:
}
template<class Real, class Policy>
void tanh_sinh_detail<Real, Policy>::init(const Real& min_complement, const mpl::int_<3>&)
void tanh_sinh_detail<Real, Policy>::init(const Real& min_complement, const std::integral_constant<int, 3>&)
{
m_inital_row_length = 9;
m_abscissas.reserve(m_max_refinements + 1);
@@ -613,7 +773,7 @@ void tanh_sinh_detail<Real, Policy>::init(const Real& min_complement, const mpl:
m_first_complements = {
1, 0, 1, 1, 3, 5, 11, 22,
};
#ifndef BOOST_MATH_NO_ATOMIC_INT
#if !defined(BOOST_MATH_NO_ATOMIC_INT) && defined(BOOST_HAS_THREADS)
m_committed_refinements = static_cast<boost::math::detail::atomic_unsigned_integer_type>(m_abscissas.size() - 1);
#else
m_committed_refinements = m_abscissas.size() - 1;
@@ -638,7 +798,7 @@ void tanh_sinh_detail<Real, Policy>::init(const Real& min_complement, const mpl:
#ifdef BOOST_HAS_FLOAT128
template<class Real, class Policy>
void tanh_sinh_detail<Real, Policy>::init(const Real& min_complement, const mpl::int_<4>&)
void tanh_sinh_detail<Real, Policy>::init(const Real& min_complement, const std::integral_constant<int, 4>&)
{
m_inital_row_length = 9;
m_abscissas.reserve(m_max_refinements + 1);
@@ -667,7 +827,7 @@ void tanh_sinh_detail<Real, Policy>::init(const Real& min_complement, const mpl:
m_first_complements = {
1, 0, 1, 1, 3, 5, 11, 22,
};
#ifndef BOOST_MATH_NO_ATOMIC_INT
#if !defined(BOOST_MATH_NO_ATOMIC_INT) && defined(BOOST_HAS_THREADS)
m_committed_refinements = static_cast<boost::math::detail::atomic_unsigned_integer_type>(m_abscissas.size() - 1);
#else
m_committed_refinements = m_abscissas.size() - 1;

22
linx64/include/boost/math/quadrature/exp_sinh.hpp
View File

@@ -18,6 +18,7 @@
#include <cmath>
#include <limits>
#include <memory>
#include <string>
#include <boost/math/quadrature/detail/exp_sinh_detail.hpp>
namespace boost{ namespace math{ namespace quadrature {
@@ -30,9 +31,9 @@ public:
: m_imp(std::make_shared<detail::exp_sinh_detail<Real, Policy>>(max_refinements)) {}
template<class F>
auto integrate(const F& f, Real a, Real b, Real tol = boost::math::tools::root_epsilon<Real>(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr)->decltype(std::declval<F>()(std::declval<Real>())) const;
auto integrate(const F& f, Real a, Real b, Real tol = boost::math::tools::root_epsilon<Real>(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) const ->decltype(std::declval<F>()(std::declval<Real>()));
template<class F>
auto integrate(const F& f, Real tol = boost::math::tools::root_epsilon<Real>(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr)->decltype(std::declval<F>()(std::declval<Real>())) const;
auto integrate(const F& f, Real tol = boost::math::tools::root_epsilon<Real>(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) const ->decltype(std::declval<F>()(std::declval<Real>()));
private:
std::shared_ptr<detail::exp_sinh_detail<Real, Policy>> m_imp;
@@ -40,9 +41,11 @@ private:
template<class Real, class Policy>
template<class F>
auto exp_sinh<Real, Policy>::integrate(const F& f, Real a, Real b, Real tolerance, Real* error, Real* L1, std::size_t* levels)->decltype(std::declval<F>()(std::declval<Real>())) const
auto exp_sinh<Real, Policy>::integrate(const F& f, Real a, Real b, Real tolerance, Real* error, Real* L1, std::size_t* levels) const ->decltype(std::declval<F>()(std::declval<Real>()))
{
typedef decltype(f(a)) K;
static_assert(!std::is_integral<K>::value,
"The return type cannot be integral, it must be either a real or complex floating point type.");
using std::abs;
using boost::math::constants::half;
using boost::math::quadrature::detail::exp_sinh_detail;
@@ -58,17 +61,17 @@ auto exp_sinh<Real, Policy>::integrate(const F& f, Real a, Real b, Real toleranc
if ((boost::math::isfinite)(a) && (b >= boost::math::tools::max_value<Real>()))
{
// If a = 0, don't use an additional level of indirection:
if (a == (Real) 0)
if (a == static_cast<Real>(0))
{
return m_imp->integrate(f, error, L1, function, tolerance, levels);
}
const auto u = [&](Real t)->Real { return f(t + a); };
const auto u = [&](Real t)->K { return f(t + a); };
return m_imp->integrate(u, error, L1, function, tolerance, levels);
}
if ((boost::math::isfinite)(b) && a <= -boost::math::tools::max_value<Real>())
{
const auto u = [&](Real t)->Real { return f(b-t);};
const auto u = [&](Real t)->K { return f(b-t);};
return m_imp->integrate(u, error, L1, function, tolerance, levels);
}
@@ -83,9 +86,14 @@ auto exp_sinh<Real, Policy>::integrate(const F& f, Real a, Real b, Real toleranc
template<class Real, class Policy>
template<class F>
auto exp_sinh<Real, Policy>::integrate(const F& f, Real tolerance, Real* error, Real* L1, std::size_t* levels)->decltype(std::declval<F>()(std::declval<Real>())) const
auto exp_sinh<Real, Policy>::integrate(const F& f, Real tolerance, Real* error, Real* L1, std::size_t* levels) const ->decltype(std::declval<F>()(std::declval<Real>()))
{
static const char* function = "boost::math::quadrature::exp_sinh<%1%>::integrate";
using std::abs;
if (abs(tolerance) > 1) {
std::string msg = std::string(__FILE__) + ":" + std::to_string(__LINE__) + ":" + std::string(function) + ": The tolerance provided is unusually large; did you confuse it with a domain bound?";
throw std::domain_error(msg);
}
return m_imp->integrate(f, error, L1, function, tolerance, levels);
}

91
linx64/include/boost/math/quadrature/gauss.hpp
View File

@@ -10,6 +10,7 @@
#pragma once
#endif
#include <limits>
#include <vector>
#include <boost/math/special_functions/legendre.hpp>
#include <boost/math/constants/constants.hpp>
@@ -28,11 +29,11 @@ struct gauss_constant_category
(std::numeric_limits<T>::is_specialized == 0) ? 999 :
(std::numeric_limits<T>::radix == 2) ?
(
(std::numeric_limits<T>::digits <= std::numeric_limits<float>::digits) && boost::is_convertible<float, T>::value ? 0 :
(std::numeric_limits<T>::digits <= std::numeric_limits<double>::digits) && boost::is_convertible<double, T>::value ? 1 :
(std::numeric_limits<T>::digits <= std::numeric_limits<long double>::digits) && boost::is_convertible<long double, T>::value ? 2 :
(std::numeric_limits<T>::digits <= std::numeric_limits<float>::digits) && std::is_convertible<float, T>::value ? 0 :
(std::numeric_limits<T>::digits <= std::numeric_limits<double>::digits) && std::is_convertible<double, T>::value ? 1 :
(std::numeric_limits<T>::digits <= std::numeric_limits<long double>::digits) && std::is_convertible<long double, T>::value ? 2 :
#ifdef BOOST_HAS_FLOAT128
(std::numeric_limits<T>::digits <= 113) && boost::is_constructible<__float128, T>::value ? 3 :
(std::numeric_limits<T>::digits <= 113) && std::is_constructible<__float128, T>::value ? 3 :
#endif
(std::numeric_limits<T>::digits10 <= 110) ? 4 : 999
) : (std::numeric_limits<T>::digits10 <= 110) ? 4 : 999;
@@ -80,7 +81,7 @@ class gauss_detail<T, 7, 0>
public:
static std::array<T, 4> const & abscissa()
{
static const std::array<T, 4> data = {
static constexpr std::array<T, 4> data = {
0.000000000e+00f,
4.058451514e-01f,
7.415311856e-01f,
@@ -90,7 +91,7 @@ public:
}
static std::array<T, 4> const & weights()
{
static const std::array<T, 4> data = {
static constexpr std::array<T, 4> data = {
4.179591837e-01f,
3.818300505e-01f,
2.797053915e-01f,
@@ -106,7 +107,7 @@ class gauss_detail<T, 7, 1>
public:
static std::array<T, 4> const & abscissa()
{
static const std::array<T, 4> data = {
static constexpr std::array<T, 4> data = {
0.00000000000000000e+00,
4.05845151377397167e-01,
7.41531185599394440e-01,
@@ -116,7 +117,7 @@ public:
}
static std::array<T, 4> const & weights()
{
static const std::array<T, 4> data = {
static constexpr std::array<T, 4> data = {
4.17959183673469388e-01,
3.81830050505118945e-01,
2.79705391489276668e-01,
@@ -132,7 +133,7 @@ class gauss_detail<T, 7, 2>
public:
static std::array<T, 4> const & abscissa()
{
static const std::array<T, 4> data = {
static constexpr std::array<T, 4> data = {
0.00000000000000000000000000000000000e+00L,
4.05845151377397166906606412076961463e-01L,
7.41531185599394439863864773280788407e-01L,
@@ -142,7 +143,7 @@ public:
}
static std::array<T, 4> const & weights()
{
static const std::array<T, 4> data = {
static constexpr std::array<T, 4> data = {
4.17959183673469387755102040816326531e-01L,
3.81830050505118944950369775488975134e-01L,
2.79705391489276667901467771423779582e-01L,
@@ -210,7 +211,7 @@ class gauss_detail<T, 10, 0>
public:
static std::array<T, 5> const & abscissa()
{
static const std::array<T, 5> data = {
static constexpr std::array<T, 5> data = {
1.488743390e-01f,
4.333953941e-01f,
6.794095683e-01f,
@@ -221,7 +222,7 @@ public:
}
static std::array<T, 5> const & weights()
{
static const std::array<T, 5> data = {
static constexpr std::array<T, 5> data = {
2.955242247e-01f,
2.692667193e-01f,
2.190863625e-01f,
@@ -238,7 +239,7 @@ class gauss_detail<T, 10, 1>
public:
static std::array<T, 5> const & abscissa()
{
static const std::array<T, 5> data = {
static constexpr std::array<T, 5> data = {
1.48874338981631211e-01,
4.33395394129247191e-01,
6.79409568299024406e-01,
@@ -249,7 +250,7 @@ public:
}
static std::array<T, 5> const & weights()
{
static const std::array<T, 5> data = {
static constexpr std::array<T, 5> data = {
2.95524224714752870e-01,
2.69266719309996355e-01,
2.19086362515982044e-01,
@@ -266,7 +267,7 @@ class gauss_detail<T, 10, 2>
public:
static std::array<T, 5> const & abscissa()
{
static const std::array<T, 5> data = {
static constexpr std::array<T, 5> data = {
1.48874338981631210884826001129719985e-01L,
4.33395394129247190799265943165784162e-01L,
6.79409568299024406234327365114873576e-01L,
@@ -277,7 +278,7 @@ public:
}
static std::array<T, 5> const & weights()
{
static const std::array<T, 5> data = {
static constexpr std::array<T, 5> data = {
2.95524224714752870173892994651338329e-01L,
2.69266719309996355091226921569469353e-01L,
2.19086362515982043995534934228163192e-01L,
@@ -350,7 +351,7 @@ class gauss_detail<T, 15, 0>
public:
static std::array<T, 8> const & abscissa()
{
static const std::array<T, 8> data = {
static constexpr std::array<T, 8> data = {
0.000000000e+00f,
2.011940940e-01f,
3.941513471e-01f,
@@ -364,7 +365,7 @@ public:
}
static std::array<T, 8> const & weights()
{
static const std::array<T, 8> data = {
static constexpr std::array<T, 8> data = {
2.025782419e-01f,
1.984314853e-01f,
1.861610000e-01f,
@@ -384,7 +385,7 @@ class gauss_detail<T, 15, 1>
public:
static std::array<T, 8> const & abscissa()
{
static const std::array<T, 8> data = {
static constexpr std::array<T, 8> data = {
0.00000000000000000e+00,
2.01194093997434522e-01,
3.94151347077563370e-01,
@@ -398,7 +399,7 @@ public:
}
static std::array<T, 8> const & weights()
{
static const std::array<T, 8> data = {
static constexpr std::array<T, 8> data = {
2.02578241925561273e-01,
1.98431485327111576e-01,
1.86161000015562211e-01,
@@ -418,7 +419,7 @@ class gauss_detail<T, 15, 2>
public:
static std::array<T, 8> const & abscissa()
{
static const std::array<T, 8> data = {
static constexpr std::array<T, 8> data = {
0.00000000000000000000000000000000000e+00L,
2.01194093997434522300628303394596208e-01L,
3.94151347077563369897207370981045468e-01L,
@@ -432,7 +433,7 @@ public:
}
static std::array<T, 8> const & weights()
{
static const std::array<T, 8> data = {
static constexpr std::array<T, 8> data = {
2.02578241925561272880620199967519315e-01L,
1.98431485327111576456118326443839325e-01L,
1.86161000015562211026800561866422825e-01L,
@@ -520,7 +521,7 @@ class gauss_detail<T, 20, 0>
public:
static std::array<T, 10> const & abscissa()
{
static const std::array<T, 10> data = {
static constexpr std::array<T, 10> data = {
7.652652113e-02f,
2.277858511e-01f,
3.737060887e-01f,
@@ -536,7 +537,7 @@ public:
}
static std::array<T, 10> const & weights()
{
static const std::array<T, 10> data = {
static constexpr std::array<T, 10> data = {
1.527533871e-01f,
1.491729865e-01f,
1.420961093e-01f,
@@ -558,7 +559,7 @@ class gauss_detail<T, 20, 1>
public:
static std::array<T, 10> const & abscissa()
{
static const std::array<T, 10> data = {
static constexpr std::array<T, 10> data = {
7.65265211334973338e-02,
2.27785851141645078e-01,
3.73706088715419561e-01,
@@ -574,7 +575,7 @@ public:
}
static std::array<T, 10> const & weights()
{
static const std::array<T, 10> data = {
static constexpr std::array<T, 10> data = {
1.52753387130725851e-01,
1.49172986472603747e-01,
1.42096109318382051e-01,
@@ -596,7 +597,7 @@ class gauss_detail<T, 20, 2>
public:
static std::array<T, 10> const & abscissa()
{
static const std::array<T, 10> data = {
static constexpr std::array<T, 10> data = {
7.65265211334973337546404093988382110e-02L,
2.27785851141645078080496195368574625e-01L,
3.73706088715419560672548177024927237e-01L,
@@ -612,7 +613,7 @@ public:
}
static std::array<T, 10> const & weights()
{
static const std::array<T, 10> data = {
static constexpr std::array<T, 10> data = {
1.52753387130725850698084331955097593e-01L,
1.49172986472603746787828737001969437e-01L,
1.42096109318382051329298325067164933e-01L,
@@ -710,7 +711,7 @@ class gauss_detail<T, 25, 0>
public:
static std::array<T, 13> const & abscissa()
{
static const std::array<T, 13> data = {
static constexpr std::array<T, 13> data = {
0.000000000e+00f,
1.228646926e-01f,
2.438668837e-01f,
@@ -729,7 +730,7 @@ public:
}
static std::array<T, 13> const & weights()
{
static const std::array<T, 13> data = {
static constexpr std::array<T, 13> data = {
1.231760537e-01f,
1.222424430e-01f,
1.194557635e-01f,
@@ -754,7 +755,7 @@ class gauss_detail<T, 25, 1>
public:
static std::array<T, 13> const & abscissa()
{
static const std::array<T, 13> data = {
static constexpr std::array<T, 13> data = {
0.00000000000000000e+00,
1.22864692610710396e-01,
2.43866883720988432e-01,
@@ -773,7 +774,7 @@ public:
}
static std::array<T, 13> const & weights()
{
static const std::array<T, 13> data = {
static constexpr std::array<T, 13> data = {
1.23176053726715451e-01,
1.22242442990310042e-01,
1.19455763535784772e-01,
@@ -798,7 +799,7 @@ class gauss_detail<T, 25, 2>
public:
static std::array<T, 13> const & abscissa()
{
static const std::array<T, 13> data = {
static constexpr std::array<T, 13> data = {
0.00000000000000000000000000000000000e+00L,
1.22864692610710396387359818808036806e-01L,
2.43866883720988432045190362797451586e-01L,
@@ -817,7 +818,7 @@ public:
}
static std::array<T, 13> const & weights()
{
static const std::array<T, 13> data = {
static constexpr std::array<T, 13> data = {
1.23176053726715451203902873079050142e-01L,
1.22242442990310041688959518945851506e-01L,
1.19455763535784772228178126512901047e-01L,
@@ -930,7 +931,7 @@ class gauss_detail<T, 30, 0>
public:
static std::array<T, 15> const & abscissa()
{
static const std::array<T, 15> data = {
static constexpr std::array<T, 15> data = {
5.147184256e-02f,
1.538699136e-01f,
2.546369262e-01f,
@@ -951,7 +952,7 @@ public:
}
static std::array<T, 15> const & weights()
{
static const std::array<T, 15> data = {
static constexpr std::array<T, 15> data = {
1.028526529e-01f,
1.017623897e-01f,
9.959342059e-02f,
@@ -978,7 +979,7 @@ class gauss_detail<T, 30, 1>
public:
static std::array<T, 15> const & abscissa()
{
static const std::array<T, 15> data = {
static constexpr std::array<T, 15> data = {
5.14718425553176958e-02,
1.53869913608583547e-01,
2.54636926167889846e-01,
@@ -999,7 +1000,7 @@ public:
}
static std::array<T, 15> const & weights()
{
static const std::array<T, 15> data = {
static constexpr std::array<T, 15> data = {
1.02852652893558840e-01,
1.01762389748405505e-01,
9.95934205867952671e-02,
@@ -1026,7 +1027,7 @@ class gauss_detail<T, 30, 2>
public:
static std::array<T, 15> const & abscissa()
{
static const std::array<T, 15> data = {
static constexpr std::array<T, 15> data = {
5.14718425553176958330252131667225737e-02L,
1.53869913608583546963794672743255920e-01L,
2.54636926167889846439805129817805108e-01L,
@@ -1047,7 +1048,7 @@ public:
}
static std::array<T, 15> const & weights()
{
static const std::array<T, 15> data = {
static constexpr std::array<T, 15> data = {
1.02852652893558840341285636705415044e-01L,
1.01762389748405504596428952168554045e-01L,
9.95934205867952670627802821035694765e-02L,
@@ -1178,6 +1179,8 @@ public:
// In many math texts, K represents the field of real or complex numbers.
// Too bad we can't put blackboard bold into C++ source!
typedef decltype(f(Real(0))) K;
static_assert(!std::is_integral<K>::value,
"The return type cannot be integral, it must be either a real or complex floating point type.");
using std::abs;
unsigned non_zero_start = 1;
K result = Real(0);
@@ -1258,9 +1261,13 @@ public:
if ((boost::math::isfinite)(a) && (boost::math::isfinite)(b))
{
if (b <= a)
if (a == b)
{
return policies::raise_domain_error(function, "Arguments to integrate are in wrong order; integration over [a,b] must have b > a.", a, Policy());
return K(0);
}
if (b < a)
{
return -integrate(f, b, a, pL1);
}
Real avg = (a + b)*constants::half<Real>();
Real scale = (b - a)*constants::half<Real>();

86
linx64/include/boost/math/quadrature/gauss_kronrod.hpp
View File

@@ -1,5 +1,5 @@
// Copyright John Maddock 2017.
// Copyright Nick Thompson 2017.
// 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)
@@ -13,7 +13,9 @@
#pragma warning(disable: 4127)
#endif
#include <array>
#include <vector>
#include <algorithm>
#include <boost/math/special_functions/legendre.hpp>
#include <boost/math/special_functions/legendre_stieltjes.hpp>
#include <boost/math/quadrature/gauss.hpp>
@@ -86,7 +88,7 @@ class gauss_kronrod_detail<T, 15, 0>
public:
static std::array<T, 8> const & abscissa()
{
static const std::array<T, 8> data = {
static constexpr std::array<T, 8> data = {
0.000000000e+00f,
2.077849550e-01f,
4.058451514e-01f,
@@ -100,7 +102,7 @@ public:
}
static std::array<T, 8> const & weights()
{
static const std::array<T, 8> data = {
static constexpr std::array<T, 8> data = {
2.094821411e-01f,
2.044329401e-01f,
1.903505781e-01f,
@@ -120,7 +122,7 @@ class gauss_kronrod_detail<T, 15, 1>
public:
static std::array<T, 8> const & abscissa()
{
static const std::array<T, 8> data = {
static constexpr std::array<T, 8> data = {
0.00000000000000000e+00,
2.07784955007898468e-01,
4.05845151377397167e-01,
@@ -134,7 +136,7 @@ public:
}
static std::array<T, 8> const & weights()
{
static const std::array<T, 8> data = {
static constexpr std::array<T, 8> data = {
2.09482141084727828e-01,
2.04432940075298892e-01,
1.90350578064785410e-01,
@@ -154,7 +156,7 @@ class gauss_kronrod_detail<T, 15, 2>
public:
static std::array<T, 8> const & abscissa()
{
static const std::array<T, 8> data = {
static constexpr std::array<T, 8> data = {
0.00000000000000000000000000000000000e+00L,
2.07784955007898467600689403773244913e-01L,
4.05845151377397166906606412076961463e-01L,
@@ -168,7 +170,7 @@ public:
}
static std::array<T, 8> const & weights()
{
static const std::array<T, 8> data = {
static constexpr std::array<T, 8> data = {
2.09482141084727828012999174891714264e-01L,
2.04432940075298892414161999234649085e-01L,
1.90350578064785409913256402421013683e-01L,
@@ -258,7 +260,7 @@ class gauss_kronrod_detail<T, 21, 0>
public:
static std::array<T, 11> const & abscissa()
{
static const std::array<T, 11> data = {
static constexpr std::array<T, 11> data = {
0.000000000e+00f,
1.488743390e-01f,
2.943928627e-01f,
@@ -275,7 +277,7 @@ public:
}
static std::array<T, 11> const & weights()
{
static const std::array<T, 11> data = {
static constexpr std::array<T, 11> data = {
1.494455540e-01f,
1.477391049e-01f,
1.427759386e-01f,
@@ -298,7 +300,7 @@ class gauss_kronrod_detail<T, 21, 1>
public:
static std::array<T, 11> const & abscissa()
{
static const std::array<T, 11> data = {
static constexpr std::array<T, 11> data = {
0.00000000000000000e+00,
1.48874338981631211e-01,
2.94392862701460198e-01,
@@ -315,7 +317,7 @@ public:
}
static std::array<T, 11> const & weights()
{
static const std::array<T, 11> data = {
static constexpr std::array<T, 11> data = {
1.49445554002916906e-01,
1.47739104901338491e-01,
1.42775938577060081e-01,
@@ -338,7 +340,7 @@ class gauss_kronrod_detail<T, 21, 2>
public:
static std::array<T, 11> const & abscissa()
{
static const std::array<T, 11> data = {
static constexpr std::array<T, 11> data = {
0.00000000000000000000000000000000000e+00L,
1.48874338981631210884826001129719985e-01L,
2.94392862701460198131126603103865566e-01L,
@@ -355,7 +357,7 @@ public:
}
static std::array<T, 11> const & weights()
{
static const std::array<T, 11> data = {
static constexpr std::array<T, 11> data = {
1.49445554002916905664936468389821204e-01L,
1.47739104901338491374841515972068046e-01L,
1.42775938577060080797094273138717061e-01L,
@@ -460,7 +462,7 @@ class gauss_kronrod_detail<T, 31, 0>
public:
static std::array<T, 16> const & abscissa()
{
static const std::array<T, 16> data = {
static constexpr std::array<T, 16> data = {
0.000000000e+00f,
1.011420669e-01f,
2.011940940e-01f,
@@ -482,7 +484,7 @@ public:
}
static std::array<T, 16> const & weights()
{
static const std::array<T, 16> data = {
static constexpr std::array<T, 16> data = {
1.013300070e-01f,
1.007698455e-01f,
9.917359872e-02f,
@@ -510,7 +512,7 @@ class gauss_kronrod_detail<T, 31, 1>
public:
static std::array<T, 16> const & abscissa()
{
static const std::array<T, 16> data = {
static constexpr std::array<T, 16> data = {
0.00000000000000000e+00,
1.01142066918717499e-01,
2.01194093997434522e-01,
@@ -532,7 +534,7 @@ public:
}
static std::array<T, 16> const & weights()
{
static const std::array<T, 16> data = {
static constexpr std::array<T, 16> data = {
1.01330007014791549e-01,
1.00769845523875595e-01,
9.91735987217919593e-02,
@@ -560,7 +562,7 @@ class gauss_kronrod_detail<T, 31, 2>
public:
static std::array<T, 16> const & abscissa()
{
static const std::array<T, 16> data = {
static constexpr std::array<T, 16> data = {
0.00000000000000000000000000000000000e+00L,
1.01142066918717499027074231447392339e-01L,
2.01194093997434522300628303394596208e-01L,
@@ -582,7 +584,7 @@ public:
}
static std::array<T, 16> const & weights()
{
static const std::array<T, 16> data = {
static constexpr std::array<T, 16> data = {
1.01330007014791549017374792767492547e-01L,
1.00769845523875595044946662617569722e-01L,
9.91735987217919593323931734846031311e-02L,
@@ -712,7 +714,7 @@ class gauss_kronrod_detail<T, 41, 0>
public:
static std::array<T, 21> const & abscissa()
{
static const std::array<T, 21> data = {
static constexpr std::array<T, 21> data = {
0.000000000e+00f,
7.652652113e-02f,
1.526054652e-01f,
@@ -739,7 +741,7 @@ public:
}
static std::array<T, 21> const & weights()
{
static const std::array<T, 21> data = {
static constexpr std::array<T, 21> data = {
7.660071192e-02f,
7.637786767e-02f,
7.570449768e-02f,
@@ -772,7 +774,7 @@ class gauss_kronrod_detail<T, 41, 1>
public:
static std::array<T, 21> const & abscissa()
{
static const std::array<T, 21> data = {
static constexpr std::array<T, 21> data = {
0.00000000000000000e+00,
7.65265211334973338e-02,
1.52605465240922676e-01,
@@ -799,7 +801,7 @@ public:
}
static std::array<T, 21> const & weights()
{
static const std::array<T, 21> data = {
static constexpr std::array<T, 21> data = {
7.66007119179996564e-02,
7.63778676720807367e-02,
7.57044976845566747e-02,
@@ -832,7 +834,7 @@ class gauss_kronrod_detail<T, 41, 2>
public:
static std::array<T, 21> const & abscissa()
{
static const std::array<T, 21> data = {
static constexpr std::array<T, 21> data = {
0.00000000000000000000000000000000000e+00L,
7.65265211334973337546404093988382110e-02L,
1.52605465240922675505220241022677528e-01L,
@@ -859,7 +861,7 @@ public:
}
static std::array<T, 21> const & weights()
{
static const std::array<T, 21> data = {
static constexpr std::array<T, 21> data = {
7.66007119179996564450499015301017408e-02L,
7.63778676720807367055028350380610018e-02L,
7.57044976845566746595427753766165583e-02L,
@@ -1014,7 +1016,7 @@ class gauss_kronrod_detail<T, 51, 0>
public:
static std::array<T, 26> const & abscissa()
{
static const std::array<T, 26> data = {
static constexpr std::array<T, 26> data = {
0.000000000e+00f,
6.154448301e-02f,
1.228646926e-01f,
@@ -1046,7 +1048,7 @@ public:
}
static std::array<T, 26> const & weights()
{
static const std::array<T, 26> data = {
static constexpr std::array<T, 26> data = {
6.158081807e-02f,
6.147118987e-02f,
6.112850972e-02f,
@@ -1084,7 +1086,7 @@ class gauss_kronrod_detail<T, 51, 1>
public:
static std::array<T, 26> const & abscissa()
{
static const std::array<T, 26> data = {
static constexpr std::array<T, 26> data = {
0.00000000000000000e+00,
6.15444830056850789e-02,
1.22864692610710396e-01,
@@ -1116,7 +1118,7 @@ public:
}
static std::array<T, 26> const & weights()
{
static const std::array<T, 26> data = {
static constexpr std::array<T, 26> data = {
6.15808180678329351e-02,
6.14711898714253167e-02,
6.11285097170530483e-02,
@@ -1154,7 +1156,7 @@ class gauss_kronrod_detail<T, 51, 2>
public:
static std::array<T, 26> const & abscissa()
{
static const std::array<T, 26> data = {
static constexpr std::array<T, 26> data = {
0.00000000000000000000000000000000000e+00L,
6.15444830056850788865463923667966313e-02L,
1.22864692610710396387359818808036806e-01L,
@@ -1186,7 +1188,7 @@ public:
}
static std::array<T, 26> const & weights()
{
static const std::array<T, 26> data = {
static constexpr std::array<T, 26> data = {
6.15808180678329350787598242400645532e-02L,
6.14711898714253166615441319652641776e-02L,
6.11285097170530483058590304162927119e-02L,
@@ -1366,7 +1368,7 @@ class gauss_kronrod_detail<T, 61, 0>
public:
static std::array<T, 31> const & abscissa()
{
static const std::array<T, 31> data = {
static constexpr std::array<T, 31> data = {
0.000000000e+00f,
5.147184256e-02f,
1.028069380e-01f,
@@ -1403,7 +1405,7 @@ public:
}
static std::array<T, 31> const & weights()
{
static const std::array<T, 31> data = {
static constexpr std::array<T, 31> data = {
5.149472943e-02f,
5.142612854e-02f,
5.122154785e-02f,
@@ -1446,7 +1448,7 @@ class gauss_kronrod_detail<T, 61, 1>
public:
static std::array<T, 31> const & abscissa()
{
static const std::array<T, 31> data = {
static constexpr std::array<T, 31> data = {
0.00000000000000000e+00,
5.14718425553176958e-02,
1.02806937966737030e-01,
@@ -1483,7 +1485,7 @@ public:
}
static std::array<T, 31> const & weights()
{
static const std::array<T, 31> data = {
static constexpr std::array<T, 31> data = {
5.14947294294515676e-02,
5.14261285374590259e-02,
5.12215478492587722e-02,
@@ -1526,7 +1528,7 @@ class gauss_kronrod_detail<T, 61, 2>
public:
static std::array<T, 31> const & abscissa()
{
static const std::array<T, 31> data = {
static constexpr std::array<T, 31> data = {
0.00000000000000000000000000000000000e+00L,
5.14718425553176958330252131667225737e-02L,
1.02806937966737030147096751318000592e-01L,
@@ -1563,7 +1565,7 @@ public:
}
static std::array<T, 31> const & weights()
{
static const std::array<T, 31> data = {
static constexpr std::array<T, 31> data = {
5.14947294294515675583404336470993075e-02L,
5.14261285374590259338628792157812598e-02L,
5.12215478492587721706562826049442083e-02L,
@@ -1869,6 +1871,8 @@ public:
static auto integrate(F f, Real a, Real b, unsigned max_depth = 15, Real tol = tools::root_epsilon<Real>(), Real* error = nullptr, Real* pL1 = nullptr)->decltype(std::declval<F>()(std::declval<Real>()))
{
typedef decltype(f(a)) K;
static_assert(!std::is_integral<K>::value,
"The return type cannot be integral, it must be either a real or complex floating point type.");
static const char* function = "boost::math::quadrature::gauss_kronrod<%1%>::integrate(f, %1%, %1%)";
if (!(boost::math::isnan)(a) && !(boost::math::isnan)(b))
{
@@ -1927,11 +1931,15 @@ public:
if ((boost::math::isfinite)(a) && (boost::math::isfinite)(b))
{
if (b <= a)
if (a==b)
{
return policies::raise_domain_error(function, "Arguments to integrate are in wrong order; integration over [a,b] must have b > a.", a, Policy());
return K(0);
}
recursive_info<F> info = { f, tol };
if (b < a)
{
return -recursive_adaptive_integrate(&info, b, a, max_depth, Real(0), error, pL1);
}
return recursive_adaptive_integrate(&info, a, b, max_depth, Real(0), error, pL1);
}
}

108
linx64/include/boost/math/quadrature/naive_monte_carlo.hpp
View File

@@ -9,7 +9,8 @@
#include <sstream>
#include <algorithm>
#include <vector>
#include <boost/atomic.hpp>
#include <atomic>
#include <memory>
#include <functional>
#include <future>
#include <thread>
@@ -18,7 +19,13 @@
#include <random>
#include <chrono>
#include <map>
#include <type_traits>
#include <boost/math/policies/error_handling.hpp>
#include <boost/math/special_functions/fpclassify.hpp>
#ifdef BOOST_NAIVE_MONTE_CARLO_DEBUG_FAILURES
# include <iostream>
#endif
namespace boost { namespace math { namespace quadrature {
@@ -29,7 +36,8 @@ namespace detail {
DOUBLE_INFINITE};
}
template<class Real, class F, class RandomNumberGenerator = std::mt19937_64, class Policy = boost::math::policies::policy<>>
template<class Real, class F, class RandomNumberGenerator = std::mt19937_64, class Policy = boost::math::policies::policy<>,
typename std::enable_if<std::is_trivially_copyable<Real>::value, bool>::type = true>
class naive_monte_carlo
{
public:
@@ -38,15 +46,17 @@ public:
Real error_goal,
bool singular = true,
uint64_t threads = std::thread::hardware_concurrency(),
uint64_t seed = 0): m_num_threads{threads}, m_seed{seed}
uint64_t seed = 0) noexcept : m_num_threads{threads}, m_seed{seed}, m_volume(1)
{
using std::numeric_limits;
using std::sqrt;
using boost::math::isinf;
uint64_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 (uint64_t i = 0; i < n; ++i)
{
@@ -55,9 +65,9 @@ public:
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 (isinf(bounds[i].first))
{
if (bounds[i].second == numeric_limits<Real>::infinity())
if (isinf(bounds[i].second))
{
m_limit_types[i] = detail::limit_classification::DOUBLE_INFINITE;
}
@@ -69,7 +79,7 @@ public:
m_dxs[i] = numeric_limits<Real>::quiet_NaN();
}
}
else if (bounds[i].second == numeric_limits<Real>::infinity())
else if (isinf(bounds[i].second))
{
m_limit_types[i] = detail::limit_classification::UPPER_BOUND_INFINITE;
if (singular)
@@ -160,10 +170,10 @@ public:
RandomNumberGenerator gen(seed);
Real inv_denom = 1/static_cast<Real>(((gen.max)()-(gen.min)()));
m_num_threads = (std::max)(m_num_threads, (uint64_t) 1);
m_thread_calls.reset(new boost::atomic<uint64_t>[threads]);
m_thread_Ss.reset(new boost::atomic<Real>[threads]);
m_thread_averages.reset(new boost::atomic<Real>[threads]);
m_num_threads = (std::max)(m_num_threads, static_cast<uint64_t>(1));
m_thread_calls.reset(new std::atomic<uint64_t>[threads]);
m_thread_Ss.reset(new std::atomic<Real>[threads]);
m_thread_averages.reset(new std::atomic<Real>[threads]);
Real avg = 0;
for (uint64_t i = 0; i < m_num_threads; ++i)
@@ -282,17 +292,17 @@ private:
m_done = false;
#ifdef BOOST_NAIVE_MONTE_CARLO_DEBUG_FAILURES
std::cout << "Failed to achieve required tolerance first time through..\n";
std::cout << " variance = " << m_variance << std::endl;
std::cout << " average = " << m_avg << std::endl;
std::cout << " total calls = " << m_total_calls << std::endl;
std::cerr << "Failed to achieve required tolerance first time through..\n";
std::cerr << " variance = " << m_variance << std::endl;
std::cerr << " average = " << m_avg << std::endl;
std::cerr << " total calls = " << m_total_calls << std::endl;
for (std::size_t i = 0; i < m_num_threads; ++i)
std::cout << " thread_calls[" << i << "] = " << m_thread_calls[i] << std::endl;
std::cerr << " thread_calls[" << i << "] = " << m_thread_calls[i] << std::endl;
for (std::size_t i = 0; i < m_num_threads; ++i)
std::cout << " thread_averages[" << i << "] = " << m_thread_averages[i] << std::endl;
std::cerr << " thread_averages[" << i << "] = " << m_thread_averages[i] << std::endl;
for (std::size_t i = 0; i < m_num_threads; ++i)
std::cout << " thread_Ss[" << i << "] = " << m_thread_Ss[i] << std::endl;
std::cerr << " thread_Ss[" << i << "] = " << m_thread_Ss[i] << std::endl;
#endif
}
@@ -306,27 +316,27 @@ private:
uint64_t total_calls = 0;
for (uint64_t i = 0; i < m_num_threads; ++i)
{
uint64_t t_calls = m_thread_calls[i].load(boost::memory_order::consume);
uint64_t t_calls = m_thread_calls[i].load(std::memory_order_consume);
total_calls += t_calls;
}
Real variance = 0;
Real avg = 0;
for (uint64_t i = 0; i < m_num_threads; ++i)
{
uint64_t t_calls = m_thread_calls[i].load(boost::memory_order::consume);
uint64_t t_calls = m_thread_calls[i].load(std::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);
avg += m_thread_averages[i].load(std::memory_order_relaxed)*(static_cast<Real>(t_calls) / static_cast<Real>(total_calls));
variance += m_thread_Ss[i].load(std::memory_order_relaxed);
}
m_avg.store(avg, boost::memory_order::release);
m_variance.store(variance / (total_calls - 1), boost::memory_order::release);
m_avg.store(avg, std::memory_order_release);
m_variance.store(variance / (total_calls - 1), std::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 (m_total_calls < 2048 || this->current_error_estimate() > m_error_goal.load(boost::memory_order::consume));
} while (m_total_calls < 2048 || this->current_error_estimate() > m_error_goal.load(std::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(),
@@ -339,7 +349,7 @@ private:
uint64_t total_calls = 0;
for (uint64_t i = 0; i < m_num_threads; ++i)
{
uint64_t t_calls = m_thread_calls[i].load(boost::memory_order::consume);
uint64_t t_calls = m_thread_calls[i].load(std::memory_order_consume);
total_calls += t_calls;
}
Real variance = 0;
@@ -347,13 +357,13 @@ private:
for (uint64_t i = 0; i < m_num_threads; ++i)
{
uint64_t t_calls = m_thread_calls[i].load(boost::memory_order::consume);
uint64_t t_calls = m_thread_calls[i].load(std::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);
avg += m_thread_averages[i].load(std::memory_order_relaxed)*(static_cast<Real>(t_calls) / static_cast<Real>(total_calls));
variance += m_thread_Ss[i].load(std::memory_order_relaxed);
}
m_avg.store(avg, boost::memory_order::release);
m_variance.store(variance / (total_calls - 1), boost::memory_order::release);
m_avg.store(avg, std::memory_order_release);
m_variance.store(variance / (total_calls - 1), std::memory_order_release);
m_total_calls = total_calls; // relaxed store, this is just user feedback
// Sometimes, the master will observe the variance at a very "good" (or bad?) moment,
@@ -362,7 +372,7 @@ private:
}
while ((--max_repeat_tries >= 0) && (this->current_error_estimate() > m_error_goal));
return m_avg.load(boost::memory_order::consume);
return m_avg.load(std::memory_order_consume);
}
void m_thread_monte(uint64_t thread_index, uint64_t seed)
@@ -372,15 +382,15 @@ private:
{
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);
Real inv_denom = static_cast<Real>(1) / static_cast<Real>(( (gen.max)() - (gen.min)() ));
Real M1 = m_thread_averages[thread_index].load(std::memory_order_consume);
Real S = m_thread_Ss[thread_index].load(std::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;
uint64_t k = m_thread_calls[thread_index].load(boost::memory_order::consume);
uint64_t k = m_thread_calls[thread_index].load(std::memory_order_consume);
while (!m_done) // relaxed load
{
int j = 0;
@@ -418,38 +428,40 @@ private:
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);
m_thread_averages[thread_index].store(M1, std::memory_order_release);
m_thread_Ss[thread_index].store(S, std::memory_order_release);
m_thread_calls[thread_index].store(k, std::memory_order_release);
}
}
catch (...)
{
// Signal the other threads that the computation is ruined:
m_done = true; // relaxed store
std::lock_guard<std::mutex> lock(m_exception_mutex); // Scoped lock to prevent race writing to m_exception
m_exception = std::current_exception();
}
}
std::function<Real(std::vector<Real> &)> m_integrand;
uint64_t m_num_threads;
uint64_t m_seed;
boost::atomic<Real> m_error_goal;
boost::atomic<bool> m_done;
std::atomic<uint64_t> m_seed;
std::atomic<Real> m_error_goal;
std::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<uint64_t> m_total_calls;
std::atomic<uint64_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<uint64_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::unique_ptr<std::atomic<uint64_t>[]> m_thread_calls;
std::atomic<Real> m_variance;
std::unique_ptr<std::atomic<Real>[]> m_thread_Ss;
std::atomic<Real> m_avg;
std::unique_ptr<std::atomic<Real>[]> m_thread_averages;
std::chrono::time_point<std::chrono::system_clock> m_start;
std::exception_ptr m_exception;
std::mutex m_exception_mutex;
};
}}}

68
linx64/include/boost/math/quadrature/ooura_fourier_integrals.hpp Normal file
View File

@@ -0,0 +1,68 @@
// Copyright Nick Thompson, 2019
// 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)
/*
* References:
* Ooura, Takuya, and Masatake Mori. "A robust double exponential formula for Fourier-type integrals." Journal of computational and applied mathematics 112.1-2 (1999): 229-241.
* http://www.kurims.kyoto-u.ac.jp/~ooura/intde.html
*/
#ifndef BOOST_MATH_QUADRATURE_OOURA_FOURIER_INTEGRALS_HPP
#define BOOST_MATH_QUADRATURE_OOURA_FOURIER_INTEGRALS_HPP
#include <memory>
#include <boost/math/quadrature/detail/ooura_fourier_integrals_detail.hpp>
namespace boost { namespace math { namespace quadrature {
template<class Real>
class ooura_fourier_sin {
public:
ooura_fourier_sin(const Real relative_error_tolerance = tools::root_epsilon<Real>(), size_t levels = sizeof(Real)) : impl_(std::make_shared<detail::ooura_fourier_sin_detail<Real>>(relative_error_tolerance, levels))
{}
template<class F>
std::pair<Real, Real> integrate(F const & f, Real omega) {
return impl_->integrate(f, omega);
}
// These are just for debugging/unit tests:
std::vector<std::vector<Real>> const & big_nodes() const {
return impl_->big_nodes();
}
std::vector<std::vector<Real>> const & weights_for_big_nodes() const {
return impl_->weights_for_big_nodes();
}
std::vector<std::vector<Real>> const & little_nodes() const {
return impl_->little_nodes();
}
std::vector<std::vector<Real>> const & weights_for_little_nodes() const {
return impl_->weights_for_little_nodes();
}
private:
std::shared_ptr<detail::ooura_fourier_sin_detail<Real>> impl_;
};
template<class Real>
class ooura_fourier_cos {
public:
ooura_fourier_cos(const Real relative_error_tolerance = tools::root_epsilon<Real>(), size_t levels = sizeof(Real)) : impl_(std::make_shared<detail::ooura_fourier_cos_detail<Real>>(relative_error_tolerance, levels))
{}
template<class F>
std::pair<Real, Real> integrate(F const & f, Real omega) {
return impl_->integrate(f, omega);
}
private:
std::shared_ptr<detail::ooura_fourier_cos_detail<Real>> impl_;
};
}}}
#endif

2
linx64/include/boost/math/quadrature/sinh_sinh.hpp
View File

@@ -30,7 +30,7 @@ public:
: m_imp(std::make_shared<detail::sinh_sinh_detail<Real, Policy> >(max_refinements)) {}
template<class F>
auto integrate(const F f, Real tol = boost::math::tools::root_epsilon<Real>(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr)->decltype(std::declval<F>()(std::declval<Real>())) const
auto integrate(const F f, Real tol = boost::math::tools::root_epsilon<Real>(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) const ->decltype(std::declval<F>()(std::declval<Real>()))
{
return m_imp->integrate(f, tol, error, L1, levels);
}

69
linx64/include/boost/math/quadrature/tanh_sinh.hpp
View File

@@ -44,14 +44,14 @@ public:
: m_imp(std::make_shared<detail::tanh_sinh_detail<Real, Policy>>(max_refinements, min_complement)) {}
template<class F>
auto integrate(const F f, Real a, Real b, Real tolerance = tools::root_epsilon<Real>(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) ->decltype(std::declval<F>()(std::declval<Real>())) const;
auto integrate(const F f, Real a, Real b, Real tolerance = tools::root_epsilon<Real>(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) const ->decltype(std::declval<F>()(std::declval<Real>()));
template<class F>
auto integrate(const F f, Real a, Real b, Real tolerance = tools::root_epsilon<Real>(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) ->decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) const;
auto integrate(const F f, Real a, Real b, Real tolerance = tools::root_epsilon<Real>(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) const ->decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>()));
template<class F>
auto integrate(const F f, Real tolerance = tools::root_epsilon<Real>(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) ->decltype(std::declval<F>()(std::declval<Real>())) const;
auto integrate(const F f, Real tolerance = tools::root_epsilon<Real>(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) const ->decltype(std::declval<F>()(std::declval<Real>()));
template<class F>
auto integrate(const F f, Real tolerance = tools::root_epsilon<Real>(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) ->decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) const;
auto integrate(const F f, Real tolerance = tools::root_epsilon<Real>(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) const ->decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>()));
private:
std::shared_ptr<detail::tanh_sinh_detail<Real, Policy>> m_imp;
@@ -59,7 +59,7 @@ private:
template<class Real, class Policy>
template<class F>
auto tanh_sinh<Real, Policy>::integrate(const F f, Real a, Real b, Real tolerance, Real* error, Real* L1, std::size_t* levels) ->decltype(std::declval<F>()(std::declval<Real>())) const
auto tanh_sinh<Real, Policy>::integrate(const F f, Real a, Real b, Real tolerance, Real* error, Real* L1, std::size_t* levels) const ->decltype(std::declval<F>()(std::declval<Real>()))
{
BOOST_MATH_STD_USING
using boost::math::constants::half;
@@ -68,7 +68,8 @@ auto tanh_sinh<Real, Policy>::integrate(const F f, Real a, Real b, Real toleranc
static const char* function = "tanh_sinh<%1%>::integrate";
typedef decltype(std::declval<F>()(std::declval<Real>())) result_type;
static_assert(!std::is_integral<result_type>::value,
"The return type cannot be integral, it must be either a real or complex floating point type.");
if (!(boost::math::isnan)(a) && !(boost::math::isnan)(b))
{
@@ -76,8 +77,8 @@ auto tanh_sinh<Real, Policy>::integrate(const F f, Real a, Real b, Real toleranc
if ((a <= -tools::max_value<Real>()) && (b >= tools::max_value<Real>()))
{
auto u = [&](const Real& t, const Real& tc)->result_type
{
Real t_sq = t*t;
{
Real t_sq = t*t;
Real inv;
if (t > 0.5f)
inv = 1 / ((2 - tc) * tc);
@@ -85,7 +86,7 @@ auto tanh_sinh<Real, Policy>::integrate(const F f, Real a, Real b, Real toleranc
inv = 1 / ((2 + tc) * -tc);
else
inv = 1 / (1 - t_sq);
return f(t*inv)*(1 + t_sq)*inv*inv;
return f(t*inv)*(1 + t_sq)*inv*inv;
};
Real limit = sqrt(tools::min_value<Real>()) * 4;
return m_imp->integrate(u, error, L1, function, limit, limit, tolerance, levels);
@@ -95,7 +96,7 @@ auto tanh_sinh<Real, Policy>::integrate(const F f, Real a, Real b, Real toleranc
if ((boost::math::isfinite)(a) && (b >= tools::max_value<Real>()))
{
auto u = [&](const Real& t, const Real& tc)->result_type
{
{
Real z, arg;
if (t > -0.5f)
z = 1 / (t + 1);
@@ -105,7 +106,7 @@ auto tanh_sinh<Real, Policy>::integrate(const F f, Real a, Real b, Real toleranc
arg = 2 * z + a - 1;
else
arg = a + tc / (2 - tc);
return f(arg)*z*z;
return f(arg)*z*z;
};
Real left_limit = sqrt(tools::min_value<Real>()) * 4;
result_type Q = Real(2) * m_imp->integrate(u, error, L1, function, left_limit, tools::min_value<Real>(), tolerance, levels);
@@ -113,6 +114,10 @@ auto tanh_sinh<Real, Policy>::integrate(const F f, Real a, Real b, Real toleranc
{
*L1 *= 2;
}
if (error)
{
*error *= 2;
}
return Q;
}
@@ -120,7 +125,7 @@ auto tanh_sinh<Real, Policy>::integrate(const F f, Real a, Real b, Real toleranc
if ((boost::math::isfinite)(b) && (a <= -tools::max_value<Real>()))
{
auto v = [&](const Real& t, const Real& tc)->result_type
{
{
Real z;
if (t > -0.5)
z = 1 / (t + 1);
@@ -140,14 +145,22 @@ auto tanh_sinh<Real, Policy>::integrate(const F f, Real a, Real b, Real toleranc
{
*L1 *= 2;
}
if (error)
{
*error *= 2;
}
return Q;
}
if ((boost::math::isfinite)(a) && (boost::math::isfinite)(b))
{
if (b <= a)
if (a == b)
{
return policies::raise_domain_error(function, "Arguments to integrate are in wrong order; integration over [a,b] must have b > a.", a, Policy());
return result_type(0);
}
if (b < a)
{
return -this->integrate(f, b, a, tolerance, error, L1, levels);
}
Real avg = (a + b)*half<Real>();
Real diff = (b - a)*half<Real>();
@@ -156,7 +169,7 @@ auto tanh_sinh<Real, Policy>::integrate(const F f, Real a, Real b, Real toleranc
bool have_small_left = fabs(a) < 0.5f;
bool have_small_right = fabs(b) < 0.5f;
Real left_min_complement = float_next(avg_over_diff_m1) - avg_over_diff_m1;
Real min_complement_limit = (std::max)(tools::min_value<Real>(), Real(tools::min_value<Real>() / diff));
Real min_complement_limit = (std::max)(tools::min_value<Real>(), float_next(Real(tools::min_value<Real>() / diff)));
if (left_min_complement < min_complement_limit)
left_min_complement = min_complement_limit;
Real right_min_complement = avg_over_diff_p1 - float_prior(avg_over_diff_p1);
@@ -170,10 +183,10 @@ auto tanh_sinh<Real, Policy>::integrate(const F f, Real a, Real b, Real toleranc
// factor to move left_min_complement and right_min_complement
// further from the end points of the range.
//
BOOST_ASSERT((left_min_complement * diff + a) > a);
BOOST_ASSERT((b - right_min_complement * diff) < b);
BOOST_MATH_ASSERT((left_min_complement * diff + a) > a);
BOOST_MATH_ASSERT((b - right_min_complement * diff) < b);
auto u = [&](Real z, Real zc)->result_type
{
{
Real position;
if (z < -0.5)
{
@@ -181,7 +194,7 @@ auto tanh_sinh<Real, Policy>::integrate(const F f, Real a, Real b, Real toleranc
return f(diff * (avg_over_diff_m1 - zc));
position = a - diff * zc;
}
if (z > 0.5)
else if (z > 0.5)
{
if(have_small_right)
return f(diff * (avg_over_diff_p1 - zc));
@@ -189,8 +202,8 @@ auto tanh_sinh<Real, Policy>::integrate(const F f, Real a, Real b, Real toleranc
}
else
position = avg + diff*z;
BOOST_ASSERT(position != a);
BOOST_ASSERT(position != b);
BOOST_MATH_ASSERT(position != a);
BOOST_MATH_ASSERT(position != b);
return f(position);
};
result_type Q = diff*m_imp->integrate(u, error, L1, function, left_min_complement, right_min_complement, tolerance, levels);
@@ -199,6 +212,10 @@ auto tanh_sinh<Real, Policy>::integrate(const F f, Real a, Real b, Real toleranc
{
*L1 *= diff;
}
if (error)
{
*error *= diff;
}
return Q;
}
}
@@ -207,7 +224,7 @@ auto tanh_sinh<Real, Policy>::integrate(const F f, Real a, Real b, Real toleranc
template<class Real, class Policy>
template<class F>
auto tanh_sinh<Real, Policy>::integrate(const F f, Real a, Real b, Real tolerance, Real* error, Real* L1, std::size_t* levels) ->decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) const
auto tanh_sinh<Real, Policy>::integrate(const F f, Real a, Real b, Real tolerance, Real* error, Real* L1, std::size_t* levels) const ->decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>()))
{
BOOST_MATH_STD_USING
using boost::math::constants::half;
@@ -237,6 +254,10 @@ auto tanh_sinh<Real, Policy>::integrate(const F f, Real a, Real b, Real toleranc
{
*L1 *= diff;
}
if (error)
{
*error *= diff;
}
return Q;
}
return policies::raise_domain_error(function, "The domain of integration is not sensible; please check the bounds.", a, Policy());
@@ -244,7 +265,7 @@ auto tanh_sinh<Real, Policy>::integrate(const F f, Real a, Real b, Real toleranc
template<class Real, class Policy>
template<class F>
auto tanh_sinh<Real, Policy>::integrate(const F f, Real tolerance, Real* error, Real* L1, std::size_t* levels) ->decltype(std::declval<F>()(std::declval<Real>())) const
auto tanh_sinh<Real, Policy>::integrate(const F f, Real tolerance, Real* error, Real* L1, std::size_t* levels) const ->decltype(std::declval<F>()(std::declval<Real>()))
{
using boost::math::quadrature::detail::tanh_sinh_detail;
static const char* function = "tanh_sinh<%1%>::integrate";
@@ -254,7 +275,7 @@ auto tanh_sinh<Real, Policy>::integrate(const F f, Real tolerance, Real* error,
template<class Real, class Policy>
template<class F>
auto tanh_sinh<Real, Policy>::integrate(const F f, Real tolerance, Real* error, Real* L1, std::size_t* levels) ->decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) const
auto tanh_sinh<Real, Policy>::integrate(const F f, Real tolerance, Real* error, Real* L1, std::size_t* levels) const ->decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>()))
{
using boost::math::quadrature::detail::tanh_sinh_detail;
static const char* function = "tanh_sinh<%1%>::integrate";

42
linx64/include/boost/math/quadrature/trapezoidal.hpp
View File

@@ -17,10 +17,12 @@
#include <cmath>
#include <limits>
#include <utility>
#include <stdexcept>
#include <boost/math/constants/constants.hpp>
#include <boost/math/special_functions/fpclassify.hpp>
#include <boost/math/policies/error_handling.hpp>
#include <boost/math/tools/cxx03_warn.hpp>
namespace boost{ namespace math{ namespace quadrature {
@@ -33,17 +35,24 @@ auto trapezoidal(F f, Real a, Real b, Real tol, std::size_t max_refinements, Rea
// In many math texts, K represents the field of real or complex numbers.
// Too bad we can't put blackboard bold into C++ source!
typedef decltype(f(a)) K;
if(a >= b)
{
return static_cast<K>(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));
}
static_assert(!std::is_integral<K>::value,
"The return type cannot be integral, it must be either a real or complex floating point type.");
if (!(boost::math::isfinite)(a))
{
return static_cast<K>(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 static_cast<K>(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));
return static_cast<K>(boost::math::policies::raise_domain_error(function, "Right endpoint of integration must be finite for adaptive trapezoidal integration but got b = %1%.\n", b, pol));
}
if (a == b)
{
return static_cast<K>(0);
}
if(a > b)
{
return -trapezoidal(f, b, a, tol, max_refinements, error_estimate, L1, pol);
}
@@ -61,11 +70,14 @@ auto trapezoidal(F f, Real a, Real b, Real tol, std::size_t max_refinements, Rea
// 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.
// We want to go through at least 5 levels so we have sampled the function at least 20 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.
// This is of course possible anyway, but 20 samples seems to be a reasonable compromise.
Real error = abs(I0 - I1);
while (k < 4 || (k < max_refinements && error > tol*IL1) )
// I take k < 5, rather than k < 4, or some other smaller minimum number,
// because I hit a truly exceptional bug where the k = 2 and k =3 refinement were bitwise equal,
// but the quadrature had not yet converged.
while (k < 5 || (k < max_refinements && error > tol*IL1) )
{
I0 = I1;
IL0 = IL1;
@@ -102,22 +114,10 @@ auto trapezoidal(F f, Real a, Real b, Real tol, std::size_t max_refinements, Rea
return static_cast<K>(I1);
}
#if BOOST_WORKAROUND(BOOST_MSVC, < 1800)
// Template argument dedcution failure otherwise:
template<class F, class Real>
auto trapezoidal(F f, Real a, Real b, Real tol = 0, std::size_t max_refinements = 12, Real* error_estimate = 0, Real* L1 = 0)->decltype(std::declval<F>()(std::declval<Real>()))
#elif !defined(BOOST_NO_CXX11_NULLPTR)
template<class F, class Real>
auto trapezoidal(F f, Real a, Real b, Real tol = boost::math::tools::root_epsilon<Real>(), std::size_t max_refinements = 12, Real* error_estimate = nullptr, Real* L1 = nullptr)->decltype(std::declval<F>()(std::declval<Real>()))
#else
template<class F, class Real>
auto trapezoidal(F f, Real a, Real b, Real tol = boost::math::tools::root_epsilon<Real>(), std::size_t max_refinements = 12, Real* error_estimate = 0, Real* L1 = 0)->decltype(std::declval<F>()(std::declval<Real>()))
#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<>());
}

47
linx64/include/boost/math/quadrature/wavelet_transforms.hpp Normal file
View File

@@ -0,0 +1,47 @@
/*
* Copyright Nick Thompson, 2020
* 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_WAVELET_TRANSFORMS_HPP
#define BOOST_MATH_QUADRATURE_WAVELET_TRANSFORMS_HPP
#include <boost/math/special_functions/daubechies_wavelet.hpp>
#include <boost/math/quadrature/trapezoidal.hpp>
namespace boost::math::quadrature {
template<class F, typename Real, int p>
class daubechies_wavelet_transform
{
public:
daubechies_wavelet_transform(F f, int grid_refinements = -1, Real tol = 100*std::numeric_limits<Real>::epsilon(),
int max_refinements = 12) : f_{f}, psi_(grid_refinements), tol_{tol}, max_refinements_{max_refinements}
{}
daubechies_wavelet_transform(F f, boost::math::daubechies_wavelet<Real, p> wavelet, Real tol = 100*std::numeric_limits<Real>::epsilon(),
int max_refinements = 12) : f_{f}, psi_{wavelet}, tol_{tol}, max_refinements_{max_refinements}
{}
auto operator()(Real s, Real t) const ->decltype(std::declval<F>()(std::declval<Real>()))
{
using std::sqrt;
using std::abs;
using boost::math::quadrature::trapezoidal;
auto g = [&] (Real u) {
return f_(s*u+t)*psi_(u);
};
auto [a,b] = psi_.support();
return sqrt(abs(s))*trapezoidal(g, a, b, tol_, max_refinements_);
}
private:
F f_;
boost::math::daubechies_wavelet<Real, p> psi_;
Real tol_;
int max_refinements_;
};
}
#endif