Added thirdparty: boost library

thirdparty/include/boost/math/quadrature/detail/ooura_fourier_integrals_detail.hpp

@@ -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

thirdparty/include/boost/math/quadrature/exp_sinh.hpp

@@ -0,0 +1,102 @@
+// Copyright Nick Thompson, 2017
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+/*
+ * This class performs exp-sinh quadrature on half infinite intervals.
+ *
+ * References:
+ *
+ * 1) Tanaka, Ken'ichiro, et al. "Function classes for double exponential integration formulas." Numerische Mathematik 111.4 (2009): 631-655.
+ */
+
+#ifndef BOOST_MATH_QUADRATURE_EXP_SINH_HPP
+#define BOOST_MATH_QUADRATURE_EXP_SINH_HPP
+
+#include <cmath>
+#include <limits>
+#include <memory>
+#include <string>
+#include <boost/math/quadrature/detail/exp_sinh_detail.hpp>
+
+namespace boost{ namespace math{ namespace quadrature {
+
+template<class Real, class Policy = policies::policy<> >
+class exp_sinh
+{
+public:
+   exp_sinh(size_t max_refinements = 9)
+      : 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) 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) const ->decltype(std::declval<F>()(std::declval<Real>()));
+
+private:
+    std::shared_ptr<detail::exp_sinh_detail<Real, Policy>> m_imp;
+};
+
+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) 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;
+
+    static const char* function = "boost::math::quadrature::exp_sinh<%1%>::integrate";
+
+    // Neither limit may be a NaN:
+    if((boost::math::isnan)(a) || (boost::math::isnan)(b))
+    {
+       return static_cast<K>(policies::raise_domain_error(function, "NaN supplied as one limit of integration - sorry I don't know what to do", a, Policy()));
+     }
+    // Right limit is infinite:
+    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 == static_cast<Real>(0))
+        {
+            return m_imp->integrate(f, error, L1, function, tolerance, levels);
+        }
+        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)->K { return f(b-t);};
+        return m_imp->integrate(u, error, L1, function, tolerance, levels);
+    }
+
+    // Infinite limits:
+    if ((a <= -boost::math::tools::max_value<Real>()) && (b >= boost::math::tools::max_value<Real>()))
+    {
+        return static_cast<K>(policies::raise_domain_error(function, "Use sinh_sinh quadrature for integration over the whole real line; exp_sinh is for half infinite integrals.", a, Policy()));
+    }
+    // If we get to here then both ends must necessarily be finite:
+    return static_cast<K>(policies::raise_domain_error(function, "Use tanh_sinh quadrature for integration over finite domains; exp_sinh is for half infinite integrals.", a, Policy()));
+}
+
+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) 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);
+}
+
+
+}}}
+#endif

thirdparty/include/boost/math/quadrature/naive_monte_carlo.hpp

@@ -0,0 +1,468 @@
+/*
+ * Copyright Nick Thompson, 2018
+ * Use, modification and distribution are subject to the
+ * Boost Software License, Version 1.0. (See accompanying file
+ * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+#ifndef BOOST_MATH_QUADRATURE_NAIVE_MONTE_CARLO_HPP
+#define BOOST_MATH_QUADRATURE_NAIVE_MONTE_CARLO_HPP
+#include <sstream>
+#include <algorithm>
+#include <vector>
+#include <atomic>
+#include <memory>
+#include <functional>
+#include <future>
+#include <thread>
+#include <initializer_list>
+#include <utility>
+#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 {
+
+namespace detail {
+  enum class limit_classification {FINITE,
+                                   LOWER_BOUND_INFINITE,
+                                   UPPER_BOUND_INFINITE,
+                                   DOUBLE_INFINITE};
+}
+
+template<class Real, class F, class RandomNumberGenerator = std::mt19937_64, class Policy = boost::math::policies::policy<>,
+         typename std::enable_if<std::is_trivially_copyable<Real>::value, bool>::type = true>
+class naive_monte_carlo
+{
+public:
+    naive_monte_carlo(const F& integrand,
+                      std::vector<std::pair<Real, Real>> const & bounds,
+                      Real error_goal,
+                      bool singular = true,
+                      uint64_t threads = std::thread::hardware_concurrency(),
+                      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);
+
+        static const char* function = "boost::math::quadrature::naive_monte_carlo<%1%>";
+        for (uint64_t i = 0; i < n; ++i)
+        {
+            if (bounds[i].second <= bounds[i].first)
+            {
+                boost::math::policies::raise_domain_error(function, "The upper bound is <= the lower bound.\n", bounds[i].second, Policy());
+                return;
+            }
+            if (isinf(bounds[i].first))
+            {
+                if (isinf(bounds[i].second))
+                {
+                    m_limit_types[i] = detail::limit_classification::DOUBLE_INFINITE;
+                }
+                else
+                {
+                    m_limit_types[i] = detail::limit_classification::LOWER_BOUND_INFINITE;
+                    // Ok ok this is bad to use the second bound as the lower limit and then reflect.
+                    m_lbs[i] = bounds[i].second;
+                    m_dxs[i] = numeric_limits<Real>::quiet_NaN();
+                }
+            }
+            else if (isinf(bounds[i].second))
+            {
+                m_limit_types[i] = detail::limit_classification::UPPER_BOUND_INFINITE;
+                if (singular)
+                {
+                    // I've found that it's easier to sample on a closed set and perturb the boundary
+                    // than to try to sample very close to the boundary.
+                    m_lbs[i] = std::nextafter(bounds[i].first, (std::numeric_limits<Real>::max)());
+                }
+                else
+                {
+                    m_lbs[i] = bounds[i].first;
+                }
+                m_dxs[i] = numeric_limits<Real>::quiet_NaN();
+            }
+            else
+            {
+                m_limit_types[i] = detail::limit_classification::FINITE;
+                if (singular)
+                {
+                    if (bounds[i].first == 0)
+                    {
+                        m_lbs[i] = std::numeric_limits<Real>::epsilon();
+                    }
+                    else
+                    {
+                        m_lbs[i] = std::nextafter(bounds[i].first, (std::numeric_limits<Real>::max)());
+                    }
+
+                    m_dxs[i] = std::nextafter(bounds[i].second, std::numeric_limits<Real>::lowest()) - m_lbs[i];
+                }
+                else
+                {
+                    m_lbs[i] = bounds[i].first;
+                    m_dxs[i] = bounds[i].second - bounds[i].first;
+                }
+                m_volume *= m_dxs[i];
+            }
+        }
+
+        m_integrand = [this, &integrand](std::vector<Real> & x)->Real
+        {
+            Real coeff = m_volume;
+            for (uint64_t i = 0; i < x.size(); ++i)
+            {
+                // Variable transformation are listed at:
+                // https://en.wikipedia.org/wiki/Numerical_integration
+                // However, we've made some changes to these so that we can evaluate on a compact domain.
+                if (m_limit_types[i] == detail::limit_classification::FINITE)
+                {
+                    x[i] = m_lbs[i] + x[i]*m_dxs[i];
+                }
+                else if (m_limit_types[i] == detail::limit_classification::UPPER_BOUND_INFINITE)
+                {
+                    Real t = x[i];
+                    Real z = 1/(1 + numeric_limits<Real>::epsilon() - t);
+                    coeff *= (z*z)*(1 + numeric_limits<Real>::epsilon());
+                    x[i] = m_lbs[i] + t*z;
+                }
+                else if (m_limit_types[i] == detail::limit_classification::LOWER_BOUND_INFINITE)
+                {
+                    Real t = x[i];
+                    Real z = 1/(t+sqrt((numeric_limits<Real>::min)()));
+                    coeff *= (z*z);
+                    x[i] = m_lbs[i] + (t-1)*z;
+                }
+                else
+                {
+                    Real t1 = 1/(1+numeric_limits<Real>::epsilon() - x[i]);
+                    Real t2 = 1/(x[i]+numeric_limits<Real>::epsilon());
+                    x[i] = (2*x[i]-1)*t1*t2/4;
+                    coeff *= (t1*t1+t2*t2)/4;
+                }
+            }
+            return coeff*integrand(x);
+        };
+
+        // If we don't do a single function call in the constructor,
+        // we can't do a restart.
+        std::vector<Real> x(m_lbs.size());
+
+        // If the seed is zero, that tells us to choose a random seed for the user:
+        if (seed == 0)
+        {
+            std::random_device rd;
+            seed = rd();
+        }
+
+        RandomNumberGenerator gen(seed);
+        Real inv_denom = 1/static_cast<Real>(((gen.max)()-(gen.min)()));
+
+        m_num_threads = (std::max)(m_num_threads, 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)
+        {
+            for (uint64_t j = 0; j < m_lbs.size(); ++j)
+            {
+                x[j] = (gen()-(gen.min)())*inv_denom;
+            }
+            Real y = m_integrand(x);
+            m_thread_averages[i] = y; // relaxed store
+            m_thread_calls[i] = 1;
+            m_thread_Ss[i] = 0;
+            avg += y;
+        }
+        avg /= m_num_threads;
+        m_avg = avg; // relaxed store
+
+        m_error_goal = error_goal; // relaxed store
+        m_start = std::chrono::system_clock::now();
+        m_done = false; // relaxed store
+        m_total_calls = m_num_threads;  // relaxed store
+        m_variance = (numeric_limits<Real>::max)();
+    }
+
+    std::future<Real> integrate()
+    {
+        // Set done to false in case we wish to restart:
+        m_done.store(false); // relaxed store, no worker threads yet
+        m_start = std::chrono::system_clock::now();
+        return std::async(std::launch::async,
+                          &naive_monte_carlo::m_integrate, this);
+    }
+
+    void cancel()
+    {
+        // If seed = 0 (meaning have the routine pick the seed), this leaves the seed the same.
+        // If seed != 0, then the seed is changed, so a restart doesn't do the exact same thing.
+        m_seed = m_seed*m_seed;
+        m_done = true; // relaxed store, worker threads will get the message eventually
+        // Make sure the error goal is infinite, because otherwise we'll loop when we do the final error goal check:
+        m_error_goal = (std::numeric_limits<Real>::max)();
+    }
+
+    Real variance() const
+    {
+        return m_variance.load();
+    }
+
+    Real current_error_estimate() const
+    {
+        using std::sqrt;
+        //
+        // There is a bug here: m_variance and m_total_calls get updated asynchronously
+        // and may be out of synch when we compute the error estimate, not sure if it matters though...
+        //
+        return sqrt(m_variance.load()/m_total_calls.load());
+    }
+
+    std::chrono::duration<Real> estimated_time_to_completion() const
+    {
+        auto now = std::chrono::system_clock::now();
+        std::chrono::duration<Real> elapsed_seconds = now - m_start;
+        Real r = this->current_error_estimate()/m_error_goal.load(); // relaxed load
+        if (r*r <= 1) {
+            return 0*elapsed_seconds;
+        }
+        return (r*r - 1)*elapsed_seconds;
+    }
+
+    void update_target_error(Real new_target_error)
+    {
+        m_error_goal = new_target_error;  // relaxed store
+    }
+
+    Real progress() const
+    {
+        Real r = m_error_goal.load()/this->current_error_estimate();  // relaxed load
+        if (r*r >= 1)
+        {
+            return 1;
+        }
+        return r*r;
+    }
+
+    Real current_estimate() const
+    {
+        return m_avg.load();
+    }
+
+    uint64_t calls() const
+    {
+        return m_total_calls.load();  // relaxed load
+    }
+
+private:
+
+   Real m_integrate()
+   {
+      uint64_t seed;
+      // If the user tells us to pick a seed, pick a seed:
+      if (m_seed == 0)
+      {
+         std::random_device rd;
+         seed = rd();
+      }
+      else // use the seed we are given:
+      {
+         seed = m_seed;
+      }
+      RandomNumberGenerator gen(seed);
+      int max_repeat_tries = 5;
+      do{
+
+         if (max_repeat_tries < 5)
+         {
+            m_done = false;
+
+#ifdef BOOST_NAIVE_MONTE_CARLO_DEBUG_FAILURES
+            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::cerr << "  thread_calls[" << i << "] = " << m_thread_calls[i] << std::endl;
+            for (std::size_t i = 0; i < m_num_threads; ++i)
+               std::cerr << "  thread_averages[" << i << "] = " << m_thread_averages[i] << std::endl;
+            for (std::size_t i = 0; i < m_num_threads; ++i)
+               std::cerr << "  thread_Ss[" << i << "] = " << m_thread_Ss[i] << std::endl;
+#endif
+         }
+
+         std::vector<std::thread> threads(m_num_threads);
+         for (uint64_t i = 0; i < threads.size(); ++i)
+         {
+            threads[i] = std::thread(&naive_monte_carlo::m_thread_monte, this, i, gen());
+         }
+         do {
+            std::this_thread::sleep_for(std::chrono::milliseconds(100));
+            uint64_t total_calls = 0;
+            for (uint64_t i = 0; i < m_num_threads; ++i)
+            {
+               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(std::memory_order_consume);
+               // Will this overflow? Not hard to remove . . .
+               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, 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(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(),
+            std::mem_fn(&std::thread::join));
+         if (m_exception)
+         {
+            std::rethrow_exception(m_exception);
+         }
+         // Incorporate their work into the final estimate:
+         uint64_t total_calls = 0;
+         for (uint64_t i = 0; i < m_num_threads; ++i)
+         {
+            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(std::memory_order_consume);
+            // Averages weighted by the number of calls the thread made:
+            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, 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,
+         // Then the threads proceed to find the variance is much greater by the time they hear the message to stop.
+         // This *WOULD* make sure that the final error estimate is within the error bounds.
+      }
+      while ((--max_repeat_tries >= 0) && (this->current_error_estimate() > m_error_goal));
+
+      return m_avg.load(std::memory_order_consume);
+    }
+
+    void m_thread_monte(uint64_t thread_index, uint64_t seed)
+    {
+        using std::numeric_limits;
+        try
+        {
+            std::vector<Real> x(m_lbs.size());
+            RandomNumberGenerator gen(seed);
+            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(std::memory_order_consume);
+            while (!m_done) // relaxed load
+            {
+                int j = 0;
+                // If we don't have a certain number of calls before an update, we can easily terminate prematurely
+                // because the variance estimate is way too low. This magic number is a reasonable compromise, as 1/sqrt(2048) = 0.02,
+                // so it should recover 2 digits if the integrand isn't poorly behaved, and if it is, it should discover that before premature termination.
+                // Of course if the user has 64 threads, then this number is probably excessive.
+                int magic_calls_before_update = 2048;
+                while (j++ < magic_calls_before_update)
+                {
+                    for (uint64_t i = 0; i < m_lbs.size(); ++i)
+                    {
+                        x[i] = (gen() - (gen.min)())*inv_denom;
+                    }
+                    Real f = m_integrand(x);
+                    using std::isfinite;
+                    if (!isfinite(f))
+                    {
+                        // The call to m_integrand transform x, so this error message states the correct node.
+                        std::stringstream os;
+                        os << "Your integrand was evaluated at {";
+                        for (uint64_t i = 0; i < x.size() -1; ++i)
+                        {
+                             os << x[i] << ", ";
+                        }
+                        os << x[x.size() -1] << "}, and returned " << f << std::endl;
+                        static const char* function = "boost::math::quadrature::naive_monte_carlo<%1%>";
+                        boost::math::policies::raise_domain_error(function, os.str().c_str(), /*this is a dummy arg to make it compile*/ 7.2, Policy());
+                    }
+                    ++k;
+                    Real term = (f - M1)/k;
+                    Real y1 = term - compensator;
+                    Real M2 = M1 + y1;
+                    compensator = (M2 - M1) - y1;
+                    S += (f - M1)*(f - M2);
+                    M1 = M2;
+                }
+                m_thread_averages[thread_index].store(M1, 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;
+    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;
+    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<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;
+};
+
+}}}
+#endif

thirdparty/include/boost/math/quadrature/ooura_fourier_integrals.hpp

@@ -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

thirdparty/include/boost/math/quadrature/sinh_sinh.hpp

@@ -0,0 +1,43 @@
+// Copyright Nick Thompson, 2017
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+/*
+ * This class performs sinh-sinh quadrature over the entire real line.
+ *
+ * References:
+ *
+ * 1) Tanaka, Ken'ichiro, et al. "Function classes for double exponential integration formulas." Numerische Mathematik 111.4 (2009): 631-655.
+ */
+
+#ifndef BOOST_MATH_QUADRATURE_SINH_SINH_HPP
+#define BOOST_MATH_QUADRATURE_SINH_SINH_HPP
+
+#include <cmath>
+#include <limits>
+#include <memory>
+#include <boost/math/quadrature/detail/sinh_sinh_detail.hpp>
+
+namespace boost{ namespace math{ namespace quadrature {
+
+template<class Real, class Policy = boost::math::policies::policy<> >
+class sinh_sinh
+{
+public:
+    sinh_sinh(size_t max_refinements = 9)
+        : 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) const ->decltype(std::declval<F>()(std::declval<Real>()))
+    {
+        return m_imp->integrate(f, tol, error, L1, levels);
+    }
+
+private:
+    std::shared_ptr<detail::sinh_sinh_detail<Real, Policy>> m_imp;
+};
+
+}}}
+#endif

thirdparty/include/boost/math/quadrature/tanh_sinh.hpp

@@ -0,0 +1,289 @@
+// Copyright Nick Thompson, 2017
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+/*
+ * This class performs tanh-sinh quadrature on the real line.
+ * Tanh-sinh quadrature is exponentially convergent for integrands in Hardy spaces,
+ * (see https://en.wikipedia.org/wiki/Hardy_space for a formal definition), and is optimal for a random function from that class.
+ *
+ * The tanh-sinh quadrature is one of a class of so called "double exponential quadratures"-there is a large family of them,
+ * but this one seems to be the most commonly used.
+ *
+ * As always, there are caveats: For instance, if the function you want to integrate is not holomorphic on the unit disk,
+ * then the rapid convergence will be spoiled. In this case, a more appropriate quadrature is (say) Romberg, which does not
+ * require the function to be holomorphic, only differentiable up to some order.
+ *
+ * In addition, if you are integrating a periodic function over a period, the trapezoidal rule is better.
+ *
+ * References:
+ *
+ * 1) Mori, Masatake. "Quadrature formulas obtained by variable transformation and the DE-rule." Journal of Computational and Applied Mathematics 12 (1985): 119-130.
+ * 2) Bailey, David H., Karthik Jeyabalan, and Xiaoye S. Li. "A comparison of three high-precision quadrature schemes." Experimental Mathematics 14.3 (2005): 317-329.
+ * 3) Press, William H., et al. "Numerical recipes third edition: the art of scientific computing." Cambridge University Press 32 (2007): 10013-2473.
+ *
+ */
+
+#ifndef BOOST_MATH_QUADRATURE_TANH_SINH_HPP
+#define BOOST_MATH_QUADRATURE_TANH_SINH_HPP
+
+#include <cmath>
+#include <limits>
+#include <memory>
+#include <boost/math/quadrature/detail/tanh_sinh_detail.hpp>
+
+namespace boost{ namespace math{ namespace quadrature {
+
+template<class Real, class Policy = policies::policy<> >
+class tanh_sinh
+{
+public:
+    tanh_sinh(size_t max_refinements = 15, const Real& min_complement = tools::min_value<Real>() * 4)
+    : 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) 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) 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) 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) const ->decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>()));
+
+private:
+    std::shared_ptr<detail::tanh_sinh_detail<Real, Policy>> m_imp;
+};
+
+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) const ->decltype(std::declval<F>()(std::declval<Real>()))
+{
+    BOOST_MATH_STD_USING
+    using boost::math::constants::half;
+    using boost::math::quadrature::detail::tanh_sinh_detail;
+
+    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))
+    {
+
+       // Infinite limits:
+       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 inv;
+             if (t > 0.5f)
+                inv = 1 / ((2 - tc) * tc);
+             else if(t < -0.5)
+                inv = 1 / ((2 + tc) * -tc);
+             else
+                inv = 1 / (1 - t_sq);
+             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);
+       }
+
+       // Right limit is infinite:
+       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);
+             else
+                z = -1 / tc;
+             if (t < 0.5)
+                arg = 2 * z + a - 1;
+             else
+                arg = a + tc / (2 - tc);
+             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);
+          if (L1)
+          {
+             *L1 *= 2;
+          }
+          if (error)
+          {
+             *error *= 2;
+          }
+
+          return Q;
+       }
+
+       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);
+             else
+                z = -1 / tc;
+             Real arg;
+             if (t < 0.5)
+                arg = 2 * z - 1;
+             else
+                arg = tc / (2 - tc);
+             return f(b - arg) * z * z;
+          };
+
+          Real left_limit = sqrt(tools::min_value<Real>()) * 4;
+          result_type Q = Real(2) * m_imp->integrate(v, error, L1, function, left_limit, tools::min_value<Real>(), tolerance, levels);
+          if (L1)
+          {
+             *L1 *= 2;
+          }
+          if (error)
+          {
+             *error *= 2;
+          }
+          return Q;
+       }
+
+       if ((boost::math::isfinite)(a) && (boost::math::isfinite)(b))
+       {
+          if (a == b)
+          {
+             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>();
+          Real avg_over_diff_m1 = a / diff;
+          Real avg_over_diff_p1 = b / diff;
+          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>(), 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);
+          if (right_min_complement < min_complement_limit)
+             right_min_complement = min_complement_limit;
+          //
+          // These asserts will fail only if rounding errors on
+          // type Real have accumulated so much error that it's
+          // broken our internal logic.  Should that prove to be
+          // a persistent issue, we might need to add a bit of fudge
+          // factor to move left_min_complement and right_min_complement
+          // further from the end points of the range.
+          //
+          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)
+             {
+                if(have_small_left)
+                  return f(diff * (avg_over_diff_m1 - zc));
+                position = a - diff * zc;
+             }
+             else if (z > 0.5)
+             {
+                if(have_small_right)
+                  return f(diff * (avg_over_diff_p1 - zc));
+                position = b - diff * zc;
+             }
+             else
+                position = avg + diff*z;
+             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);
+
+          if (L1)
+          {
+             *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());
+}
+
+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) const ->decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>()))
+{
+   BOOST_MATH_STD_USING
+      using boost::math::constants::half;
+   using boost::math::quadrature::detail::tanh_sinh_detail;
+
+   static const char* function = "tanh_sinh<%1%>::integrate";
+
+   if ((boost::math::isfinite)(a) && (boost::math::isfinite)(b))
+   {
+      if (b <= a)
+      {
+         return policies::raise_domain_error(function, "Arguments to integrate are in wrong order; integration over [a,b] must have b > a.", a, Policy());
+      }
+      auto u = [&](Real z, Real zc)->Real
+      {
+         if (z < 0)
+            return f((a - b) * zc / 2 + a, (b - a) * zc / 2);
+         else
+            return f((a - b) * zc / 2 + b, (b - a) * zc / 2);
+      };
+      Real diff = (b - a)*half<Real>();
+      Real left_min_complement = tools::min_value<Real>() * 4;
+      Real right_min_complement = tools::min_value<Real>() * 4;
+      Real Q = diff*m_imp->integrate(u, error, L1, function, left_min_complement, right_min_complement, tolerance, levels);
+
+      if (L1)
+      {
+         *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());
+}
+
+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) const ->decltype(std::declval<F>()(std::declval<Real>()))
+{
+   using boost::math::quadrature::detail::tanh_sinh_detail;
+   static const char* function = "tanh_sinh<%1%>::integrate";
+   Real min_complement = tools::epsilon<Real>();
+   return m_imp->integrate([&](const Real& arg, const Real&) { return f(arg); }, error, L1, function, min_complement, min_complement, tolerance, levels);
+}
+
+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) 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";
+   Real min_complement = tools::min_value<Real>() * 4;
+   return m_imp->integrate(f, error, L1, function, min_complement, min_complement, tolerance, levels);
+}
+
+}
+}
+}
+#endif

thirdparty/include/boost/math/quadrature/trapezoidal.hpp

@@ -0,0 +1,125 @@
+/*
+ * Copyright Nick Thompson, 2017
+ * Use, modification and distribution are subject to the
+ * Boost Software License, Version 1.0. (See accompanying file
+ * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ *
+ * Use the adaptive trapezoidal rule to estimate the integral of periodic functions over a period,
+ * or to integrate a function whose derivative vanishes at the endpoints.
+ *
+ * If your function does not satisfy these conditions, and instead is simply continuous and bounded
+ * over the whole interval, then this routine will still converge, albeit slowly. However, there
+ * are much more efficient methods in this case, including Romberg, Simpson, and double exponential quadrature.
+ */
+
+#ifndef BOOST_MATH_QUADRATURE_TRAPEZOIDAL_HPP
+#define BOOST_MATH_QUADRATURE_TRAPEZOIDAL_HPP
+
+#include <cmath>
+#include <limits>
+#include <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 {
+
+template<class F, class Real, class Policy>
+auto trapezoidal(F f, Real a, Real b, Real tol, std::size_t max_refinements, Real* error_estimate, Real* L1, const Policy& pol)->decltype(std::declval<F>()(std::declval<Real>()))
+{
+    static const char* function = "boost::math::quadrature::trapezoidal<%1%>(F, %1%, %1%, %1%)";
+    using std::abs;
+    using boost::math::constants::half;
+    // 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;
+    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 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);
+    }
+
+
+    K ya = f(a);
+    K yb = f(b);
+    Real h = (b - a)*half<Real>();
+    K I0 = (ya + yb)*h;
+    Real IL0 = (abs(ya) + abs(yb))*h;
+
+    K yh = f(a + h);
+    K I1;
+    I1 = I0*half<Real>() + yh*h;
+    Real IL1 = IL0*half<Real>() + abs(yh)*h;
+
+    // The recursion is:
+    // I_k = 1/2 I_{k-1} + 1/2^k \sum_{j=1; j odd, j < 2^k} f(a + j(b-a)/2^k)
+    std::size_t k = 2;
+    // We want to go through at least 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 20 samples seems to be a reasonable compromise.
+    Real error = abs(I0 - I1);
+    // 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;
+
+        I1 = I0*half<Real>();
+        IL1 = IL0*half<Real>();
+        std::size_t p = static_cast<std::size_t>(1u) << k;
+        h *= half<Real>();
+        K sum = 0;
+        Real absum = 0;
+
+        for(std::size_t j = 1; j < p; j += 2)
+        {
+            K y = f(a + j*h);
+            sum += y;
+            absum += abs(y);
+        }
+
+        I1 += sum*h;
+        IL1 += absum*h;
+        ++k;
+        error = abs(I0 - I1);
+    }
+
+    if (error_estimate)
+    {
+        *error_estimate = error;
+    }
+
+    if (L1)
+    {
+        *L1 = IL1;
+    }
+
+    return static_cast<K>(I1);
+}
+
+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>()))
+{
+   return trapezoidal(f, a, b, tol, max_refinements, error_estimate, L1, boost::math::policies::policy<>());
+}
+
+}}}
+#endif

thirdparty/include/boost/math/quadrature/wavelet_transforms.hpp

@@ -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