File indexing completed on 2025-01-18 09:39:53
0001
0002
0003
0004
0005
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
0025
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
0035 if (abs(x) > 0.125) {
0036 eta = 2*x - alpha*(1/expx - 1) + (expx - 1)/4;
0037 }
0038 else {
0039 eta = 2*x - alpha*expm1(-x) + expm1(x)/4;
0040 }
0041 return {eta, eta_prime};
0042 }
0043
0044
0045
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
0067
0068
0069
0070
0071
0072
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
0093 Real phi_prime = -(expm1_meta + x*exp_meta*eta_prime)/(expm1_meta*expm1_meta);
0094
0095
0096
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
0155
0156
0157
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
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
0243
0244
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
0261
0262
0263
0264
0265
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
0304
0305
0306 PreciseReal h = PreciseReal(1)/PreciseReal(1<<i);
0307
0308 std::vector<Real> bnode_row;
0309 std::vector<Real> bweight_row;
0310
0311
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
0342 } while(abs(w) > unit_roundoff*max_weight);
0343
0344
0345 bnode_row.shrink_to_fit();
0346 bweight_row.shrink_to_fit();
0347
0348
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
0361 break;
0362 }
0363 if (lnode_row.size() > 0) {
0364 if (lnode_row[lnode_row.size()-1] == node) {
0365
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
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
0387 std::lock_guard<std::mutex> lock(node_weight_mutex_);
0388 #endif
0389
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
0402
0403
0404 Real I0 = 0;
0405 auto const & b_nodes = big_nodes_[i];
0406 auto const & b_weights = bweights_[i];
0407
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
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
0426 std::vector<std::vector<Real>> big_nodes_;
0427
0428
0429 std::vector<std::vector<Real>> bweights_;
0430
0431
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
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
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
0592
0593
0594
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
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
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