Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:39:53

0001 // Copyright Nick Thompson, 2019
0002 // Use, modification and distribution are subject to the
0003 // Boost Software License, Version 1.0.
0004 // (See accompanying file LICENSE_1_0.txt
0005 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0006 #ifndef BOOST_MATH_QUADRATURE_DETAIL_OOURA_FOURIER_INTEGRALS_DETAIL_HPP
0007 #define BOOST_MATH_QUADRATURE_DETAIL_OOURA_FOURIER_INTEGRALS_DETAIL_HPP
0008 #include <utility> // for std::pair.
0009 #include <vector>
0010 #include <iostream>
0011 #include <boost/math/special_functions/expm1.hpp>
0012 #include <boost/math/special_functions/sin_pi.hpp>
0013 #include <boost/math/special_functions/cos_pi.hpp>
0014 #include <boost/math/constants/constants.hpp>
0015 #include <boost/math/tools/config.hpp>
0016 
0017 #ifdef BOOST_HAS_THREADS
0018 #include <mutex>
0019 #include <atomic>
0020 #endif
0021 
0022 namespace boost { namespace math { namespace quadrature { namespace detail {
0023 
0024 // Ooura and Mori, A robust double exponential formula for Fourier-type integrals,
0025 // eta is the argument to the exponential in equation 3.3:
0026 template<class Real>
0027 std::pair<Real, Real> ooura_eta(Real x, Real alpha) {
0028     using std::expm1;
0029     using std::exp;
0030     using std::abs;
0031     Real expx = exp(x);
0032     Real eta_prime = 2 + alpha/expx + expx/4;
0033     Real eta;
0034     // This is the fast branch:
0035     if (abs(x) > 0.125) {
0036         eta = 2*x - alpha*(1/expx - 1) + (expx - 1)/4;
0037     }
0038     else {// this is the slow branch using expm1 for small x:
0039         eta = 2*x - alpha*expm1(-x) + expm1(x)/4;
0040     }
0041     return {eta, eta_prime};
0042 }
0043 
0044 // Ooura and Mori, A robust double exponential formula for Fourier-type integrals,
0045 // equation 3.6:
0046 template<class Real>
0047 Real calculate_ooura_alpha(Real h)
0048 {
0049     using boost::math::constants::pi;
0050     using std::log1p;
0051     using std::sqrt;
0052     Real x = sqrt(16 + 4*log1p(pi<Real>()/h)/h);
0053     return 1/x;
0054 }
0055 
0056 template<class Real>
0057 std::pair<Real, Real> ooura_sin_node_and_weight(long n, Real h, Real alpha)
0058 {
0059     using std::expm1;
0060     using std::exp;
0061     using std::abs;
0062     using boost::math::constants::pi;
0063     using std::isnan;
0064 
0065     if (n == 0) {
0066         // Equation 44 of https://arxiv.org/pdf/0911.4796.pdf
0067         // Fourier Transform of the Stretched Exponential Function: Analytic Error Bounds,
0068         // Double Exponential Transform, and Open-Source Implementation,
0069         // Joachim Wuttke, 
0070         // The C library libkww provides functions to compute the Kohlrausch-Williams-Watts function, 
0071         // the Laplace-Fourier transform of the stretched (or compressed) exponential function exp(-t^beta)
0072         // for exponent beta between 0.1 and 1.9 with sixteen decimal digits accuracy.
0073 
0074         Real eta_prime_0 = Real(2) + alpha + Real(1)/Real(4);
0075         Real node = pi<Real>()/(eta_prime_0*h);
0076         Real weight = pi<Real>()*boost::math::sin_pi(1/(eta_prime_0*h));
0077         Real eta_dbl_prime = -alpha + Real(1)/Real(4);
0078         Real phi_prime_0 = (1 - eta_dbl_prime/(eta_prime_0*eta_prime_0))/2;
0079         weight *= phi_prime_0;
0080         return {node, weight};
0081     }
0082     Real x = n*h;
0083     auto p = ooura_eta(x, alpha);
0084     auto eta = p.first;
0085     auto eta_prime = p.second;
0086 
0087     Real expm1_meta = expm1(-eta);
0088     Real exp_meta = exp(-eta);
0089     Real node = -n*pi<Real>()/expm1_meta;
0090 
0091 
0092     // I have verified that this is not a significant source of inaccuracy in the weight computation:
0093     Real phi_prime = -(expm1_meta + x*exp_meta*eta_prime)/(expm1_meta*expm1_meta);
0094 
0095     // The main source of inaccuracy is in computation of sin_pi.
0096     // But I've agonized over this, and I think it's as good as it can get:
0097     Real s = pi<Real>();
0098     Real arg;
0099     if(eta > 1) {
0100         arg = n/( 1/exp_meta - 1 );
0101         s *= boost::math::sin_pi(arg);
0102         if (n&1) {
0103             s *= -1;
0104         }
0105     }
0106     else if (eta < -1) {
0107         arg = n/(1-exp_meta);
0108         s *= boost::math::sin_pi(arg);
0109     }
0110     else {
0111         arg = -n*exp_meta/expm1_meta;
0112         s *= boost::math::sin_pi(arg);
0113         if (n&1) {
0114             s *= -1;
0115         }
0116     }
0117 
0118     Real weight = s*phi_prime;
0119     return {node, weight};
0120 }
0121 
0122 #ifdef BOOST_MATH_INSTRUMENT_OOURA
0123 template<class Real>
0124 void print_ooura_estimate(size_t i, Real I0, Real I1, Real omega) {
0125     using std::abs;
0126     std::cout << std::defaultfloat
0127               << std::setprecision(std::numeric_limits<Real>::digits10)
0128               << std::fixed;
0129     std::cout << "h = " << Real(1)/Real(1<<i) << ", I_h = " << I0/omega
0130               << " = " << std::hexfloat << I0/omega << ", absolute error estimate = "
0131               << std::defaultfloat << std::scientific << abs(I0-I1)  << std::endl;
0132 }
0133 #endif
0134 
0135 
0136 template<class Real>
0137 std::pair<Real, Real> ooura_cos_node_and_weight(long n, Real h, Real alpha)
0138 {
0139     using std::expm1;
0140     using std::exp;
0141     using std::abs;
0142     using boost::math::constants::pi;
0143 
0144     Real x = h*(n-Real(1)/Real(2));
0145     auto p = ooura_eta(x, alpha);
0146     auto eta = p.first;
0147     auto eta_prime = p.second;
0148     Real expm1_meta = expm1(-eta);
0149     Real exp_meta = exp(-eta);
0150     Real node = pi<Real>()*(Real(1)/Real(2)-n)/expm1_meta;
0151 
0152     Real phi_prime = -(expm1_meta + x*exp_meta*eta_prime)/(expm1_meta*expm1_meta);
0153 
0154     // Takuya Ooura and Masatake Mori,
0155     // Journal of Computational and Applied Mathematics, 112 (1999) 229-241.
0156     // A robust double exponential formula for Fourier-type integrals.
0157     // Equation 4.6
0158     Real s = pi<Real>();
0159     Real arg;
0160     if (eta < -1) {
0161         arg = -(n-Real(1)/Real(2))/expm1_meta;
0162         s *= boost::math::cos_pi(arg);
0163     }
0164     else {
0165         arg = -(n-Real(1)/Real(2))*exp_meta/expm1_meta;
0166         s *= boost::math::sin_pi(arg);
0167         if (n&1) {
0168             s *= -1;
0169         }
0170     }
0171 
0172     Real weight = s*phi_prime;
0173     return {node, weight};
0174 }
0175 
0176 
0177 template<class Real>
0178 class ooura_fourier_sin_detail {
0179 public:
0180     ooura_fourier_sin_detail(const Real relative_error_goal, size_t levels) {
0181 #ifdef BOOST_MATH_INSTRUMENT_OOURA
0182       std::cout << "ooura_fourier_sin with relative error goal " << relative_error_goal 
0183         << " & " << levels << " levels." << std::endl;
0184 #endif // BOOST_MATH_INSTRUMENT_OOURA
0185         if (relative_error_goal < std::numeric_limits<Real>::epsilon() * 2) {
0186             throw std::domain_error("The relative error goal cannot be smaller than the unit roundoff.");
0187         }
0188         using std::abs;
0189         requested_levels_ = levels;
0190         starting_level_ = 0;
0191         rel_err_goal_ = relative_error_goal;
0192         big_nodes_.reserve(levels);
0193         bweights_.reserve(levels);
0194         little_nodes_.reserve(levels);
0195         lweights_.reserve(levels);
0196 
0197         for (size_t i = 0; i < levels; ++i) {
0198             if (std::is_same<Real, float>::value) {
0199                 add_level<double>(i);
0200             }
0201             else if (std::is_same<Real, double>::value) {
0202                 add_level<long double>(i);
0203             }
0204             else {
0205                 add_level<Real>(i);
0206             }
0207         }
0208     }
0209 
0210     std::vector<std::vector<Real>> const & big_nodes() const {
0211         return big_nodes_;
0212     }
0213 
0214     std::vector<std::vector<Real>> const & weights_for_big_nodes() const {
0215         return bweights_;
0216     }
0217 
0218     std::vector<std::vector<Real>> const & little_nodes() const {
0219         return little_nodes_;
0220     }
0221 
0222     std::vector<std::vector<Real>> const & weights_for_little_nodes() const {
0223         return lweights_;
0224     }
0225 
0226     template<class F>
0227     std::pair<Real,Real> integrate(F const & f, Real omega) {
0228         using std::abs;
0229         using std::max;
0230         using boost::math::constants::pi;
0231 
0232         if (omega == 0) {
0233             return {Real(0), Real(0)};
0234         }
0235         if (omega < 0) {
0236             auto p = this->integrate(f, -omega);
0237             return {-p.first, p.second};
0238         }
0239 
0240         Real I1 = std::numeric_limits<Real>::quiet_NaN();
0241         Real relative_error_estimate = std::numeric_limits<Real>::quiet_NaN();
0242         // As we compute integrals, we learn about their structure.
0243         // Assuming we compute f(t)sin(wt) for many different omega, this gives some
0244         // a posteriori ability to choose a refinement level that is roughly appropriate.
0245         size_t i = starting_level_;
0246         do {
0247             Real I0 = estimate_integral(f, omega, i);
0248 #ifdef BOOST_MATH_INSTRUMENT_OOURA
0249             print_ooura_estimate(i, I0, I1, omega);
0250 #endif
0251             Real absolute_error_estimate = abs(I0-I1);
0252             Real scale = (max)(abs(I0), abs(I1));
0253             if (!isnan(I1) && absolute_error_estimate <= rel_err_goal_*scale) {
0254                 starting_level_ = (max)(long(i) - 1, long(0));
0255                 return {I0/omega, absolute_error_estimate/scale};
0256             }
0257             I1 = I0;
0258         } while(++i < big_nodes_.size());
0259 
0260         // We've used up all our precomputed levels.
0261         // Now we need to add more.
0262         // It might seems reasonable to just keep adding levels indefinitely, if that's what the user wants.
0263         // But in fact the nodes and weights just merge into each other and the error gets worse after a certain number.
0264         // This value for max_additional_levels was chosen by observation of a slowly converging oscillatory integral:
0265         // f(x) := cos(7cos(x))sin(x)/x
0266         size_t max_additional_levels = 4;
0267         while (big_nodes_.size() < requested_levels_ + max_additional_levels) {
0268             size_t ii = big_nodes_.size();
0269             if (std::is_same<Real, float>::value) {
0270                 add_level<double>(ii);
0271             }
0272             else if (std::is_same<Real, double>::value) {
0273                 add_level<long double>(ii);
0274             }
0275             else {
0276                 add_level<Real>(ii);
0277             }
0278             Real I0 = estimate_integral(f, omega, ii);
0279             Real absolute_error_estimate = abs(I0-I1);
0280             Real scale = (max)(abs(I0), abs(I1));
0281 #ifdef BOOST_MATH_INSTRUMENT_OOURA
0282             print_ooura_estimate(ii, I0, I1, omega);
0283 #endif
0284             if (absolute_error_estimate <= rel_err_goal_*scale) {
0285                 starting_level_ = (max)(long(ii) - 1, long(0));
0286                 return {I0/omega, absolute_error_estimate/scale};
0287             }
0288             I1 = I0;
0289             ++ii;
0290         }
0291 
0292         starting_level_ = static_cast<long>(big_nodes_.size() - 2);
0293         return {I1/omega, relative_error_estimate};
0294     }
0295 
0296 private:
0297 
0298     template<class PreciseReal>
0299     void add_level(size_t i) {
0300         using std::abs;
0301         size_t current_num_levels = big_nodes_.size();
0302         Real unit_roundoff = std::numeric_limits<Real>::epsilon()/2;
0303         // h0 = 1. Then all further levels have h_i = 1/2^i.
0304         // Since the nodes don't nest, we could conceivably divide h by (say) 1.5, or 3.
0305         // It's not clear how much benefit (or loss) would be obtained from this.
0306         PreciseReal h = PreciseReal(1)/PreciseReal(1<<i);
0307 
0308         std::vector<Real> bnode_row;
0309         std::vector<Real> bweight_row;
0310 
0311         // This is a pretty good estimate for how many elements will be placed in the vector:
0312         bnode_row.reserve((static_cast<size_t>(1)<<i)*sizeof(Real));
0313         bweight_row.reserve((static_cast<size_t>(1)<<i)*sizeof(Real));
0314 
0315         std::vector<Real> lnode_row;
0316         std::vector<Real> lweight_row;
0317 
0318         lnode_row.reserve((static_cast<size_t>(1)<<i)*sizeof(Real));
0319         lweight_row.reserve((static_cast<size_t>(1)<<i)*sizeof(Real));
0320 
0321         Real max_weight = 1;
0322         auto alpha = calculate_ooura_alpha(h);
0323         long n = 0;
0324         Real w;
0325         do {
0326             auto precise_nw = ooura_sin_node_and_weight(n, h, alpha);
0327             Real node = static_cast<Real>(precise_nw.first);
0328             Real weight = static_cast<Real>(precise_nw.second);
0329             w = weight;
0330             if (bnode_row.size() == bnode_row.capacity()) {
0331                 bnode_row.reserve(2*bnode_row.size());
0332                 bweight_row.reserve(2*bnode_row.size());
0333             }
0334 
0335             bnode_row.push_back(node);
0336             bweight_row.push_back(weight);
0337             if (abs(weight) > max_weight) {
0338                 max_weight = abs(weight);
0339             }
0340             ++n;
0341             // f(t)->0 as t->infty, which is why the weights are computed up to the unit roundoff.
0342         } while(abs(w) > unit_roundoff*max_weight);
0343 
0344         // This class tends to consume a lot of memory; shrink the vectors back down to size:
0345         bnode_row.shrink_to_fit();
0346         bweight_row.shrink_to_fit();
0347         // Why we are splitting the nodes into regimes where t_n >> 1 and t_n << 1?
0348         // It will create the opportunity to sensibly truncate the quadrature sum to significant terms.
0349         n = -1;
0350         do {
0351             auto precise_nw = ooura_sin_node_and_weight(n, h, alpha);
0352             Real node = static_cast<Real>(precise_nw.first);
0353             if (node <= 0) {
0354                 break;
0355             }
0356             Real weight = static_cast<Real>(precise_nw.second);
0357             w = weight;
0358             using std::isnan;
0359             if (isnan(node)) {
0360                 // This occurs at n = -11 in quad precision:
0361                 break;
0362             }
0363             if (lnode_row.size() > 0) {
0364                 if (lnode_row[lnode_row.size()-1] == node) {
0365                     // The nodes have fused into each other:
0366                     break;
0367                 }
0368             }
0369             if (lnode_row.size() == lnode_row.capacity()) {
0370                 lnode_row.reserve(2*lnode_row.size());
0371                 lweight_row.reserve(2*lnode_row.size());
0372             }
0373             lnode_row.push_back(node);
0374             lweight_row.push_back(weight);
0375             if (abs(weight) > max_weight) {
0376                 max_weight = abs(weight);
0377             }
0378             --n;
0379             // f(t)->infty is possible as t->0, hence compute up to the min.
0380         } while(abs(w) > (std::numeric_limits<Real>::min)()*max_weight);
0381 
0382         lnode_row.shrink_to_fit();
0383         lweight_row.shrink_to_fit();
0384 
0385         #ifdef BOOST_HAS_THREADS
0386         // std::scoped_lock once C++17 is more common?
0387         std::lock_guard<std::mutex> lock(node_weight_mutex_);
0388         #endif 
0389         // Another thread might have already finished this calculation and appended it to the nodes/weights:
0390         if (current_num_levels == big_nodes_.size()) {
0391             big_nodes_.push_back(bnode_row);
0392             bweights_.push_back(bweight_row);
0393 
0394             little_nodes_.push_back(lnode_row);
0395             lweights_.push_back(lweight_row);
0396         }
0397     }
0398 
0399     template<class F>
0400     Real estimate_integral(F const & f, Real omega, size_t i) {
0401         // Because so few function evaluations are required to get high accuracy on the integrals in the tests,
0402         // Kahan summation doesn't really help.
0403         //auto cond = boost::math::tools::summation_condition_number<Real, true>(0);
0404         Real I0 = 0;
0405         auto const & b_nodes = big_nodes_[i];
0406         auto const & b_weights = bweights_[i];
0407         // Will benchmark if this is helpful:
0408         Real inv_omega = 1/omega;
0409         for(size_t j = 0 ; j < b_nodes.size(); ++j) {
0410             I0 += f(b_nodes[j]*inv_omega)*b_weights[j];
0411         }
0412 
0413         auto const & l_nodes = little_nodes_[i];
0414         auto const & l_weights = lweights_[i];
0415         // If f decays rapidly as |t|->infty, not all of these calls are necessary.
0416         for (size_t j = 0; j < l_nodes.size(); ++j) {
0417             I0 += f(l_nodes[j]*inv_omega)*l_weights[j];
0418         }
0419         return I0;
0420     }
0421 
0422     #ifdef BOOST_HAS_THREADS
0423     std::mutex node_weight_mutex_;
0424     #endif
0425     // Nodes for n >= 0, giving t_n = pi*phi(nh)/h. Generally t_n >> 1.
0426     std::vector<std::vector<Real>> big_nodes_;
0427     // The term bweights_ will indicate that these are weights corresponding
0428     // to the big nodes:
0429     std::vector<std::vector<Real>> bweights_;
0430 
0431     // Nodes for n < 0: Generally t_n << 1, and an invariant is that t_n > 0.
0432     std::vector<std::vector<Real>> little_nodes_;
0433     std::vector<std::vector<Real>> lweights_;
0434     Real rel_err_goal_;
0435 
0436     #ifdef BOOST_HAS_THREADS
0437     std::atomic<long> starting_level_{};
0438     #else
0439     long starting_level_;
0440     #endif
0441     size_t requested_levels_;
0442 };
0443 
0444 template<class Real>
0445 class ooura_fourier_cos_detail {
0446 public:
0447     ooura_fourier_cos_detail(const Real relative_error_goal, size_t levels) {
0448 #ifdef BOOST_MATH_INSTRUMENT_OOURA
0449       std::cout << "ooura_fourier_cos with relative error goal " << relative_error_goal
0450         << " & " << levels << " levels." << std::endl;
0451       std::cout << "epsilon for type = " << std::numeric_limits<Real>::epsilon() << std::endl;
0452 #endif // BOOST_MATH_INSTRUMENT_OOURA
0453         if (relative_error_goal < std::numeric_limits<Real>::epsilon() * 2) {
0454             throw std::domain_error("The relative error goal cannot be smaller than the unit roundoff!");
0455         }
0456 
0457         using std::abs;
0458         requested_levels_ = levels;
0459         starting_level_ = 0;
0460         rel_err_goal_ = relative_error_goal;
0461         big_nodes_.reserve(levels);
0462         bweights_.reserve(levels);
0463         little_nodes_.reserve(levels);
0464         lweights_.reserve(levels);
0465 
0466         for (size_t i = 0; i < levels; ++i) {
0467             if (std::is_same<Real, float>::value) {
0468                 add_level<double>(i);
0469             }
0470             else if (std::is_same<Real, double>::value) {
0471                 add_level<long double>(i);
0472             }
0473             else {
0474                 add_level<Real>(i);
0475             }
0476         }
0477 
0478     }
0479 
0480     template<class F>
0481     std::pair<Real,Real> integrate(F const & f, Real omega) {
0482         using std::abs;
0483         using std::max;
0484         using boost::math::constants::pi;
0485 
0486         if (omega == 0) {
0487             throw std::domain_error("At omega = 0, the integral is not oscillatory. The user must choose an appropriate method for this case.\n");
0488         }
0489 
0490         if (omega < 0) {
0491             return this->integrate(f, -omega);
0492         }
0493 
0494         Real I1 = std::numeric_limits<Real>::quiet_NaN();
0495         Real absolute_error_estimate = std::numeric_limits<Real>::quiet_NaN();
0496         Real scale = std::numeric_limits<Real>::quiet_NaN();
0497         size_t i = starting_level_;
0498         do {
0499             Real I0 = estimate_integral(f, omega, i);
0500 #ifdef BOOST_MATH_INSTRUMENT_OOURA
0501             print_ooura_estimate(i, I0, I1, omega);
0502 #endif
0503             absolute_error_estimate = abs(I0-I1);
0504             scale = (max)(abs(I0), abs(I1));
0505             if (!isnan(I1) && absolute_error_estimate <= rel_err_goal_*scale) {
0506                 starting_level_ = (max)(long(i) - 1, long(0));
0507                 return {I0/omega, absolute_error_estimate/scale};
0508             }
0509             I1 = I0;
0510         } while(++i < big_nodes_.size());
0511 
0512         size_t max_additional_levels = 4;
0513         while (big_nodes_.size() < requested_levels_ + max_additional_levels) {
0514             size_t ii = big_nodes_.size();
0515             if (std::is_same<Real, float>::value) {
0516                 add_level<double>(ii);
0517             }
0518             else if (std::is_same<Real, double>::value) {
0519                 add_level<long double>(ii);
0520             }
0521             else {
0522                 add_level<Real>(ii);
0523             }
0524             Real I0 = estimate_integral(f, omega, ii);
0525 #ifdef BOOST_MATH_INSTRUMENT_OOURA
0526             print_ooura_estimate(ii, I0, I1, omega);
0527 #endif
0528             absolute_error_estimate = abs(I0-I1);
0529             scale = (max)(abs(I0), abs(I1));
0530             if (absolute_error_estimate <= rel_err_goal_*scale) {
0531                 starting_level_ = (max)(long(ii) - 1, long(0));
0532                 return {I0/omega, absolute_error_estimate/scale};
0533             }
0534             I1 = I0;
0535             ++ii;
0536         }
0537 
0538         starting_level_ = static_cast<long>(big_nodes_.size() - 2);
0539         return {I1/omega, absolute_error_estimate/scale};
0540     }
0541 
0542 private:
0543 
0544     template<class PreciseReal>
0545     void add_level(size_t i) {
0546         using std::abs;
0547         size_t current_num_levels = big_nodes_.size();
0548         Real unit_roundoff = std::numeric_limits<Real>::epsilon()/2;
0549         PreciseReal h = PreciseReal(1)/PreciseReal(1<<i);
0550 
0551         std::vector<Real> bnode_row;
0552         std::vector<Real> bweight_row;
0553         bnode_row.reserve((static_cast<size_t>(1)<<i)*sizeof(Real));
0554         bweight_row.reserve((static_cast<size_t>(1)<<i)*sizeof(Real));
0555 
0556         std::vector<Real> lnode_row;
0557         std::vector<Real> lweight_row;
0558 
0559         lnode_row.reserve((static_cast<size_t>(1)<<i)*sizeof(Real));
0560         lweight_row.reserve((static_cast<size_t>(1)<<i)*sizeof(Real));
0561 
0562         Real max_weight = 1;
0563         auto alpha = calculate_ooura_alpha(h);
0564         long n = 0;
0565         Real w;
0566         do {
0567             auto precise_nw = ooura_cos_node_and_weight(n, h, alpha);
0568             Real node = static_cast<Real>(precise_nw.first);
0569             Real weight = static_cast<Real>(precise_nw.second);
0570             w = weight;
0571             if (bnode_row.size() == bnode_row.capacity()) {
0572                 bnode_row.reserve(2*bnode_row.size());
0573                 bweight_row.reserve(2*bnode_row.size());
0574             }
0575 
0576             bnode_row.push_back(node);
0577             bweight_row.push_back(weight);
0578             if (abs(weight) > max_weight) {
0579                 max_weight = abs(weight);
0580             }
0581             ++n;
0582             // f(t)->0 as t->infty, which is why the weights are computed up to the unit roundoff.
0583         } while(abs(w) > unit_roundoff*max_weight);
0584 
0585         bnode_row.shrink_to_fit();
0586         bweight_row.shrink_to_fit();
0587         n = -1;
0588         do {
0589             auto precise_nw = ooura_cos_node_and_weight(n, h, alpha);
0590             Real node = static_cast<Real>(precise_nw.first);
0591             // The function cannot be singular at zero,
0592             // so zero is not a unreasonable node,
0593             // unlike in the case of the Fourier Sine.
0594             // Hence only break if the node is negative.
0595             if (node < 0) {
0596                 break;
0597             }
0598             Real weight = static_cast<Real>(precise_nw.second);
0599             w = weight;
0600             if (lnode_row.size() > 0) {
0601                 if (lnode_row.back() == node) {
0602                     // The nodes have fused into each other:
0603                     break;
0604                 }
0605             }
0606             if (lnode_row.size() == lnode_row.capacity()) {
0607                 lnode_row.reserve(2*lnode_row.size());
0608                 lweight_row.reserve(2*lnode_row.size());
0609             }
0610 
0611             lnode_row.push_back(node);
0612             lweight_row.push_back(weight);
0613             if (abs(weight) > max_weight) {
0614                 max_weight = abs(weight);
0615             }
0616             --n;
0617         } while(abs(w) > (std::numeric_limits<Real>::min)()*max_weight);
0618 
0619         lnode_row.shrink_to_fit();
0620         lweight_row.shrink_to_fit();
0621 
0622         #ifdef BOOST_HAS_THREADS
0623         std::lock_guard<std::mutex> lock(node_weight_mutex_);
0624         #endif
0625 
0626         // Another thread might have already finished this calculation and appended it to the nodes/weights:
0627         if (current_num_levels == big_nodes_.size()) {
0628             big_nodes_.push_back(bnode_row);
0629             bweights_.push_back(bweight_row);
0630 
0631             little_nodes_.push_back(lnode_row);
0632             lweights_.push_back(lweight_row);
0633         }
0634     }
0635 
0636     template<class F>
0637     Real estimate_integral(F const & f, Real omega, size_t i) {
0638         Real I0 = 0;
0639         auto const & b_nodes = big_nodes_[i];
0640         auto const & b_weights = bweights_[i];
0641         Real inv_omega = 1/omega;
0642         for(size_t j = 0 ; j < b_nodes.size(); ++j) {
0643             I0 += f(b_nodes[j]*inv_omega)*b_weights[j];
0644         }
0645 
0646         auto const & l_nodes = little_nodes_[i];
0647         auto const & l_weights = lweights_[i];
0648         for (size_t j = 0; j < l_nodes.size(); ++j) {
0649             I0 += f(l_nodes[j]*inv_omega)*l_weights[j];
0650         }
0651         return I0;
0652     }
0653 
0654     #ifdef BOOST_HAS_THREADS
0655     std::mutex node_weight_mutex_;
0656     #endif 
0657 
0658     std::vector<std::vector<Real>> big_nodes_;
0659     std::vector<std::vector<Real>> bweights_;
0660 
0661     std::vector<std::vector<Real>> little_nodes_;
0662     std::vector<std::vector<Real>> lweights_;
0663     Real rel_err_goal_;
0664 
0665     #ifdef BOOST_HAS_THREADS
0666     std::atomic<long> starting_level_{};
0667     #else
0668     long starting_level_;
0669     #endif
0670 
0671     size_t requested_levels_;
0672 };
0673 
0674 
0675 }}}}
0676 #endif