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 #define BOOST_TEST_MODULE test_ooura_fourier_transform
10 #include <boost/type_index.hpp>
11 #include <boost/test/included/unit_test.hpp>
12 #include <boost/test/tools/floating_point_comparison.hpp>
13 #include <boost/math/quadrature/ooura_fourier_integrals.hpp>
14 #include <boost/multiprecision/cpp_bin_float.hpp>
16 using boost::math::quadrature::ooura_fourier_sin
;
17 using boost::math::quadrature::ooura_fourier_cos
;
18 using boost::math::constants::pi
;
21 float float_tol
= 10*std::numeric_limits
<float>::epsilon();
22 ooura_fourier_sin
<float> float_sin_integrator(float_tol
);
24 double double_tol
= 10*std::numeric_limits
<double>::epsilon();
25 ooura_fourier_sin
<double> double_sin_integrator(double_tol
);
27 long double long_double_tol
= 10*std::numeric_limits
<long double>::epsilon();
28 ooura_fourier_sin
<long double> long_double_sin_integrator(long_double_tol
);
31 auto get_sin_integrator() {
32 if constexpr (std::is_same_v
<Real
, float>) {
33 return float_sin_integrator
;
35 if constexpr (std::is_same_v
<Real
, double>) {
36 return double_sin_integrator
;
38 if constexpr (std::is_same_v
<Real
, long double>) {
39 return long_double_sin_integrator
;
43 ooura_fourier_cos
<float> float_cos_integrator(float_tol
);
44 ooura_fourier_cos
<double> double_cos_integrator(double_tol
);
45 ooura_fourier_cos
<long double> long_double_cos_integrator(long_double_tol
);
48 auto get_cos_integrator() {
49 if constexpr (std::is_same_v
<Real
, float>) {
50 return float_cos_integrator
;
52 if constexpr (std::is_same_v
<Real
, double>) {
53 return double_cos_integrator
;
55 if constexpr (std::is_same_v
<Real
, long double>) {
56 return long_double_cos_integrator
;
64 using boost::math::quadrature::detail::ooura_eta
;
65 std::cout
<< "Testing eta function on type " << boost::typeindex::type_id
<Real
>().pretty_name() << "\n";
69 auto [eta
, eta_prime
] = ooura_eta(x
, alpha
);
70 BOOST_CHECK_SMALL(eta
, (std::numeric_limits
<Real
>::min
)());
71 BOOST_CHECK_CLOSE_FRACTION(eta_prime
, 2 + alpha
+ Real(1)/Real(4), 10*std::numeric_limits
<Real
>::epsilon());
76 for (Real z
= 0.125; z
< 500; z
+= 0.125) {
78 auto [eta
, eta_prime
] = ooura_eta(x
, alpha
);
79 BOOST_CHECK_CLOSE_FRACTION(eta
, 2*x
+ alpha
*(1-1/z
) + (z
-1)/4, 10*std::numeric_limits
<Real
>::epsilon());
80 BOOST_CHECK_CLOSE_FRACTION(eta_prime
, 2 + alpha
/z
+ z
/4, 10*std::numeric_limits
<Real
>::epsilon());
86 void test_ooura_sin_nodes_and_weights()
88 using boost::math::quadrature::detail::ooura_sin_node_and_weight
;
89 using boost::math::quadrature::detail::ooura_eta
;
90 std::cout
<< "Testing nodes and weights on type " << boost::typeindex::type_id
<Real
>().pretty_name() << "\n";
95 auto [node
, weight
] = ooura_sin_node_and_weight(n
, h
, alpha
);
96 Real expected_node
= pi
<Real
>()/(1-exp(-ooura_eta(n
*h
, alpha
).first
));
97 BOOST_CHECK_CLOSE_FRACTION(node
, expected_node
,10*std::numeric_limits
<Real
>::epsilon());
102 void test_ooura_alpha() {
103 std::cout
<< "Testing Ooura alpha on type " << boost::typeindex::type_id
<Real
>().pretty_name() << "\n";
106 using boost::math::quadrature::detail::calculate_ooura_alpha
;
107 Real alpha
= calculate_ooura_alpha(Real(1));
108 Real expected
= 1/sqrt(16 + 4*log1p(pi
<Real
>()));
109 BOOST_CHECK_CLOSE_FRACTION(alpha
, expected
, 10*std::numeric_limits
<Real
>::epsilon());
112 void test_node_weight_precision_agreement()
115 using boost::math::quadrature::detail::ooura_sin_node_and_weight
;
116 using boost::math::quadrature::detail::ooura_eta
;
117 using boost::multiprecision::cpp_bin_float_quad
;
118 std::cout
<< "Testing agreement in two different precisions of nodes and weights\n";
119 cpp_bin_float_quad alpha_quad
= 1;
121 cpp_bin_float_quad h_quad
= 1/cpp_bin_float_quad(int_max
);
122 double alpha_dbl
= 1;
123 double h_dbl
= static_cast<double>(h_quad
);
124 std::cout
<< std::fixed
;
125 for (long n
= -1; n
> -6*int_max
; --n
) {
126 auto [node_dbl
, weight_dbl
] = ooura_sin_node_and_weight(n
, h_dbl
, alpha_dbl
);
127 auto p
= ooura_sin_node_and_weight(n
, h_quad
, alpha_quad
);
128 double node_quad
= static_cast<double>(p
.first
);
129 double weight_quad
= static_cast<double>(p
.second
);
130 auto node_dist
= abs(boost::math::float_distance(node_quad
, node_dbl
));
131 if ( (weight_quad
< 0 && weight_dbl
> 0) || (weight_dbl
< 0 && weight_quad
> 0) ){
132 std::cout
<< "Weights at different precisions have different signs!\n";
134 auto weight_dist
= abs(boost::math::float_distance(weight_quad
, weight_dbl
));
135 if (weight_dist
> 100) {
136 std::cout
<< std::fixed
;
137 std::cout
<<"n =" << n
<< ", x = " << n
*h_dbl
<< ", node distance = " << node_dist
<< ", weight distance = " << weight_dist
<< "\n";
138 std::cout
<< std::scientific
;
139 std::cout
<< "computed weight = " << weight_dbl
<< ", actual weight = " << weight_quad
<< "\n";
149 std::cout
<< "Testing sinc integral on type " << boost::typeindex::type_id
<Real
>().pretty_name() << "\n";
150 using std::numeric_limits
;
151 Real tol
= 50*numeric_limits
<Real
>::epsilon();
152 auto integrator
= get_sin_integrator
<Real
>();
153 auto f
= [](Real x
)->Real
{ return 1/x
; };
157 auto [Is
, err
] = integrator
.integrate(f
, omega
);
158 BOOST_CHECK_CLOSE_FRACTION(Is
, pi
<Real
>()/2, tol
);
160 auto [Isn
, errn
] = integrator
.integrate(f
, -omega
);
161 BOOST_CHECK_CLOSE_FRACTION(Isn
, -pi
<Real
>()/2, tol
);
170 std::cout
<< "Testing exponential integral on type " << boost::typeindex::type_id
<Real
>().pretty_name() << "\n";
172 using std::numeric_limits
;
173 Real tol
= 50*numeric_limits
<Real
>::epsilon();
174 auto integrator
= get_sin_integrator
<Real
>();
175 auto f
= [](Real x
)->Real
{return exp(-x
);};
179 auto [Is
, err
] = integrator
.integrate(f
, omega
);
180 Real exact
= omega
/(1+omega
*omega
);
181 BOOST_CHECK_CLOSE_FRACTION(Is
, exact
, tol
);
190 std::cout
<< "Testing integral of sin(kx)/sqrt(x) on type " << boost::typeindex::type_id
<Real
>().pretty_name() << "\n";
192 using std::numeric_limits
;
193 Real tol
= 10*numeric_limits
<Real
>::epsilon();
194 auto integrator
= get_sin_integrator
<Real
>();
195 auto f
= [](Real x
)->Real
{ return 1/sqrt(x
);};
198 auto [Is
, err
] = integrator
.integrate(f
, omega
);
199 Real exact
= sqrt(pi
<Real
>()/(2*omega
));
200 BOOST_CHECK_CLOSE_FRACTION(Is
, exact
, 10*tol
);
205 // See: https://scicomp.stackexchange.com/questions/32790/numerical-evaluation-of-highly-oscillatory-integral/32799#32799
207 Real
asymptotic(Real lambda
) {
210 using boost::math::constants::pi
;
211 Real I1
= cos(lambda
- pi
<Real
>()/4)*sqrt(2*pi
<Real
>()/lambda
);
212 Real I2
= sin(lambda
- pi
<Real
>()/4)*sqrt(2*pi
<Real
>()/(lambda
*lambda
*lambda
))/8;
217 void test_double_osc()
219 std::cout
<< "Testing double oscillation on type " << boost::typeindex::type_id
<Real
>().pretty_name() << "\n";
221 using std::numeric_limits
;
222 auto integrator
= get_sin_integrator
<Real
>();
224 auto f
= [&lambda
](Real x
)->Real
{ return cos(lambda
*cos(x
))/x
; };
226 auto [Is
, err
] = integrator
.integrate(f
, omega
);
227 Real exact
= asymptotic(lambda
);
228 BOOST_CHECK_CLOSE_FRACTION(2*Is
, exact
, 0.05);
232 void test_zero_integrand()
234 // Make sure relative error tolerance doesn't break on zero integrand:
235 std::cout
<< "Testing zero integrand on type " << boost::typeindex::type_id
<Real
>().pretty_name() << "\n";
237 using std::numeric_limits
;
238 auto integrator
= get_sin_integrator
<Real
>();
239 auto f
= [](Real
/* x */)->Real
{ return Real(0); };
241 auto [Is
, err
] = integrator
.integrate(f
, omega
);
243 BOOST_CHECK_EQUAL(Is
, exact
);
247 // This works, but doesn't recover the precision you want in a unit test:
248 // template<class Real>
251 // std::cout << "Testing integral of log(x)sin(x) on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
254 // using std::numeric_limits;
255 // using boost::math::constants::euler;
256 // Real tol = 1000*numeric_limits<Real>::epsilon();
257 // auto f = [](Real x)->Real { return exp(-100*numeric_limits<Real>::epsilon()*x)*log(x);};
259 // Real Is = ooura_fourier_sin<decltype(f), Real>(f, omega, sqrt(numeric_limits<Real>::epsilon())/100);
260 // BOOST_CHECK_CLOSE_FRACTION(Is, -euler<Real>(), tol);
265 void test_cos_integral1()
267 std::cout
<< "Testing integral of cos(x)/(x*x+1) on type " << boost::typeindex::type_id
<Real
>().pretty_name() << "\n";
269 using boost::math::constants::half_pi
;
270 using boost::math::constants::e
;
271 using std::numeric_limits
;
272 Real tol
= 10*numeric_limits
<Real
>::epsilon();
274 auto integrator
= get_cos_integrator
<Real
>();
275 auto f
= [](Real x
)->Real
{ return 1/(x
*x
+1);};
277 auto [Is
, err
] = integrator
.integrate(f
, omega
);
278 Real exact
= half_pi
<Real
>()/e
<Real
>();
279 BOOST_CHECK_CLOSE_FRACTION(Is
, exact
, tol
);
283 void test_cos_integral2()
285 std::cout
<< "Testing integral of exp(-a*x) on type " << boost::typeindex::type_id
<Real
>().pretty_name() << "\n";
287 using boost::math::constants::half_pi
;
288 using boost::math::constants::e
;
289 using std::numeric_limits
;
290 Real tol
= 10*numeric_limits
<Real
>::epsilon();
292 auto integrator
= get_cos_integrator
<Real
>();
293 for (Real a
= 1; a
< 5; ++a
) {
294 auto f
= [&a
](Real x
)->Real
{ return exp(-a
*x
);};
295 for(Real omega
= 1; omega
< 5; ++omega
) {
296 auto [Is
, err
] = integrator
.integrate(f
, omega
);
297 Real exact
= a
/(a
*a
+omega
*omega
);
298 BOOST_CHECK_CLOSE_FRACTION(Is
, exact
, 10*tol
);
306 std::cout
<< "Testing nodes and weights on type " << boost::typeindex::type_id
<Real
>().pretty_name() << "\n";
307 auto sin_integrator
= get_sin_integrator
<Real
>();
309 auto const & big_nodes
= sin_integrator
.big_nodes();
310 for (auto & node_row
: big_nodes
) {
311 Real t0
= node_row
[0];
312 for (size_t i
= 1; i
< node_row
.size(); ++i
) {
313 Real t1
= node_row
[i
];
314 BOOST_CHECK(t1
> t0
);
319 auto const & little_nodes
= sin_integrator
.little_nodes();
320 for (auto & node_row
: little_nodes
) {
321 Real t0
= node_row
[0];
322 for (size_t i
= 1; i
< node_row
.size(); ++i
) {
323 Real t1
= node_row
[i
];
324 BOOST_CHECK(t1
< t0
);
331 BOOST_AUTO_TEST_CASE(ooura_fourier_transform_test
)
333 test_cos_integral1
<float>();
334 test_cos_integral1
<double>();
335 test_cos_integral1
<long double>();
337 test_cos_integral2
<float>();
338 test_cos_integral2
<double>();
339 test_cos_integral2
<long double>();
341 //test_node_weight_precision_agreement();
342 test_zero_integrand
<float>();
343 test_zero_integrand
<double>();
345 test_ooura_eta
<float>();
346 test_ooura_eta
<double>();
347 test_ooura_eta
<long double>();
349 test_ooura_sin_nodes_and_weights
<float>();
350 test_ooura_sin_nodes_and_weights
<double>();
351 test_ooura_sin_nodes_and_weights
<long double>();
353 test_ooura_alpha
<float>();
354 test_ooura_alpha
<double>();
355 test_ooura_alpha
<long double>();
359 test_sinc
<long double>();
363 test_exp
<long double>();
368 test_double_osc
<float>();
369 test_double_osc
<double>();
371 //test_double_osc<long double>();
373 // This test should be last:
375 test_nodes
<double>();
376 test_nodes
<long double>();