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>
ooura_eta(Real x,Real alpha)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>
calculate_ooura_alpha(Real h)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>
ooura_sin_node_and_weight(long n,Real h,Real alpha)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>
print_ooura_estimate(size_t i,Real I0,Real I1,Real omega)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>
ooura_cos_node_and_weight(long n,Real h,Real alpha)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:
ooura_fourier_sin_detail(const Real relative_error_goal,size_t levels)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
big_nodes() const206 std::vector<std::vector<Real>> const & big_nodes() const {
207 return big_nodes_;
208 }
209
weights_for_big_nodes() const210 std::vector<std::vector<Real>> const & weights_for_big_nodes() const {
211 return bweights_;
212 }
213
little_nodes() const214 std::vector<std::vector<Real>> const & little_nodes() const {
215 return little_nodes_;
216 }
217
weights_for_little_nodes() const218 std::vector<std::vector<Real>> const & weights_for_little_nodes() const {
219 return lweights_;
220 }
221
222 template<class F>
integrate(F const & f,Real omega)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>
add_level(size_t i)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>
estimate_integral(F const & f,Real omega,size_t i)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:
ooura_fourier_cos_detail(const Real relative_error_goal,size_t levels)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>
integrate(F const & f,Real omega)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>
add_level(size_t i)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>
estimate_integral(F const & f,Real omega,size_t i)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
652