1 // Copyright Paul A. Bristow, 2019
2 // Copyright Nick Thompson, 2019
3
4 // Use, modification and distribution are subject to the
5 // Boost Software License, Version 1.0.
6 // (See accompanying file LICENSE_1_0.txt
7 // or copy at http://www.boost.org/LICENSE_1_0.txt)
8
9 #if !defined(__cpp_structured_bindings) || (__cpp_structured_bindings < 201606L)
10 # error "This example requires a C++17 compiler that supports 'structured bindings'. Try /std:c++17 or -std=c++17 or later."
11 #endif
12
13 //#define BOOST_MATH_INSTRUMENT_OOURA // or -DBOOST_MATH_INSTRUMENT_OOURA etc for diagnostic output.
14
15 #include <boost/math/quadrature/ooura_fourier_integrals.hpp>
16 #include <boost/multiprecision/cpp_bin_float.hpp> // for cpp_bin_float_quad, cpp_bin_float_50...
17 #include <boost/math/constants/constants.hpp> // For pi (including for multiprecision types, if used.)
18
19 #include <cmath>
20 #include <iostream>
21 #include <limits>
22 #include <iostream>
23 #include <exception>
24
main()25 int main()
26 {
27 try
28 {
29 typedef boost::multiprecision::cpp_bin_float_quad Real;
30
31 std::cout.precision(std::numeric_limits<Real>::max_digits10); // Show all potentially significant digits.
32
33 using boost::math::quadrature::ooura_fourier_cos;
34 using boost::math::constants::half_pi;
35 using boost::math::constants::e;
36
37 //[ooura_fourier_integrals_multiprecision_example_1
38
39 // Use the default parameters for tolerance root_epsilon and eight levels for a type of 8 bytes.
40 //auto integrator = ooura_fourier_cos<Real>();
41 // Decide on a (tight) tolerance.
42 const Real tol = 2 * std::numeric_limits<Real>::epsilon();
43 auto integrator = ooura_fourier_cos<Real>(tol, 8); // Loops or gets worse for more than 8.
44
45 auto f = [](Real x)
46 { // More complex example function.
47 return 1 / (x * x + 1);
48 };
49
50 double omega = 1;
51 auto [result, relative_error] = integrator.integrate(f, omega);
52
53 //] [/ooura_fourier_integrals_multiprecision_example_1]
54
55 //[ooura_fourier_integrals_multiprecision_example_2
56 std::cout << "Integral = " << result << ", relative error estimate " << relative_error << std::endl;
57
58 const Real expected = half_pi<Real>() / e<Real>(); // Expect integral = 1/(2e)
59 std::cout << "pi/(2e) = " << expected << ", difference " << result - expected << std::endl;
60 //] [/ooura_fourier_integrals_multiprecision_example_2]
61 }
62 catch (std::exception const & ex)
63 {
64 // Lacking try&catch blocks, the program will abort after any throw, whereas the
65 // message below from the thrown exception will give some helpful clues as to the cause of the problem.
66 std::cout << "\n""Message from thrown exception was:\n " << ex.what() << std::endl;
67 }
68 } // int main()
69
70 /*
71
72 //[ooura_fourier_integrals_example_multiprecision_output_1
73 ``
74 Integral = 0.5778636748954608589550465916563501587, relative error estimate 4.609814684522163895264277312610830278e-17
75 pi/(2e) = 0.5778636748954608659545328919193707407, difference -6.999486300263020581921171645255733758e-18
76 ``
77 //] [/ooura_fourier_integrals_example_multiprecision_output_1]
78
79
80 //[ooura_fourier_integrals_example_multiprecision_diagnostic_output_1
81 ``
82 ooura_fourier_cos with relative error goal 3.851859888774471706111955885169854637e-34 & 15 levels.
83 epsilon for type = 1.925929944387235853055977942584927319e-34
84 h = 1.000000000000000000000000000000000, I_h = 0.588268622591776615359568690603776 = 0.5882686225917766153595686906037760, absolute error estimate = nan
85 h = 0.500000000000000000000000000000000, I_h = 0.577871642184837461311756940493259 = 0.5778716421848374613117569404932595, absolute error estimate = 1.039698040693915404781175011051656e-02
86 h = 0.250000000000000000000000000000000, I_h = 0.577863671186882539559996800783122 = 0.5778636711868825395599968007831220, absolute error estimate = 7.970997954921751760139710137450075e-06
87 h = 0.125000000000000000000000000000000, I_h = 0.577863674895460885593491133506723 = 0.5778636748954608855934911335067232, absolute error estimate = 3.708578346033494332723601147051768e-09
88 h = 0.062500000000000000000000000000000, I_h = 0.577863674895460858955046591656350 = 0.5778636748954608589550465916563502, absolute error estimate = 2.663844454185037302771663314961535e-17
89 h = 0.031250000000000000000000000000000, I_h = 0.577863674895460858955046591656348 = 0.5778636748954608589550465916563484, absolute error estimate = 1.733336949948512267750380148326435e-33
90 h = 0.015625000000000000000000000000000, I_h = 0.577863674895460858955046591656348 = 0.5778636748954608589550465916563479, absolute error estimate = 4.814824860968089632639944856462318e-34
91 h = 0.007812500000000000000000000000000, I_h = 0.577863674895460858955046591656347 = 0.5778636748954608589550465916563473, absolute error estimate = 6.740754805355325485695922799047246e-34
92 h = 0.003906250000000000000000000000000, I_h = 0.577863674895460858955046591656347 = 0.5778636748954608589550465916563475, absolute error estimate = 1.925929944387235853055977942584927e-34
93 Integral = 5.778636748954608589550465916563475e-01, relative error estimate 3.332844800697411177051445985473052e-34
94 pi/(2e) = 5.778636748954608589550465916563481e-01, difference -6.740754805355325485695922799047246e-34
95 ``
96 //] [/ooura_fourier_integrals_example_multiprecision_diagnostic_output_1]
97
98 */
99