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