boost/math/quadrature/detail/ooura_fourier_integrals_detail.hpp

   1 // Copyright Nick Thompson, 2019
   2 // Use, modification and distribution are subject to the
   3 // Boost Software License, Version 1.0.
   4 // (See accompanying file LICENSE_1_0.txt
   5 // or copy at http://www.boost.org/LICENSE_1_0.txt)
   6 #ifndef BOOST_MATH_QUADRATURE_DETAIL_OOURA_FOURIER_INTEGRALS_DETAIL_HPP
   7 #define BOOST_MATH_QUADRATURE_DETAIL_OOURA_FOURIER_INTEGRALS_DETAIL_HPP
   8 #include <utility> // for std::pair.
   9 #include <mutex>
  10 #include <atomic>
  11 #include <vector>
  12 #include <iostream>
  13 #include <boost/math/special_functions/expm1.hpp>
  14 #include <boost/math/special_functions/sin_pi.hpp>
  15 #include <boost/math/special_functions/cos_pi.hpp>
  16 #include <boost/math/constants/constants.hpp>
  17
  18 namespace boost { namespace math { namespace quadrature { namespace detail {
  19
  20 // Ooura and Mori, A robust double exponential formula for Fourier-type integrals,
  21 // eta is the argument to the exponential in equation 3.3:
  22 template<class Real>
  23 std::pair<Real, Real> ooura_eta(Real x, Real alpha) {
  24     using std::expm1;
  25     using std::exp;
  26     using std::abs;
  27     Real expx = exp(x);
  28     Real eta_prime = 2 + alpha/expx + expx/4;
  29     Real eta;
  30     // This is the fast branch:
  31     if (abs(x) > 0.125) {
  32         eta = 2*x - alpha*(1/expx - 1) + (expx - 1)/4;
  33     }
  34     else {// this is the slow branch using expm1 for small x:
  35         eta = 2*x - alpha*expm1(-x) + expm1(x)/4;
  36     }
  37     return {eta, eta_prime};
  38 }
  39
  40 // Ooura and Mori, A robust double exponential formula for Fourier-type integrals,
  41 // equation 3.6:
  42 template<class Real>
  43 Real calculate_ooura_alpha(Real h)
  44 {
  45     using boost::math::constants::pi;
  46     using std::log1p;
  47     using std::sqrt;
  48     Real x = sqrt(16 + 4*log1p(pi<Real>()/h)/h);
  49     return 1/x;
  50 }
  51
  52 template<class Real>
  53 std::pair<Real, Real> ooura_sin_node_and_weight(long n, Real h, Real alpha)
  54 {
  55     using std::expm1;
  56     using std::exp;
  57     using std::abs;
  58     using boost::math::constants::pi;
  59     using std::isnan;
  60
  61     if (n == 0) {
  62         // Equation 44 of https://arxiv.org/pdf/0911.4796.pdf
  63         // Fourier Transform of the Stretched Exponential Function: Analytic Error Bounds,
  64         // Double Exponential Transform, and Open-Source Implementation,
  65         // Joachim Wuttke,
  66         // The C library libkww provides functions to compute the Kohlrausch-Williams-Watts function,
  67         // the Laplace-Fourier transform of the stretched (or compressed) exponential function exp(-t^beta)
  68         // for exponent beta between 0.1 and 1.9 with sixteen decimal digits accuracy.
  69
  70         Real eta_prime_0 = Real(2) + alpha + Real(1)/Real(4);
  71         Real node = pi<Real>()/(eta_prime_0*h);
  72         Real weight = pi<Real>()*boost::math::sin_pi(1/(eta_prime_0*h));
  73         Real eta_dbl_prime = -alpha + Real(1)/Real(4);
  74         Real phi_prime_0 = (1 - eta_dbl_prime/(eta_prime_0*eta_prime_0))/2;
  75         weight *= phi_prime_0;
  76         return {node, weight};
  77     }
  78     Real x = n*h;
  79     auto p = ooura_eta(x, alpha);
  80     auto eta = p.first;
  81     auto eta_prime = p.second;
  82
  83     Real expm1_meta = expm1(-eta);
  84     Real exp_meta = exp(-eta);
  85     Real node = -n*pi<Real>()/expm1_meta;
  86
  87
  88     // I have verified that this is not a significant source of inaccuracy in the weight computation:
  89     Real phi_prime = -(expm1_meta + x*exp_meta*eta_prime)/(expm1_meta*expm1_meta);
  90
  91     // The main source of inaccuracy is in computation of sin_pi.
  92     // But I've agonized over this, and I think it's as good as it can get:
  93     Real s = pi<Real>();
  94     Real arg;
  95     if(eta > 1) {
  96         arg = n/( 1/exp_meta - 1 );
  97         s *= boost::math::sin_pi(arg);
  98         if (n&1) {
  99             s *= -1;
 100         }
 101     }
 102     else if (eta < -1) {
 103         arg = n/(1-exp_meta);
 104         s *= boost::math::sin_pi(arg);
 105     }
 106     else {
 107         arg = -n*exp_meta/expm1_meta;
 108         s *= boost::math::sin_pi(arg);
 109         if (n&1) {
 110             s *= -1;
 111         }
 112     }
 113
 114     Real weight = s*phi_prime;
 115     return {node, weight};
 116 }
 117
 118 #ifdef BOOST_MATH_INSTRUMENT_OOURA
 119 template<class Real>
 120 void print_ooura_estimate(size_t i, Real I0, Real I1, Real omega) {
 121     using std::abs;
 122     std::cout << std::defaultfloat
 123               << std::setprecision(std::numeric_limits<Real>::digits10)
 124               << std::fixed;
 125     std::cout << "h = " << Real(1)/Real(1<<i) << ", I_h = " << I0/omega
 126               << " = " << std::hexfloat << I0/omega << ", absolute error estimate = "
 127               << std::defaultfloat << std::scientific << abs(I0-I1)  << std::endl;
 128 }
 129 #endif
 130
 131
 132 template<class Real>
 133 std::pair<Real, Real> ooura_cos_node_and_weight(long n, Real h, Real alpha)
 134 {
 135     using std::expm1;
 136     using std::exp;
 137     using std::abs;
 138     using boost::math::constants::pi;
 139
 140     Real x = h*(n-Real(1)/Real(2));
 141     auto p = ooura_eta(x, alpha);
 142     auto eta = p.first;
 143     auto eta_prime = p.second;
 144     Real expm1_meta = expm1(-eta);
 145     Real exp_meta = exp(-eta);
 146     Real node = pi<Real>()*(Real(1)/Real(2)-n)/expm1_meta;
 147
 148     Real phi_prime = -(expm1_meta + x*exp_meta*eta_prime)/(expm1_meta*expm1_meta);
 149
 150     // Takuya Ooura and Masatake Mori,
 151     // Journal of Computational and Applied Mathematics, 112 (1999) 229-241.
 152     // A robust double exponential formula for Fourier-type integrals.
 153     // Equation 4.6
 154     Real s = pi<Real>();
 155     Real arg;
 156     if (eta < -1) {
 157         arg = -(n-Real(1)/Real(2))/expm1_meta;
 158         s *= boost::math::cos_pi(arg);
 159     }
 160     else {
 161         arg = -(n-Real(1)/Real(2))*exp_meta/expm1_meta;
 162         s *= boost::math::sin_pi(arg);
 163         if (n&1) {
 164             s *= -1;
 165         }
 166     }
 167
 168     Real weight = s*phi_prime;
 169     return {node, weight};
 170 }
 171
 172
 173 template<class Real>
 174 class ooura_fourier_sin_detail {
 175 public:
 176     ooura_fourier_sin_detail(const Real relative_error_goal, size_t levels) {
 177 #ifdef BOOST_MATH_INSTRUMENT_OOURA
 178       std::cout << "ooura_fourier_sin with relative error goal " << relative_error_goal
 179         << " & " << levels << " levels." << std::endl;
 180 #endif // BOOST_MATH_INSTRUMENT_OOURA
 181         if (relative_error_goal < std::numeric_limits<Real>::epsilon() * 2) {
 182             throw std::domain_error("The relative error goal cannot be smaller than the unit roundoff.");
 183         }
 184         using std::abs;
 185         requested_levels_ = levels;
 186         starting_level_ = 0;
 187         rel_err_goal_ = relative_error_goal;
 188         big_nodes_.reserve(levels);
 189         bweights_.reserve(levels);
 190         little_nodes_.reserve(levels);
 191         lweights_.reserve(levels);
 192
 193         for (size_t i = 0; i < levels; ++i) {
 194             if (std::is_same<Real, float>::value) {
 195                 add_level<double>(i);
 196             }
 197             else if (std::is_same<Real, double>::value) {
 198                 add_level<long double>(i);
 199             }
 200             else {
 201                 add_level<Real>(i);
 202             }
 203         }
 204     }
 205
 206     std::vector<std::vector<Real>> const & big_nodes() const {
 207         return big_nodes_;
 208     }
 209
 210     std::vector<std::vector<Real>> const & weights_for_big_nodes() const {
 211         return bweights_;
 212     }
 213
 214     std::vector<std::vector<Real>> const & little_nodes() const {
 215         return little_nodes_;
 216     }
 217
 218     std::vector<std::vector<Real>> const & weights_for_little_nodes() const {
 219         return lweights_;
 220     }
 221
 222     template<class F>
 223     std::pair<Real,Real> integrate(F const & f, Real omega) {
 224         using std::abs;
 225         using std::max;
 226         using boost::math::constants::pi;
 227
 228         if (omega == 0) {
 229             return {Real(0), Real(0)};
 230         }
 231         if (omega < 0) {
 232             auto p = this->integrate(f, -omega);
 233             return {-p.first, p.second};
 234         }
 235
 236         Real I1 = std::numeric_limits<Real>::quiet_NaN();
 237         Real relative_error_estimate = std::numeric_limits<Real>::quiet_NaN();
 238         // As we compute integrals, we learn about their structure.
 239         // Assuming we compute f(t)sin(wt) for many different omega, this gives some
 240         // a posteriori ability to choose a refinement level that is roughly appropriate.
 241         size_t i = starting_level_;
 242         do {
 243             Real I0 = estimate_integral(f, omega, i);
 244 #ifdef BOOST_MATH_INSTRUMENT_OOURA
 245             print_ooura_estimate(i, I0, I1, omega);
 246 #endif
 247             Real absolute_error_estimate = abs(I0-I1);
 248             Real scale = (max)(abs(I0), abs(I1));
 249             if (!isnan(I1) && absolute_error_estimate <= rel_err_goal_*scale) {
 250                 starting_level_ = (max)(long(i) - 1, long(0));
 251                 return {I0/omega, absolute_error_estimate/scale};
 252             }
 253             I1 = I0;
 254         } while(++i < big_nodes_.size());
 255
 256         // We've used up all our precomputed levels.
 257         // Now we need to add more.
 258         // It might seems reasonable to just keep adding levels indefinitely, if that's what the user wants.
 259         // But in fact the nodes and weights just merge into each other and the error gets worse after a certain number.
 260         // This value for max_additional_levels was chosen by observation of a slowly converging oscillatory integral:
 261         // f(x) := cos(7cos(x))sin(x)/x
 262         size_t max_additional_levels = 4;
 263         while (big_nodes_.size() < requested_levels_ + max_additional_levels) {
 264             size_t ii = big_nodes_.size();
 265             if (std::is_same<Real, float>::value) {
 266                 add_level<double>(ii);
 267             }
 268             else if (std::is_same<Real, double>::value) {
 269                 add_level<long double>(ii);
 270             }
 271             else {
 272                 add_level<Real>(ii);
 273             }
 274             Real I0 = estimate_integral(f, omega, ii);
 275             Real absolute_error_estimate = abs(I0-I1);
 276             Real scale = (max)(abs(I0), abs(I1));
 277 #ifdef BOOST_MATH_INSTRUMENT_OOURA
 278             print_ooura_estimate(ii, I0, I1, omega);
 279 #endif
 280             if (absolute_error_estimate <= rel_err_goal_*scale) {
 281                 starting_level_ = (max)(long(ii) - 1, long(0));
 282                 return {I0/omega, absolute_error_estimate/scale};
 283             }
 284             I1 = I0;
 285             ++ii;
 286         }
 287
 288         starting_level_ = static_cast<long>(big_nodes_.size() - 2);
 289         return {I1/omega, relative_error_estimate};
 290     }
 291
 292 private:
 293
 294     template<class PreciseReal>
 295     void add_level(size_t i) {
 296         using std::abs;
 297         size_t current_num_levels = big_nodes_.size();
 298         Real unit_roundoff = std::numeric_limits<Real>::epsilon()/2;
 299         // h0 = 1. Then all further levels have h_i = 1/2^i.
 300         // Since the nodes don't nest, we could conceivably divide h by (say) 1.5, or 3.
 301         // It's not clear how much benefit (or loss) would be obtained from this.
 302         PreciseReal h = PreciseReal(1)/PreciseReal(1<<i);
 303
 304         std::vector<Real> bnode_row;
 305         std::vector<Real> bweight_row;
 306
 307         // This is a pretty good estimate for how many elements will be placed in the vector:
 308         bnode_row.reserve((static_cast<size_t>(1)<<i)*sizeof(Real));
 309         bweight_row.reserve((static_cast<size_t>(1)<<i)*sizeof(Real));
 310
 311         std::vector<Real> lnode_row;
 312         std::vector<Real> lweight_row;
 313
 314         lnode_row.reserve((static_cast<size_t>(1)<<i)*sizeof(Real));
 315         lweight_row.reserve((static_cast<size_t>(1)<<i)*sizeof(Real));
 316
 317         Real max_weight = 1;
 318         auto alpha = calculate_ooura_alpha(h);
 319         long n = 0;
 320         Real w;
 321         do {
 322             auto precise_nw = ooura_sin_node_and_weight(n, h, alpha);
 323             Real node = static_cast<Real>(precise_nw.first);
 324             Real weight = static_cast<Real>(precise_nw.second);
 325             w = weight;
 326             if (bnode_row.size() == bnode_row.capacity()) {
 327                 bnode_row.reserve(2*bnode_row.size());
 328                 bweight_row.reserve(2*bnode_row.size());
 329             }
 330
 331             bnode_row.push_back(node);
 332             bweight_row.push_back(weight);
 333             if (abs(weight) > max_weight) {
 334                 max_weight = abs(weight);
 335             }
 336             ++n;
 337             // f(t)->0 as t->infty, which is why the weights are computed up to the unit roundoff.
 338         } while(abs(w) > unit_roundoff*max_weight);
 339
 340         // This class tends to consume a lot of memory; shrink the vectors back down to size:
 341         bnode_row.shrink_to_fit();
 342         bweight_row.shrink_to_fit();
 343         // Why we are splitting the nodes into regimes where t_n >> 1 and t_n << 1?
 344         // It will create the opportunity to sensibly truncate the quadrature sum to significant terms.
 345         n = -1;
 346         do {
 347             auto precise_nw = ooura_sin_node_and_weight(n, h, alpha);
 348             Real node = static_cast<Real>(precise_nw.first);
 349             if (node <= 0) {
 350                 break;
 351             }
 352             Real weight = static_cast<Real>(precise_nw.second);
 353             w = weight;
 354             using std::isnan;
 355             if (isnan(node)) {
 356                 // This occurs at n = -11 in quad precision:
 357                 break;
 358             }
 359             if (lnode_row.size() > 0) {
 360                 if (lnode_row[lnode_row.size()-1] == node) {
 361                     // The nodes have fused into each other:
 362                     break;
 363                 }
 364             }
 365             if (lnode_row.size() == lnode_row.capacity()) {
 366                 lnode_row.reserve(2*lnode_row.size());
 367                 lweight_row.reserve(2*lnode_row.size());
 368             }
 369             lnode_row.push_back(node);
 370             lweight_row.push_back(weight);
 371             if (abs(weight) > max_weight) {
 372                 max_weight = abs(weight);
 373             }
 374             --n;
 375             // f(t)->infty is possible as t->0, hence compute up to the min.
 376         } while(abs(w) > (std::numeric_limits<Real>::min)()*max_weight);
 377
 378         lnode_row.shrink_to_fit();
 379         lweight_row.shrink_to_fit();
 380
 381         // std::scoped_lock once C++17 is more common?
 382         std::lock_guard<std::mutex> lock(node_weight_mutex_);
 383         // Another thread might have already finished this calculation and appended it to the nodes/weights:
 384         if (current_num_levels == big_nodes_.size()) {
 385             big_nodes_.push_back(bnode_row);
 386             bweights_.push_back(bweight_row);
 387
 388             little_nodes_.push_back(lnode_row);
 389             lweights_.push_back(lweight_row);
 390         }
 391     }
 392
 393     template<class F>
 394     Real estimate_integral(F const & f, Real omega, size_t i) {
 395         // Because so few function evaluations are required to get high accuracy on the integrals in the tests,
 396         // Kahan summation doesn't really help.
 397         //auto cond = boost::math::tools::summation_condition_number<Real, true>(0);
 398         Real I0 = 0;
 399         auto const & b_nodes = big_nodes_[i];
 400         auto const & b_weights = bweights_[i];
 401         // Will benchmark if this is helpful:
 402         Real inv_omega = 1/omega;
 403         for(size_t j = 0 ; j < b_nodes.size(); ++j) {
 404             I0 += f(b_nodes[j]*inv_omega)*b_weights[j];
 405         }
 406
 407         auto const & l_nodes = little_nodes_[i];
 408         auto const & l_weights = lweights_[i];
 409         // If f decays rapidly as |t|->infty, not all of these calls are necessary.
 410         for (size_t j = 0; j < l_nodes.size(); ++j) {
 411             I0 += f(l_nodes[j]*inv_omega)*l_weights[j];
 412         }
 413         return I0;
 414     }
 415
 416     std::mutex node_weight_mutex_;
 417     // Nodes for n >= 0, giving t_n = pi*phi(nh)/h. Generally t_n >> 1.
 418     std::vector<std::vector<Real>> big_nodes_;
 419     // The term bweights_ will indicate that these are weights corresponding
 420     // to the big nodes:
 421     std::vector<std::vector<Real>> bweights_;
 422
 423     // Nodes for n < 0: Generally t_n << 1, and an invariant is that t_n > 0.
 424     std::vector<std::vector<Real>> little_nodes_;
 425     std::vector<std::vector<Real>> lweights_;
 426     Real rel_err_goal_;
 427     std::atomic<long> starting_level_;
 428     size_t requested_levels_;
 429 };
 430
 431 template<class Real>
 432 class ooura_fourier_cos_detail {
 433 public:
 434     ooura_fourier_cos_detail(const Real relative_error_goal, size_t levels) {
 435 #ifdef BOOST_MATH_INSTRUMENT_OOURA
 436       std::cout << "ooura_fourier_cos with relative error goal " << relative_error_goal
 437         << " & " << levels << " levels." << std::endl;
 438       std::cout << "epsilon for type = " << std::numeric_limits<Real>::epsilon() << std::endl;
 439 #endif // BOOST_MATH_INSTRUMENT_OOURA
 440         if (relative_error_goal < std::numeric_limits<Real>::epsilon() * 2) {
 441             throw std::domain_error("The relative error goal cannot be smaller than the unit roundoff!");
 442         }
 443
 444         using std::abs;
 445         requested_levels_ = levels;
 446         starting_level_ = 0;
 447         rel_err_goal_ = relative_error_goal;
 448         big_nodes_.reserve(levels);
 449         bweights_.reserve(levels);
 450         little_nodes_.reserve(levels);
 451         lweights_.reserve(levels);
 452
 453         for (size_t i = 0; i < levels; ++i) {
 454             if (std::is_same<Real, float>::value) {
 455                 add_level<double>(i);
 456             }
 457             else if (std::is_same<Real, double>::value) {
 458                 add_level<long double>(i);
 459             }
 460             else {
 461                 add_level<Real>(i);
 462             }
 463         }
 464
 465     }
 466
 467     template<class F>
 468     std::pair<Real,Real> integrate(F const & f, Real omega) {
 469         using std::abs;
 470         using std::max;
 471         using boost::math::constants::pi;
 472
 473         if (omega == 0) {
 474             throw std::domain_error("At omega = 0, the integral is not oscillatory. The user must choose an appropriate method for this case.\n");
 475         }
 476
 477         if (omega < 0) {
 478             return this->integrate(f, -omega);
 479         }
 480
 481         Real I1 = std::numeric_limits<Real>::quiet_NaN();
 482         Real absolute_error_estimate = std::numeric_limits<Real>::quiet_NaN();
 483         Real scale = std::numeric_limits<Real>::quiet_NaN();
 484         size_t i = starting_level_;
 485         do {
 486             Real I0 = estimate_integral(f, omega, i);
 487 #ifdef BOOST_MATH_INSTRUMENT_OOURA
 488             print_ooura_estimate(i, I0, I1, omega);
 489 #endif
 490             absolute_error_estimate = abs(I0-I1);
 491             scale = (max)(abs(I0), abs(I1));
 492             if (!isnan(I1) && absolute_error_estimate <= rel_err_goal_*scale) {
 493                 starting_level_ = (max)(long(i) - 1, long(0));
 494                 return {I0/omega, absolute_error_estimate/scale};
 495             }
 496             I1 = I0;
 497         } while(++i < big_nodes_.size());
 498
 499         size_t max_additional_levels = 4;
 500         while (big_nodes_.size() < requested_levels_ + max_additional_levels) {
 501             size_t ii = big_nodes_.size();
 502             if (std::is_same<Real, float>::value) {
 503                 add_level<double>(ii);
 504             }
 505             else if (std::is_same<Real, double>::value) {
 506                 add_level<long double>(ii);
 507             }
 508             else {
 509                 add_level<Real>(ii);
 510             }
 511             Real I0 = estimate_integral(f, omega, ii);
 512 #ifdef BOOST_MATH_INSTRUMENT_OOURA
 513             print_ooura_estimate(ii, I0, I1, omega);
 514 #endif
 515             absolute_error_estimate = abs(I0-I1);
 516             scale = (max)(abs(I0), abs(I1));
 517             if (absolute_error_estimate <= rel_err_goal_*scale) {
 518                 starting_level_ = (max)(long(ii) - 1, long(0));
 519                 return {I0/omega, absolute_error_estimate/scale};
 520             }
 521             I1 = I0;
 522             ++ii;
 523         }
 524
 525         starting_level_ = static_cast<long>(big_nodes_.size() - 2);
 526         return {I1/omega, absolute_error_estimate/scale};
 527     }
 528
 529 private:
 530
 531     template<class PreciseReal>
 532     void add_level(size_t i) {
 533         using std::abs;
 534         size_t current_num_levels = big_nodes_.size();
 535         Real unit_roundoff = std::numeric_limits<Real>::epsilon()/2;
 536         PreciseReal h = PreciseReal(1)/PreciseReal(1<<i);
 537
 538         std::vector<Real> bnode_row;
 539         std::vector<Real> bweight_row;
 540         bnode_row.reserve((static_cast<size_t>(1)<<i)*sizeof(Real));
 541         bweight_row.reserve((static_cast<size_t>(1)<<i)*sizeof(Real));
 542
 543         std::vector<Real> lnode_row;
 544         std::vector<Real> lweight_row;
 545
 546         lnode_row.reserve((static_cast<size_t>(1)<<i)*sizeof(Real));
 547         lweight_row.reserve((static_cast<size_t>(1)<<i)*sizeof(Real));
 548
 549         Real max_weight = 1;
 550         auto alpha = calculate_ooura_alpha(h);
 551         long n = 0;
 552         Real w;
 553         do {
 554             auto precise_nw = ooura_cos_node_and_weight(n, h, alpha);
 555             Real node = static_cast<Real>(precise_nw.first);
 556             Real weight = static_cast<Real>(precise_nw.second);
 557             w = weight;
 558             if (bnode_row.size() == bnode_row.capacity()) {
 559                 bnode_row.reserve(2*bnode_row.size());
 560                 bweight_row.reserve(2*bnode_row.size());
 561             }
 562
 563             bnode_row.push_back(node);
 564             bweight_row.push_back(weight);
 565             if (abs(weight) > max_weight) {
 566                 max_weight = abs(weight);
 567             }
 568             ++n;
 569             // f(t)->0 as t->infty, which is why the weights are computed up to the unit roundoff.
 570         } while(abs(w) > unit_roundoff*max_weight);
 571
 572         bnode_row.shrink_to_fit();
 573         bweight_row.shrink_to_fit();
 574         n = -1;
 575         do {
 576             auto precise_nw = ooura_cos_node_and_weight(n, h, alpha);
 577             Real node = static_cast<Real>(precise_nw.first);
 578             // The function cannot be singular at zero,
 579             // so zero is not a unreasonable node,
 580             // unlike in the case of the Fourier Sine.
 581             // Hence only break if the node is negative.
 582             if (node < 0) {
 583                 break;
 584             }
 585             Real weight = static_cast<Real>(precise_nw.second);
 586             w = weight;
 587             if (lnode_row.size() > 0) {
 588                 if (lnode_row.back() == node) {
 589                     // The nodes have fused into each other:
 590                     break;
 591                 }
 592             }
 593             if (lnode_row.size() == lnode_row.capacity()) {
 594                 lnode_row.reserve(2*lnode_row.size());
 595                 lweight_row.reserve(2*lnode_row.size());
 596             }
 597
 598             lnode_row.push_back(node);
 599             lweight_row.push_back(weight);
 600             if (abs(weight) > max_weight) {
 601                 max_weight = abs(weight);
 602             }
 603             --n;
 604         } while(abs(w) > (std::numeric_limits<Real>::min)()*max_weight);
 605
 606         lnode_row.shrink_to_fit();
 607         lweight_row.shrink_to_fit();
 608
 609         std::lock_guard<std::mutex> lock(node_weight_mutex_);
 610         // Another thread might have already finished this calculation and appended it to the nodes/weights:
 611         if (current_num_levels == big_nodes_.size()) {
 612             big_nodes_.push_back(bnode_row);
 613             bweights_.push_back(bweight_row);
 614
 615             little_nodes_.push_back(lnode_row);
 616             lweights_.push_back(lweight_row);
 617         }
 618     }
 619
 620     template<class F>
 621     Real estimate_integral(F const & f, Real omega, size_t i) {
 622         Real I0 = 0;
 623         auto const & b_nodes = big_nodes_[i];
 624         auto const & b_weights = bweights_[i];
 625         Real inv_omega = 1/omega;
 626         for(size_t j = 0 ; j < b_nodes.size(); ++j) {
 627             I0 += f(b_nodes[j]*inv_omega)*b_weights[j];
 628         }
 629
 630         auto const & l_nodes = little_nodes_[i];
 631         auto const & l_weights = lweights_[i];
 632         for (size_t j = 0; j < l_nodes.size(); ++j) {
 633             I0 += f(l_nodes[j]*inv_omega)*l_weights[j];
 634         }
 635         return I0;
 636     }
 637
 638     std::mutex node_weight_mutex_;
 639     std::vector<std::vector<Real>> big_nodes_;
 640     std::vector<std::vector<Real>> bweights_;
 641
 642     std::vector<std::vector<Real>> little_nodes_;
 643     std::vector<std::vector<Real>> lweights_;
 644     Real rel_err_goal_;
 645     std::atomic<long> starting_level_;
 646     size_t requested_levels_;
 647 };
 648
 649
 650 }}}}
 651 #endif