]>
Commit | Line | Data |
---|---|---|
b32b8144 FG |
1 | // Copyright Nick Thompson, 2017 |
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 | ||
7 | #define BOOST_TEST_MODULE exp_sinh_quadrature_test | |
8 | ||
92f5a8d4 TL |
9 | #include <complex> |
10 | #include <boost/multiprecision/cpp_complex.hpp> | |
b32b8144 FG |
11 | #include <boost/math/concepts/real_concept.hpp> |
12 | #include <boost/test/included/unit_test.hpp> | |
92f5a8d4 | 13 | #include <boost/test/tools/floating_point_comparison.hpp> |
b32b8144 FG |
14 | #include <boost/math/quadrature/exp_sinh.hpp> |
15 | #include <boost/math/special_functions/sinc.hpp> | |
16 | #include <boost/math/special_functions/bessel.hpp> | |
17 | #include <boost/multiprecision/cpp_bin_float.hpp> | |
18 | #include <boost/multiprecision/cpp_dec_float.hpp> | |
19 | #include <boost/math/special_functions/next.hpp> | |
20 | #include <boost/math/special_functions/gamma.hpp> | |
21 | #include <boost/math/special_functions/sinc.hpp> | |
22 | #include <boost/type_traits/is_class.hpp> | |
23 | ||
92f5a8d4 TL |
24 | #ifdef BOOST_HAS_FLOAT128 |
25 | #include <boost/multiprecision/complex128.hpp> | |
26 | #endif | |
27 | ||
b32b8144 FG |
28 | using std::exp; |
29 | using std::cos; | |
30 | using std::tan; | |
31 | using std::log; | |
32 | using std::sqrt; | |
33 | using std::abs; | |
34 | using std::sinh; | |
35 | using std::cosh; | |
36 | using std::pow; | |
37 | using std::atan; | |
38 | using boost::multiprecision::cpp_bin_float_50; | |
39 | using boost::multiprecision::cpp_bin_float_100; | |
40 | using boost::multiprecision::cpp_bin_float_quad; | |
41 | using boost::math::constants::pi; | |
42 | using boost::math::constants::half_pi; | |
43 | using boost::math::constants::two_div_pi; | |
44 | using boost::math::constants::half; | |
45 | using boost::math::constants::third; | |
46 | using boost::math::constants::half; | |
47 | using boost::math::constants::third; | |
48 | using boost::math::constants::catalan; | |
49 | using boost::math::constants::ln_two; | |
50 | using boost::math::constants::root_two; | |
51 | using boost::math::constants::root_two_pi; | |
52 | using boost::math::constants::root_pi; | |
53 | using boost::math::quadrature::exp_sinh; | |
54 | ||
92f5a8d4 | 55 | #if !defined(TEST1) && !defined(TEST2) && !defined(TEST3) && !defined(TEST4) && !defined(TEST5) && !defined(TEST6) && !defined(TEST7) && !defined(TEST8) |
b32b8144 FG |
56 | # define TEST1 |
57 | # define TEST2 | |
58 | # define TEST3 | |
59 | # define TEST4 | |
60 | # define TEST5 | |
61 | # define TEST6 | |
62 | # define TEST7 | |
92f5a8d4 | 63 | # define TEST8 |
b32b8144 FG |
64 | #endif |
65 | ||
66 | #ifdef BOOST_MSVC | |
67 | #pragma warning (disable:4127) | |
68 | #endif | |
69 | ||
70 | // | |
71 | // Coefficient generation code: | |
72 | // | |
73 | template <class T> | |
74 | void print_levels(const T& v, const char* suffix) | |
75 | { | |
76 | std::cout << "{\n"; | |
77 | for (unsigned i = 0; i < v.size(); ++i) | |
78 | { | |
79 | std::cout << " { "; | |
80 | for (unsigned j = 0; j < v[i].size(); ++j) | |
81 | { | |
82 | std::cout << v[i][j] << suffix << ", "; | |
83 | } | |
84 | std::cout << "},\n"; | |
85 | } | |
86 | std::cout << " };\n"; | |
87 | } | |
88 | ||
89 | template <class T> | |
90 | void print_levels(const std::pair<T, T>& p, const char* suffix = "") | |
91 | { | |
92 | std::cout << " static const std::vector<std::vector<Real> > abscissa = "; | |
93 | print_levels(p.first, suffix); | |
94 | std::cout << " static const std::vector<std::vector<Real> > weights = "; | |
95 | print_levels(p.second, suffix); | |
96 | } | |
97 | ||
98 | template <class Real, class TargetType> | |
99 | std::pair<std::vector<std::vector<Real>>, std::vector<std::vector<Real>> > generate_constants(unsigned max_rows) | |
100 | { | |
101 | using boost::math::constants::half_pi; | |
102 | using boost::math::constants::two_div_pi; | |
103 | using boost::math::constants::pi; | |
92f5a8d4 TL |
104 | auto g = [](Real t)->Real { return exp(half_pi<Real>()*sinh(t)); }; |
105 | auto w = [](Real t)->Real { return cosh(t)*half_pi<Real>()*exp(half_pi<Real>()*sinh(t)); }; | |
b32b8144 FG |
106 | |
107 | std::vector<std::vector<Real>> abscissa, weights; | |
108 | ||
109 | std::vector<Real> temp; | |
110 | ||
111 | Real tmp = (Real(boost::math::tools::log_min_value<TargetType>()) + log(Real(boost::math::tools::epsilon<TargetType>())))*0.5f; | |
112 | Real t_min = asinh(two_div_pi<Real>()*tmp); | |
113 | // truncate t_min to an exact binary value: | |
114 | t_min = floor(t_min * 128) / 128; | |
115 | ||
116 | std::cout << "m_t_min = " << t_min << ";\n"; | |
117 | ||
118 | // t_max is chosen to make g'(t_max) ~ sqrt(max) (g' grows faster than g). | |
119 | // This will allow some flexibility on the users part; they can at least square a number function without overflow. | |
120 | // But there is no unique choice; the further out we can evaluate the function, the better we can do on slowly decaying integrands. | |
121 | const Real t_max = log(2 * two_div_pi<Real>()*log(2 * two_div_pi<Real>()*sqrt(Real(boost::math::tools::max_value<TargetType>())))); | |
122 | ||
123 | Real h = 1; | |
124 | for (Real t = t_min; t < t_max; t += h) | |
125 | { | |
126 | temp.push_back(g(t)); | |
127 | } | |
128 | abscissa.push_back(temp); | |
129 | temp.clear(); | |
130 | ||
131 | for (Real t = t_min; t < t_max; t += h) | |
132 | { | |
133 | temp.push_back(w(t * h)); | |
134 | } | |
135 | weights.push_back(temp); | |
136 | temp.clear(); | |
137 | ||
138 | for (unsigned row = 1; row < max_rows; ++row) | |
139 | { | |
140 | h /= 2; | |
141 | for (Real t = t_min + h; t < t_max; t += 2 * h) | |
142 | temp.push_back(g(t)); | |
143 | abscissa.push_back(temp); | |
144 | temp.clear(); | |
145 | } | |
146 | h = 1; | |
147 | for (unsigned row = 1; row < max_rows; ++row) | |
148 | { | |
149 | h /= 2; | |
150 | for (Real t = t_min + h; t < t_max; t += 2 * h) | |
151 | temp.push_back(w(t)); | |
152 | weights.push_back(temp); | |
153 | temp.clear(); | |
154 | } | |
155 | ||
156 | return std::make_pair(abscissa, weights); | |
157 | } | |
158 | ||
159 | ||
160 | template <class Real> | |
161 | const exp_sinh<Real>& get_integrator() | |
162 | { | |
163 | static const exp_sinh<Real> integrator(14); | |
164 | return integrator; | |
165 | } | |
166 | ||
167 | template <class Real> | |
168 | Real get_convergence_tolerance() | |
169 | { | |
170 | return boost::math::tools::root_epsilon<Real>(); | |
171 | } | |
172 | ||
173 | template<class Real> | |
174 | void test_right_limit_infinite() | |
175 | { | |
176 | std::cout << "Testing right limit infinite for tanh_sinh in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n"; | |
177 | Real tol = 10 * boost::math::tools::epsilon<Real>(); | |
178 | Real Q; | |
179 | Real Q_expected; | |
180 | Real error; | |
181 | Real L1; | |
182 | auto integrator = get_integrator<Real>(); | |
183 | ||
184 | // Example 12 | |
92f5a8d4 | 185 | const auto f2 = [](const Real& t)->Real { return exp(-t)/sqrt(t); }; |
b32b8144 FG |
186 | Q = integrator.integrate(f2, get_convergence_tolerance<Real>(), &error, &L1); |
187 | Q_expected = root_pi<Real>(); | |
188 | Real tol_mult = 1; | |
189 | // Multiprecision type have higher error rates, probably evaluation of f() is less accurate: | |
190 | if (std::numeric_limits<Real>::digits10 > std::numeric_limits<long double>::digits10) | |
191 | tol_mult = 12; | |
192 | else if (std::numeric_limits<Real>::digits10 > std::numeric_limits<double>::digits10) | |
193 | tol_mult = 5; | |
194 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol * tol_mult); | |
195 | // The integrand is strictly positive, so it coincides with the value of the integral: | |
196 | BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol * tol_mult); | |
197 | ||
198 | auto f3 = [](Real t)->Real { Real z = exp(-t); if (z == 0) { return z; } return z*cos(t); }; | |
199 | Q = integrator.integrate(f3, get_convergence_tolerance<Real>(), &error, &L1); | |
200 | Q_expected = half<Real>(); | |
201 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); | |
202 | Q = integrator.integrate(f3, 10, std::numeric_limits<Real>::has_infinity ? std::numeric_limits<Real>::infinity() : boost::math::tools::max_value<Real>(), get_convergence_tolerance<Real>(), &error, &L1); | |
203 | Q_expected = boost::lexical_cast<Real>("-6.6976341310426674140007086979326069121526743314567805278252392932e-6"); | |
204 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 10 * tol); | |
205 | // Integrating through zero risks precision loss: | |
206 | Q = integrator.integrate(f3, -10, std::numeric_limits<Real>::has_infinity ? std::numeric_limits<Real>::infinity() : boost::math::tools::max_value<Real>(), get_convergence_tolerance<Real>(), &error, &L1); | |
207 | Q_expected = boost::lexical_cast<Real>("-15232.3213626280525704332288302799653087046646639974940243044623285817777006"); | |
208 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, std::numeric_limits<Real>::digits10 > 30 ? 1000 * tol : tol); | |
209 | ||
210 | auto f4 = [](Real t)->Real { return 1/(1+t*t); }; | |
211 | Q = integrator.integrate(f4, 1, std::numeric_limits<Real>::has_infinity ? std::numeric_limits<Real>::infinity() : boost::math::tools::max_value<Real>(), get_convergence_tolerance<Real>(), &error, &L1); | |
212 | Q_expected = pi<Real>()/4; | |
213 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); | |
214 | BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol); | |
215 | Q = integrator.integrate(f4, 20, std::numeric_limits<Real>::has_infinity ? std::numeric_limits<Real>::infinity() : boost::math::tools::max_value<Real>(), get_convergence_tolerance<Real>(), &error, &L1); | |
216 | Q_expected = boost::lexical_cast<Real>("0.0499583957219427614100062870348448814912770804235071744108534548299835954767"); | |
217 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); | |
218 | BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol); | |
219 | Q = integrator.integrate(f4, 500, std::numeric_limits<Real>::has_infinity ? std::numeric_limits<Real>::infinity() : boost::math::tools::max_value<Real>(), get_convergence_tolerance<Real>(), &error, &L1); | |
220 | Q_expected = boost::lexical_cast<Real>("0.0019999973333397333150476759363217553199063513829126652556286269630"); | |
221 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); | |
222 | BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol); | |
223 | } | |
224 | ||
225 | template<class Real> | |
226 | void test_left_limit_infinite() | |
227 | { | |
228 | std::cout << "Testing left limit infinite for 1/(1+t^2) on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n"; | |
229 | Real tol = 10 * boost::math::tools::epsilon<Real>(); | |
230 | Real Q; | |
231 | Real Q_expected; | |
232 | Real error; | |
233 | Real L1; | |
234 | auto integrator = get_integrator<Real>(); | |
235 | ||
236 | // Example 11: | |
92f5a8d4 | 237 | auto f1 = [](const Real& t)->Real { return 1/(1+t*t);}; |
b32b8144 FG |
238 | Q = integrator.integrate(f1, std::numeric_limits<Real>::has_infinity ? -std::numeric_limits<Real>::infinity() : -boost::math::tools::max_value<Real>(), 0, get_convergence_tolerance<Real>(), &error, &L1); |
239 | Q_expected = half_pi<Real>(); | |
240 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); | |
241 | BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol); | |
242 | Q = integrator.integrate(f1, std::numeric_limits<Real>::has_infinity ? -std::numeric_limits<Real>::infinity() : -boost::math::tools::max_value<Real>(), -20, get_convergence_tolerance<Real>(), &error, &L1); | |
243 | Q_expected = boost::lexical_cast<Real>("0.0499583957219427614100062870348448814912770804235071744108534548299835954767"); | |
244 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); | |
245 | BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol); | |
246 | Q = integrator.integrate(f1, std::numeric_limits<Real>::has_infinity ? -std::numeric_limits<Real>::infinity() : -boost::math::tools::max_value<Real>(), -500, get_convergence_tolerance<Real>(), &error, &L1); | |
247 | Q_expected = boost::lexical_cast<Real>("0.0019999973333397333150476759363217553199063513829126652556286269630"); | |
248 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); | |
249 | BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol); | |
250 | } | |
251 | ||
252 | ||
253 | // Some examples of tough integrals from NR, section 4.5.4: | |
254 | template<class Real> | |
255 | void test_nr_examples() | |
256 | { | |
257 | using std::sin; | |
258 | using std::cos; | |
259 | using std::pow; | |
260 | using std::exp; | |
261 | using std::sqrt; | |
262 | std::cout << "Testing type " << boost::typeindex::type_id<Real>().pretty_name() << "\n"; | |
263 | Real tol = 10 * boost::math::tools::epsilon<Real>(); | |
264 | std::cout << std::setprecision(std::numeric_limits<Real>::digits10); | |
265 | Real Q; | |
266 | Real Q_expected; | |
267 | Real L1; | |
268 | Real error; | |
269 | auto integrator = get_integrator<Real>(); | |
270 | ||
271 | auto f0 = [] (Real)->Real { return (Real) 0; }; | |
272 | Q = integrator.integrate(f0, get_convergence_tolerance<Real>(), &error, &L1); | |
273 | Q_expected = 0; | |
274 | BOOST_CHECK_CLOSE_FRACTION(Q, 0.0f, tol); | |
275 | BOOST_CHECK_CLOSE_FRACTION(L1, 0.0f, tol); | |
276 | ||
92f5a8d4 | 277 | auto f = [](const Real& x)->Real { return 1/(1+x*x); }; |
b32b8144 FG |
278 | Q = integrator.integrate(f, get_convergence_tolerance<Real>(), &error, &L1); |
279 | Q_expected = half_pi<Real>(); | |
280 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); | |
281 | BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol); | |
282 | ||
283 | auto f1 = [](Real x)->Real { | |
284 | Real z1 = exp(-x); | |
285 | if (z1 == 0) | |
286 | { | |
287 | return (Real) 0; | |
288 | } | |
289 | Real z2 = pow(x, -3*half<Real>())*z1; | |
290 | if (z2 == 0) | |
291 | { | |
292 | return (Real) 0; | |
293 | } | |
294 | return sin(x*half<Real>())*z2; | |
295 | }; | |
296 | ||
297 | Q = integrator.integrate(f1, get_convergence_tolerance<Real>(), &error, &L1); | |
298 | Q_expected = sqrt(pi<Real>()*(sqrt((Real) 5) - 2)); | |
299 | ||
300 | // The integrand is oscillatory; the accuracy is low. | |
301 | Real tol_mul = 1; | |
302 | if (std::numeric_limits<Real>::digits10 > 40) | |
303 | tol_mul = 500000; | |
304 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol_mul * tol); | |
305 | ||
306 | auto f2 = [](Real x)->Real { return x > boost::math::tools::log_max_value<Real>() ? Real(0) : Real(pow(x, -(Real) 2/(Real) 7)*exp(-x*x)); }; | |
307 | Q = integrator.integrate(f2, get_convergence_tolerance<Real>(), &error, &L1); | |
308 | Q_expected = half<Real>()*boost::math::tgamma((Real) 5/ (Real) 14); | |
309 | tol_mul = 1; | |
310 | if (std::numeric_limits<Real>::is_specialized == false) | |
92f5a8d4 TL |
311 | tol_mul = 6; |
312 | else if (std::numeric_limits<Real>::digits10 > 40) | |
b32b8144 | 313 | tol_mul = 100; |
92f5a8d4 TL |
314 | else |
315 | tol_mul = 3; | |
b32b8144 FG |
316 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol_mul * tol); |
317 | BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol_mul * tol); | |
318 | ||
319 | auto f3 = [](Real x)->Real { return (Real) 1/ (sqrt(x)*(1+x)); }; | |
320 | Q = integrator.integrate(f3, get_convergence_tolerance<Real>(), &error, &L1); | |
321 | Q_expected = pi<Real>(); | |
322 | ||
323 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 10*boost::math::tools::epsilon<Real>()); | |
324 | BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, 10*boost::math::tools::epsilon<Real>()); | |
325 | ||
92f5a8d4 | 326 | auto f4 = [](const Real& t)->Real { return t > boost::math::tools::log_max_value<Real>() ? Real(0) : Real(exp(-t*t*half<Real>())); }; |
b32b8144 FG |
327 | Q = integrator.integrate(f4, get_convergence_tolerance<Real>(), &error, &L1); |
328 | Q_expected = root_two_pi<Real>()/2; | |
329 | tol_mul = 1; | |
330 | if (std::numeric_limits<Real>::digits10 > 40) | |
331 | tol_mul = 5000; | |
332 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol_mul * tol); | |
333 | BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol_mul * tol); | |
334 | ||
92f5a8d4 | 335 | auto f5 = [](const Real& t)->Real { return 1/cosh(t);}; |
b32b8144 FG |
336 | Q = integrator.integrate(f5, get_convergence_tolerance<Real>(), &error, &L1); |
337 | Q_expected = half_pi<Real>(); | |
338 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol * 12); // Fails at float precision without higher error rate | |
339 | BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol * 12); | |
340 | } | |
341 | ||
342 | // Definite integrals found in the CRC Handbook of Mathematical Formulas | |
343 | template<class Real> | |
344 | void test_crc() | |
345 | { | |
346 | using std::sin; | |
347 | using std::pow; | |
348 | using std::exp; | |
349 | using std::sqrt; | |
350 | using std::log; | |
351 | using std::cos; | |
352 | std::cout << "Testing integral from CRC handbook on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n"; | |
353 | Real tol = 10 * boost::math::tools::epsilon<Real>(); | |
354 | std::cout << std::setprecision(std::numeric_limits<Real>::digits10); | |
355 | Real Q; | |
356 | Real Q_expected; | |
357 | Real L1; | |
358 | Real error; | |
359 | auto integrator = get_integrator<Real>(); | |
360 | ||
92f5a8d4 | 361 | auto f0 = [](const Real& x)->Real { return x > boost::math::tools::log_max_value<Real>() ? Real(0) : Real(log(x)*exp(-x)); }; |
b32b8144 FG |
362 | Q = integrator.integrate(f0, get_convergence_tolerance<Real>(), &error, &L1); |
363 | Q_expected = -boost::math::constants::euler<Real>(); | |
364 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); | |
365 | ||
366 | // Test the integral representation of the gamma function: | |
367 | auto f1 = [](Real t)->Real { Real x = exp(-t); | |
368 | if(x == 0) | |
369 | { | |
370 | return (Real) 0; | |
371 | } | |
372 | return pow(t, (Real) 12 - 1)*x; | |
373 | }; | |
374 | ||
375 | Q = integrator.integrate(f1, get_convergence_tolerance<Real>(), &error, &L1); | |
376 | Q_expected = boost::math::tgamma(12.0f); | |
377 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); | |
378 | ||
379 | // Integral representation of the modified bessel function: | |
380 | // K_5(12) | |
381 | auto f2 = [](Real t)->Real { | |
382 | Real x = 12*cosh(t); | |
383 | if (x > boost::math::tools::log_max_value<Real>()) | |
384 | { | |
385 | return (Real) 0; | |
386 | } | |
387 | return exp(-x)*cosh(5*t); | |
388 | }; | |
389 | Q = integrator.integrate(f2, get_convergence_tolerance<Real>(), &error, &L1); | |
390 | Q_expected = boost::math::cyl_bessel_k<int, Real>(5, 12); | |
391 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); | |
392 | // Laplace transform of cos(at) | |
393 | Real a = 20; | |
394 | Real s = 1; | |
395 | auto f3 = [&](Real t)->Real { | |
396 | Real x = s * t; | |
397 | if (x > boost::math::tools::log_max_value<Real>()) | |
398 | { | |
399 | return (Real) 0; | |
400 | } | |
401 | return cos(a * t) * exp(-x); | |
402 | }; | |
403 | ||
b32b8144 FG |
404 | // Since the integrand is oscillatory, we increase the tolerance: |
405 | Real tol_mult = 10; | |
406 | // Multiprecision type have higher error rates, probably evaluation of f() is less accurate: | |
407 | if (!boost::is_class<Real>::value) | |
408 | { | |
92f5a8d4 TL |
409 | // For high oscillation frequency, the quadrature sum is ill-conditioned. |
410 | Q = integrator.integrate(f3, get_convergence_tolerance<Real>(), &error, &L1); | |
411 | Q_expected = s/(a*a+s*s); | |
b32b8144 FG |
412 | if (std::numeric_limits<Real>::digits10 > std::numeric_limits<double>::digits10) |
413 | tol_mult = 5000; // we should really investigate this more?? | |
414 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol_mult*tol); | |
415 | } | |
416 | ||
417 | // | |
418 | // This one doesn't pass for real_concept.. | |
419 | // | |
420 | if (std::numeric_limits<Real>::is_specialized) | |
421 | { | |
422 | // Laplace transform of J_0(t): | |
423 | auto f4 = [&](Real t)->Real { | |
424 | Real x = s * t; | |
425 | if (x > boost::math::tools::log_max_value<Real>()) | |
426 | { | |
427 | return (Real)0; | |
428 | } | |
429 | return boost::math::cyl_bessel_j(0, t) * exp(-x); | |
430 | }; | |
431 | ||
432 | Q = integrator.integrate(f4, get_convergence_tolerance<Real>(), &error, &L1); | |
433 | Q_expected = 1 / sqrt(1 + s*s); | |
434 | tol_mult = 3; | |
435 | // Multiprecision type have higher error rates, probably evaluation of f() is less accurate: | |
436 | if (std::numeric_limits<Real>::digits10 > std::numeric_limits<long double>::digits10) | |
437 | tol_mult = 750; | |
438 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol_mult * tol); | |
439 | } | |
92f5a8d4 | 440 | auto f6 = [](const Real& t)->Real { return t > boost::math::tools::log_max_value<Real>() ? Real(0) : Real(exp(-t*t)*log(t));}; |
b32b8144 FG |
441 | Q = integrator.integrate(f6, get_convergence_tolerance<Real>(), &error, &L1); |
442 | Q_expected = -boost::math::constants::root_pi<Real>()*(boost::math::constants::euler<Real>() + 2*ln_two<Real>())/4; | |
443 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); | |
444 | ||
445 | // CRC Section 5.5, integral 591 | |
446 | // The parameter p allows us to control the strength of the singularity. | |
447 | // Rapid convergence is not guaranteed for this function, as the branch cut makes it non-analytic on a disk. | |
448 | // This converges only when our test type has an extended exponent range as all the area of the integral | |
449 | // occurs so close to 0 (or 1) that we need abscissa values exceptionally small to find it. | |
450 | // "There's a lot of room at the bottom". | |
451 | // This version is transformed via argument substitution (exp(-x) for x) so that the integral is spread | |
452 | // over (0, INF). | |
453 | tol *= boost::math::tools::digits<Real>() > 100 ? 100000 : 75; | |
454 | for (Real pn = 99; pn > 0; pn -= 10) { | |
455 | Real p = pn / 100; | |
456 | auto f = [&](Real x)->Real | |
457 | { | |
458 | return x > 1000 * boost::math::tools::log_max_value<Real>() ? Real(0) : Real(exp(-x * (1 - p) + p * log(-boost::math::expm1(-x)))); | |
459 | }; | |
460 | Q = integrator.integrate(f, get_convergence_tolerance<Real>(), &error, &L1); | |
461 | Q_expected = 1 / boost::math::sinc_pi(p*pi<Real>()); | |
462 | /* | |
463 | std::cout << std::setprecision(std::numeric_limits<Real>::max_digits10) << p << std::endl; | |
464 | std::cout << std::setprecision(std::numeric_limits<Real>::max_digits10) << Q << std::endl; | |
465 | std::cout << std::setprecision(std::numeric_limits<Real>::max_digits10) << Q_expected << std::endl; | |
466 | std::cout << fabs((Q - Q_expected) / Q_expected) << std::endl; | |
467 | */ | |
468 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); | |
469 | } | |
470 | // and for p < 1: | |
471 | for (Real p = -0.99; p < 0; p += 0.1) { | |
472 | auto f = [&](Real x)->Real | |
473 | { | |
474 | return x > 1000 * boost::math::tools::log_max_value<Real>() ? Real(0) : Real(exp(-p * log(-boost::math::expm1(-x)) - (1 + p) * x)); | |
475 | }; | |
476 | Q = integrator.integrate(f, get_convergence_tolerance<Real>(), &error, &L1); | |
477 | Q_expected = 1 / boost::math::sinc_pi(p*pi<Real>()); | |
478 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); | |
479 | } | |
480 | } | |
481 | ||
92f5a8d4 TL |
482 | template<class Complex> |
483 | void test_complex_modified_bessel() | |
484 | { | |
485 | std::cout << "Testing complex modified Bessel function on type " << boost::typeindex::type_id<Complex>().pretty_name() << "\n"; | |
486 | typedef typename Complex::value_type Real; | |
487 | Real tol = 100 * boost::math::tools::epsilon<Real>(); | |
488 | Real error; | |
489 | Real L1; | |
490 | auto integrator = get_integrator<Real>(); | |
491 | ||
492 | // Integral Representation of Modified Complex Bessel function: | |
493 | // https://en.wikipedia.org/wiki/Bessel_function#Modified_Bessel_functions | |
494 | Complex z{2, 3}; | |
495 | const auto f = [&z](const Real& t)->Complex | |
496 | { | |
497 | using std::cosh; | |
498 | using std::exp; | |
499 | Real cosht = cosh(t); | |
500 | if (cosht > boost::math::tools::log_max_value<Real>()) | |
501 | { | |
502 | return Complex{0, 0}; | |
503 | } | |
504 | Complex arg = -z*cosht; | |
505 | Complex res = exp(arg); | |
506 | return res; | |
507 | }; | |
508 | ||
509 | Complex K0 = integrator.integrate(f, get_convergence_tolerance<Real>(), &error, &L1); | |
510 | ||
511 | // Mathematica code: N[BesselK[0, 2 + 3 I], 140] | |
512 | Real K0_x_expected = boost::lexical_cast<Real>("-0.08296852656762551490517953520589186885781541203818846830385526187936132191822538822296497597191327722262903004145527496422090506197776994"); | |
513 | Real K0_y_expected = boost::lexical_cast<Real>("0.027949603635183423629723306332336002340909030265538548521150904238352846705644065168365102147901993976999717171115546662967229050834575193041"); | |
514 | BOOST_CHECK_CLOSE_FRACTION(K0.real(), K0_x_expected, tol); | |
515 | BOOST_CHECK_CLOSE_FRACTION(K0.imag(), K0_y_expected, tol); | |
516 | } | |
517 | ||
b32b8144 FG |
518 | |
519 | BOOST_AUTO_TEST_CASE(exp_sinh_quadrature_test) | |
520 | { | |
521 | // | |
522 | // Uncomment to generate the coefficients: | |
523 | // | |
524 | ||
525 | /* | |
526 | std::cout << std::scientific << std::setprecision(8); | |
527 | print_levels(generate_constants<cpp_bin_float_100, float>(8), "f"); | |
528 | std::cout << std::setprecision(18); | |
529 | print_levels(generate_constants<cpp_bin_float_100, double>(8), ""); | |
530 | std::cout << std::setprecision(35); | |
531 | print_levels(generate_constants<cpp_bin_float_100, cpp_bin_float_quad>(8), "L"); | |
532 | */ | |
533 | ||
b32b8144 FG |
534 | #ifdef TEST1 |
535 | test_left_limit_infinite<float>(); | |
536 | test_right_limit_infinite<float>(); | |
537 | test_nr_examples<float>(); | |
538 | test_crc<float>(); | |
539 | #endif | |
540 | #ifdef TEST2 | |
541 | test_left_limit_infinite<double>(); | |
542 | test_right_limit_infinite<double>(); | |
543 | test_nr_examples<double>(); | |
544 | test_crc<double>(); | |
545 | #endif | |
11fdf7f2 | 546 | #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS |
b32b8144 FG |
547 | #ifdef TEST3 |
548 | test_left_limit_infinite<long double>(); | |
549 | test_right_limit_infinite<long double>(); | |
550 | test_nr_examples<long double>(); | |
551 | test_crc<long double>(); | |
552 | #endif | |
11fdf7f2 | 553 | #endif |
b32b8144 FG |
554 | #ifdef TEST4 |
555 | test_left_limit_infinite<cpp_bin_float_quad>(); | |
556 | test_right_limit_infinite<cpp_bin_float_quad>(); | |
557 | test_nr_examples<cpp_bin_float_quad>(); | |
558 | test_crc<cpp_bin_float_quad>(); | |
559 | #endif | |
560 | ||
11fdf7f2 | 561 | #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS |
b32b8144 FG |
562 | #ifdef TEST5 |
563 | test_left_limit_infinite<boost::math::concepts::real_concept>(); | |
564 | test_right_limit_infinite<boost::math::concepts::real_concept>(); | |
565 | test_nr_examples<boost::math::concepts::real_concept>(); | |
566 | test_crc<boost::math::concepts::real_concept>(); | |
567 | #endif | |
11fdf7f2 | 568 | #endif |
b32b8144 FG |
569 | #ifdef TEST6 |
570 | test_left_limit_infinite<boost::multiprecision::cpp_bin_float_50>(); | |
571 | test_right_limit_infinite<boost::multiprecision::cpp_bin_float_50>(); | |
572 | test_nr_examples<boost::multiprecision::cpp_bin_float_50>(); | |
573 | test_crc<boost::multiprecision::cpp_bin_float_50>(); | |
574 | #endif | |
575 | #ifdef TEST7 | |
576 | test_left_limit_infinite<boost::multiprecision::cpp_dec_float_50>(); | |
577 | test_right_limit_infinite<boost::multiprecision::cpp_dec_float_50>(); | |
578 | test_nr_examples<boost::multiprecision::cpp_dec_float_50>(); | |
92f5a8d4 TL |
579 | // |
580 | // This one causes stack overflows on the CI machine, but not locally, | |
581 | // assume it's due to resticted resources on the server, and <shrug> for now... | |
582 | // | |
583 | #if ! BOOST_WORKAROUND(BOOST_MSVC, == 1900) | |
b32b8144 FG |
584 | test_crc<boost::multiprecision::cpp_dec_float_50>(); |
585 | #endif | |
92f5a8d4 TL |
586 | #endif |
587 | #ifdef TEST8 | |
588 | test_complex_modified_bessel<std::complex<float>>(); | |
589 | test_complex_modified_bessel<std::complex<double>>(); | |
590 | test_complex_modified_bessel<std::complex<long double>>(); | |
591 | #ifdef BOOST_HAS_FLOAT128 | |
592 | test_complex_modified_bessel<boost::multiprecision::complex128>(); | |
593 | #endif | |
594 | test_complex_modified_bessel<boost::multiprecision::cpp_complex_quad>(); | |
595 | #endif | |
b32b8144 | 596 | } |