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)
7 #define BOOST_TEST_MODULE tanh_sinh_quadrature_test
9 #include <boost/config.hpp>
10 #include <boost/detail/workaround.hpp>
12 #if !defined(BOOST_NO_CXX11_DECLTYPE) && !defined(BOOST_NO_CXX11_TRAILING_RESULT_TYPES) && !defined(BOOST_NO_SFINAE_EXPR)
14 #include <boost/math/concepts/real_concept.hpp>
15 #include <boost/test/included/unit_test.hpp>
16 #include <boost/test/floating_point_comparison.hpp>
17 #include <boost/math/quadrature/tanh_sinh.hpp>
18 #include <boost/math/special_functions/sinc.hpp>
19 #include <boost/multiprecision/cpp_bin_float.hpp>
20 #include <boost/multiprecision/cpp_dec_float.hpp>
21 #include <boost/math/special_functions/next.hpp>
22 #include <boost/math/special_functions/gamma.hpp>
23 #include <boost/math/special_functions/beta.hpp>
24 #include <boost/math/special_functions/ellint_rc.hpp>
25 #include <boost/math/special_functions/ellint_rj.hpp>
27 #ifdef BOOST_HAS_FLOAT128
28 #include <boost/multiprecision/float128.hpp>
32 #pragma warning(disable:4127) // Conditional expression is constant
35 #if !defined(TEST1) && !defined(TEST2) && !defined(TEST3) && !defined(TEST4) && !defined(TEST5) && !defined(TEST6) && !defined(TEST7) && !defined(TEST8)\
36 && !defined(TEST1A) && !defined(TEST2A) && !defined(TEST3A) && !defined(TEST6A)
69 using boost::multiprecision::cpp_bin_float_50
;
70 using boost::multiprecision::cpp_bin_float_100
;
71 using boost::multiprecision::cpp_dec_float_50
;
72 using boost::multiprecision::cpp_dec_float_100
;
73 using boost::multiprecision::cpp_bin_float_quad
;
74 using boost::math::sinc_pi
;
75 using boost::math::quadrature::tanh_sinh
;
76 using boost::math::quadrature::detail::tanh_sinh_detail
;
77 using boost::math::constants::pi
;
78 using boost::math::constants::half_pi
;
79 using boost::math::constants::two_div_pi
;
80 using boost::math::constants::two_pi
;
81 using boost::math::constants::half
;
82 using boost::math::constants::third
;
83 using boost::math::constants::half
;
84 using boost::math::constants::third
;
85 using boost::math::constants::catalan
;
86 using boost::math::constants::ln_two
;
87 using boost::math::constants::root_two
;
88 using boost::math::constants::root_two_pi
;
89 using boost::math::constants::root_pi
;
92 void print_levels(const T
& v
, const char* suffix
)
95 for (unsigned i
= 0; i
< v
.size(); ++i
)
98 for (unsigned j
= 0; j
< v
[i
].size(); ++j
)
100 std::cout
<< v
[i
][j
] << suffix
<< ", ";
104 std::cout
<< " };\n";
108 void print_complement_indexes(const T
& v
)
111 for (unsigned i
= 0; i
< v
.size(); ++i
)
114 while (v
[i
][index
] >= 0)
116 std::cout
<< index
<< ", ";
122 void print_levels(const std::pair
<T
, T
>& p
, const char* suffix
= "")
124 std::cout
<< " static const std::vector<std::vector<Real> > abscissa = ";
125 print_levels(p
.first
, suffix
);
126 std::cout
<< " static const std::vector<std::vector<Real> > weights = ";
127 print_levels(p
.second
, suffix
);
128 std::cout
<< " static const std::vector<unsigned > indexes = ";
129 print_complement_indexes(p
.first
);
132 template <class Real
>
133 std::pair
<std::vector
<std::vector
<Real
>>, std::vector
<std::vector
<Real
>> > generate_constants(unsigned max_index
, unsigned max_rows
)
135 using boost::math::constants::half_pi
;
136 using boost::math::constants::two_div_pi
;
137 using boost::math::constants::pi
;
138 auto g
= [](Real t
) { return tanh(half_pi
<Real
>()*sinh(t
)); };
139 auto w
= [](Real t
) { Real cs
= cosh(half_pi
<Real
>() * sinh(t
)); return half_pi
<Real
>() * cosh(t
) / (cs
* cs
); };
140 auto gc
= [](Real t
) { Real u2
= half_pi
<Real
>() * sinh(t
); return 1 / (exp(u2
) *cosh(u2
)); };
141 auto g_inv
= [](float x
)->float { return asinh(two_div_pi
<float>()*atanh(x
)); };
142 auto gc_inv
= [](float x
)
144 float l
= log(sqrt((2 - x
) / x
));
145 return log((sqrt(4 * l
* l
+ pi
<float>() * pi
<float>()) + 2 * l
) / pi
<float>());
148 std::vector
<std::vector
<Real
>> abscissa
, weights
;
150 std::vector
<Real
> temp
;
152 float t_crossover
= gc_inv(0.5f
);
155 for (unsigned i
= 0; i
< max_index
; ++i
)
157 temp
.push_back((i
< t_crossover
) ? g(i
* h
) : -gc(i
* h
));
159 abscissa
.push_back(temp
);
162 for (unsigned i
= 0; i
< max_index
; ++i
)
164 temp
.push_back(w(i
* h
));
166 weights
.push_back(temp
);
169 for (unsigned row
= 1; row
< max_rows
; ++row
)
172 for (Real t
= h
; t
< max_index
- 1; t
+= 2 * h
)
173 temp
.push_back((t
< t_crossover
) ? g(t
) : -gc(t
));
174 abscissa
.push_back(temp
);
178 for (unsigned row
= 1; row
< max_rows
; ++row
)
181 for (Real t
= h
; t
< max_index
- 1; t
+= 2 * h
)
182 temp
.push_back(w(t
));
183 weights
.push_back(temp
);
187 return std::make_pair(abscissa
, weights
);
190 template <class Real
>
191 const tanh_sinh
<Real
>& get_integrator()
193 static const tanh_sinh
<Real
> integrator(15);
197 template <class Real
>
198 Real
get_convergence_tolerance()
200 return boost::math::tools::root_epsilon
<Real
>();
207 std::cout
<< "Testing linear functions are integrated properly by tanh_sinh on type " << boost::typeindex::type_id
<Real
>().pretty_name() << "\n";
208 Real tol
= 10*boost::math::tools::epsilon
<Real
>();
209 auto integrator
= get_integrator
<Real
>();
210 auto f
= [](const Real
& x
)
216 Real Q
= integrator
.integrate(f
, (Real
) 0, (Real
) 1, get_convergence_tolerance
<Real
>(), &error
, &L1
);
217 BOOST_CHECK_CLOSE_FRACTION(Q
, 9.5, tol
);
218 BOOST_CHECK_CLOSE_FRACTION(L1
, 9.5, tol
);
223 void test_quadratic()
225 std::cout
<< "Testing quadratic functions are integrated properly by tanh_sinh on type " << boost::typeindex::type_id
<Real
>().pretty_name() << "\n";
226 Real tol
= 10*boost::math::tools::epsilon
<Real
>();
227 auto integrator
= get_integrator
<Real
>();
228 auto f
= [](const Real
& x
) { return 5*x
*x
+ 7*x
+ 12; };
231 Real Q
= integrator
.integrate(f
, (Real
) 0, (Real
) 1, get_convergence_tolerance
<Real
>(), &error
, &L1
);
232 BOOST_CHECK_CLOSE_FRACTION(Q
, (Real
) 17 + half
<Real
>()*third
<Real
>(), tol
);
233 BOOST_CHECK_CLOSE_FRACTION(L1
, (Real
) 17 + half
<Real
>()*third
<Real
>(), tol
);
240 std::cout
<< "Testing integration of endpoint singularities on type " << boost::typeindex::type_id
<Real
>().pretty_name() << "\n";
241 Real tol
= 10*boost::math::tools::epsilon
<Real
>();
244 auto integrator
= get_integrator
<Real
>();
245 auto f
= [](const Real
& x
) { return log(x
)*log(1-x
); };
246 Real Q
= integrator
.integrate(f
, (Real
) 0, (Real
) 1, get_convergence_tolerance
<Real
>(), &error
, &L1
);
247 Real Q_expected
= 2 - pi
<Real
>()*pi
<Real
>()*half
<Real
>()*third
<Real
>();
249 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
250 BOOST_CHECK_CLOSE_FRACTION(L1
, Q_expected
, tol
);
254 // Examples taken from
255 //http://crd-legacy.lbl.gov/~dhbailey/dhbpapers/quadrature.pdf
259 std::cout
<< "Testing integration of C(a) on type " << boost::typeindex::type_id
<Real
>().pretty_name() << "\n";
260 Real tol
= 10 * boost::math::tools::epsilon
<Real
>();
264 auto integrator
= get_integrator
<Real
>();
265 auto f1
= [](const Real
& x
) { return atan(x
)/(x
*(x
*x
+ 1)) ; };
266 Real Q
= integrator
.integrate(f1
, (Real
) 0, (Real
) 1, get_convergence_tolerance
<Real
>(), &error
, &L1
);
267 Real Q_expected
= pi
<Real
>()*ln_two
<Real
>()/8 + catalan
<Real
>()*half
<Real
>();
268 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
269 BOOST_CHECK_CLOSE_FRACTION(L1
, Q_expected
, tol
);
271 auto f2
= [](Real x
)->Real
{ Real t0
= x
*x
+ 1; Real t1
= sqrt(t0
); return atan(t1
)/(t0
*t1
); };
272 Q
= integrator
.integrate(f2
, (Real
) 0 , (Real
) 1, get_convergence_tolerance
<Real
>(), &error
, &L1
);
273 Q_expected
= pi
<Real
>()/4 - pi
<Real
>()/root_two
<Real
>() + 3*atan(root_two
<Real
>())/root_two
<Real
>();
274 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
275 BOOST_CHECK_CLOSE_FRACTION(L1
, Q_expected
, tol
);
277 auto f5
= [](Real t
)->Real
{ return t
*t
*log(t
)/((t
*t
- 1)*(t
*t
*t
*t
+ 1)); };
278 Q
= integrator
.integrate(f5
, (Real
) 0 , (Real
) 1);
279 Q_expected
= pi
<Real
>()*pi
<Real
>()*(2 - root_two
<Real
>())/32;
280 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
283 // Oh it suffers on this one.
284 auto f6
= [](Real t
)->Real
{ Real x
= log(t
); return x
*x
; };
285 Q
= integrator
.integrate(f6
, (Real
) 0 , (Real
) 1, get_convergence_tolerance
<Real
>(), &error
, &L1
);
288 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, 50*tol
);
289 BOOST_CHECK_CLOSE_FRACTION(L1
, Q_expected
, 50*tol
);
292 // Although it doesn't get to the requested tolerance on this integral, the error bounds can be queried and are reasonable:
293 tol
= sqrt(boost::math::tools::epsilon
<Real
>());
294 auto f7
= [](const Real
& t
) { return sqrt(tan(t
)); };
295 Q
= integrator
.integrate(f7
, (Real
) 0 , (Real
) half_pi
<Real
>(), get_convergence_tolerance
<Real
>(), &error
, &L1
);
296 Q_expected
= pi
<Real
>()/root_two
<Real
>();
297 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
298 BOOST_CHECK_CLOSE_FRACTION(L1
, Q_expected
, tol
);
300 auto f8
= [](const Real
& t
) { return log(cos(t
)); };
301 Q
= integrator
.integrate(f8
, (Real
) 0 , half_pi
<Real
>(), get_convergence_tolerance
<Real
>(), &error
, &L1
);
302 Q_expected
= -pi
<Real
>()*ln_two
<Real
>()*half
<Real
>();
303 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
304 BOOST_CHECK_CLOSE_FRACTION(L1
, -Q_expected
, tol
);
309 void test_three_quadrature_schemes_examples()
311 std::cout
<< "Testing integral in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id
<Real
>().pretty_name() << "\n";
312 Real tol
= 10 * boost::math::tools::epsilon
<Real
>();
316 auto integrator
= get_integrator
<Real
>();
318 auto f1
= [](const Real
& t
) { return t
*boost::math::log1p(t
); };
319 Q
= integrator
.integrate(f1
, (Real
) 0 , (Real
) 1);
320 Q_expected
= half
<Real
>()*half
<Real
>();
321 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
325 auto f2
= [](const Real
& t
) { return t
*t
*atan(t
); };
326 Q
= integrator
.integrate(f2
, (Real
) 0 , (Real
) 1);
327 Q_expected
= (pi
<Real
>() -2 + 2*ln_two
<Real
>())/12;
328 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, 2 * tol
);
331 auto f3
= [](const Real
& t
) { return exp(t
)*cos(t
); };
332 Q
= integrator
.integrate(f3
, (Real
) 0, half_pi
<Real
>());
333 Q_expected
= boost::math::expm1(half_pi
<Real
>())*half
<Real
>();
334 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
337 auto f4
= [](Real x
)->Real
{ Real t0
= sqrt(x
*x
+ 2); return atan(t0
)/(t0
*(x
*x
+1)); };
338 Q
= integrator
.integrate(f4
, (Real
) 0 , (Real
) 1);
339 Q_expected
= 5*pi
<Real
>()*pi
<Real
>()/96;
340 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
343 auto f5
= [](const Real
& t
) { return sqrt(t
)*log(t
); };
344 Q
= integrator
.integrate(f5
, (Real
) 0 , (Real
) 1);
345 Q_expected
= -4/ (Real
) 9;
346 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
349 auto f6
= [](const Real
& t
) { return sqrt(1 - t
*t
); };
350 Q
= integrator
.integrate(f6
, (Real
) 0 , (Real
) 1);
351 Q_expected
= pi
<Real
>()/4;
352 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
357 void test_integration_over_real_line()
359 std::cout
<< "Testing integrals over entire real line in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id
<Real
>().pretty_name() << "\n";
360 Real tol
= 10 * boost::math::tools::epsilon
<Real
>();
365 auto integrator
= get_integrator
<Real
>();
367 auto f1
= [](const Real
& t
) { return 1/(1+t
*t
);};
368 Q
= integrator
.integrate(f1
, std::numeric_limits
<Real
>::has_infinity
? -std::numeric_limits
<Real
>::infinity() : -boost::math::tools::max_value
<Real
>(), std::numeric_limits
<Real
>::has_infinity
? std::numeric_limits
<Real
>::infinity() : boost::math::tools::max_value
<Real
>(), get_convergence_tolerance
<Real
>(), &error
, &L1
);
369 Q_expected
= pi
<Real
>();
370 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
371 BOOST_CHECK_CLOSE_FRACTION(L1
, Q_expected
, tol
);
373 auto f2
= [](const Real
& t
) { return exp(-t
*t
*half
<Real
>()); };
374 Q
= integrator
.integrate(f2
, std::numeric_limits
<Real
>::has_infinity
? -std::numeric_limits
<Real
>::infinity() : -boost::math::tools::max_value
<Real
>(), std::numeric_limits
<Real
>::has_infinity
? std::numeric_limits
<Real
>::infinity() : boost::math::tools::max_value
<Real
>(), get_convergence_tolerance
<Real
>(), &error
, &L1
);
375 Q_expected
= root_two_pi
<Real
>();
376 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
* 2);
377 BOOST_CHECK_CLOSE_FRACTION(L1
, Q_expected
, tol
* 2);
379 // This test shows how oscillatory integrals are approximated very poorly by this method:
380 //std::cout << "Testing sinc integral: \n";
381 //Q = integrator.integrate(boost::math::sinc_pi<Real>, -std::numeric_limits<Real>::infinity(), std::numeric_limits<Real>::infinity(), &error, &L1);
382 //std::cout << "Error estimate of sinc integral: " << error << std::endl;
383 //std::cout << "L1 norm of sinc integral " << L1 << std::endl;
384 //Q_expected = pi<Real>();
385 //BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
387 auto f4
= [](const Real
& t
) { return 1/cosh(t
);};
388 Q
= integrator
.integrate(f4
, std::numeric_limits
<Real
>::has_infinity
? -std::numeric_limits
<Real
>::infinity() : -boost::math::tools::max_value
<Real
>(), std::numeric_limits
<Real
>::has_infinity
? std::numeric_limits
<Real
>::infinity() : boost::math::tools::max_value
<Real
>(), get_convergence_tolerance
<Real
>(), &error
, &L1
);
389 Q_expected
= pi
<Real
>();
390 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
391 BOOST_CHECK_CLOSE_FRACTION(L1
, Q_expected
, tol
);
396 void test_right_limit_infinite()
398 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";
399 Real tol
= 10 * boost::math::tools::epsilon
<Real
>();
404 auto integrator
= get_integrator
<Real
>();
407 auto f1
= [](const Real
& t
) { return 1/(1+t
*t
);};
408 Q
= integrator
.integrate(f1
, 0, std::numeric_limits
<Real
>::has_infinity
? std::numeric_limits
<Real
>::infinity() : boost::math::tools::max_value
<Real
>(), get_convergence_tolerance
<Real
>(), &error
, &L1
);
409 Q_expected
= half_pi
<Real
>();
410 BOOST_CHECK_CLOSE(Q
, Q_expected
, 100*tol
);
413 auto f2
= [](const Real
& t
) { return exp(-t
)/sqrt(t
); };
414 Q
= integrator
.integrate(f2
, 0, std::numeric_limits
<Real
>::has_infinity
? std::numeric_limits
<Real
>::infinity() : boost::math::tools::max_value
<Real
>(), get_convergence_tolerance
<Real
>(), &error
, &L1
);
415 Q_expected
= root_pi
<Real
>();
416 BOOST_CHECK_CLOSE(Q
, Q_expected
, 1000*tol
);
418 auto f3
= [](const Real
& t
) { return exp(-t
)*cos(t
); };
419 Q
= integrator
.integrate(f3
, 0, std::numeric_limits
<Real
>::has_infinity
? std::numeric_limits
<Real
>::infinity() : boost::math::tools::max_value
<Real
>(), get_convergence_tolerance
<Real
>(), &error
, &L1
);
420 Q_expected
= half
<Real
>();
421 BOOST_CHECK_CLOSE(Q
, Q_expected
, 100*tol
);
423 auto f4
= [](const Real
& t
) { return 1/(1+t
*t
); };
424 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
);
425 Q_expected
= pi
<Real
>()/4;
426 BOOST_CHECK_CLOSE(Q
, Q_expected
, 100*tol
);
430 void test_left_limit_infinite()
432 std::cout
<< "Testing left limit infinite for tanh_sinh in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id
<Real
>().pretty_name() << "\n";
433 Real tol
= 10 * boost::math::tools::epsilon
<Real
>();
436 auto integrator
= get_integrator
<Real
>();
439 auto f1
= [](const Real
& t
) { return 1/(1+t
*t
);};
440 Q
= integrator
.integrate(f1
, std::numeric_limits
<Real
>::has_infinity
? -std::numeric_limits
<Real
>::infinity() : -boost::math::tools::max_value
<Real
>(), Real(0));
441 Q_expected
= half_pi
<Real
>();
442 BOOST_CHECK_CLOSE(Q
, Q_expected
, 100*tol
);
446 // A horrible function taken from
447 // http://www.chebfun.org/examples/quad/GaussClenCurt.html
451 std::cout
<< "Testing Trefenthen's horrible integral on type " << boost::typeindex::type_id
<Real
>().pretty_name() << "\n";
452 // We only know the integral to double precision, so requesting a higher tolerance doesn't make sense.
453 Real tol
= 10 * std::numeric_limits
<float>::epsilon();
458 auto integrator
= get_integrator
<Real
>();
460 auto f
= [](Real x
)->Real
{ return x
*sin(2*exp(2*sin(2*exp(2*x
) ) ) ); };
461 Q
= integrator
.integrate(f
, (Real
) -1, (Real
) 1, get_convergence_tolerance
<Real
>(), &error
, &L1
);
462 // NIntegrate[x*Sin[2*Exp[2*Sin[2*Exp[2*x]]]], {x, -1, 1}, WorkingPrecision -> 130, MaxRecursion -> 100]
463 Q_expected
= boost::lexical_cast
<Real
>("0.33673283478172753598559003181355241139806404130031017259552729882281");
464 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
465 // Over again without specifying the bounds:
466 Q
= integrator
.integrate(f
, get_convergence_tolerance
<Real
>(), &error
, &L1
);
467 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
470 // Some examples of tough integrals from NR, section 4.5.4:
472 void test_nr_examples()
474 std::cout
<< "Testing singular integrals from NR 4.5.4 on type " << boost::typeindex::type_id
<Real
>().pretty_name() << "\n";
475 Real tol
= 10 * boost::math::tools::epsilon
<Real
>();
480 auto integrator
= get_integrator
<Real
>();
482 auto f1
= [](Real x
)->Real
484 return (sin(x
* half
<Real
>()) * exp(-x
) / x
) / sqrt(x
);
486 Q
= integrator
.integrate(f1
, 0, std::numeric_limits
<Real
>::has_infinity
? std::numeric_limits
<Real
>::infinity() : boost::math::tools::max_value
<Real
>(), get_convergence_tolerance
<Real
>(), &error
, &L1
);
487 Q_expected
= sqrt(pi
<Real
>()*(sqrt((Real
) 5) - 2));
488 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, 25*tol
);
490 auto f2
= [](Real x
)->Real
{ return pow(x
, -(Real
) 2/(Real
) 7)*exp(-x
*x
); };
491 Q
= integrator
.integrate(f2
, 0, std::numeric_limits
<Real
>::has_infinity
? std::numeric_limits
<Real
>::infinity() : boost::math::tools::max_value
<Real
>());
492 Q_expected
= half
<Real
>()*boost::math::tgamma((Real
) 5/ (Real
) 14);
493 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
* 3);
497 // Test integrand known to fool some termination schemes:
499 void test_early_termination()
501 std::cout
<< "Testing Clenshaw & Curtis's example of integrand which fools termination schemes on type " << boost::typeindex::type_id
<Real
>().pretty_name() << "\n";
502 Real tol
= 10 * boost::math::tools::epsilon
<Real
>();
507 auto integrator
= get_integrator
<Real
>();
509 auto f1
= [](Real x
)->Real
{ return 23*cosh(x
)/ (Real
) 25 - cos(x
) ; };
510 Q
= integrator
.integrate(f1
, (Real
) -1, (Real
) 1, get_convergence_tolerance
<Real
>(), &error
, &L1
);
511 Q_expected
= 46*sinh((Real
) 1)/(Real
) 25 - 2*sin((Real
) 1);
512 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
513 // Over again with no bounds:
514 Q
= integrator
.integrate(f1
);
515 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
518 // Test some definite integrals from the CRC handbook, 32nd edition:
522 std::cout
<< "Testing CRC formulas on type " << boost::typeindex::type_id
<Real
>().pretty_name() << "\n";
523 Real tol
= 10 * boost::math::tools::epsilon
<Real
>();
528 auto integrator
= get_integrator
<Real
>();
530 // CRC Definite integral 585
531 auto f1
= [](Real x
)->Real
{ Real t
= log(1/x
); return x
*x
*t
*t
*t
; };
532 Q
= integrator
.integrate(f1
, (Real
) 0, (Real
) 1, get_convergence_tolerance
<Real
>(), &error
, &L1
);
533 Q_expected
= (Real
) 2/ (Real
) 27;
534 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
537 auto f2
= [](Real x
)->Real
{ return sqrt(cos(x
)); };
538 Q
= integrator
.integrate(f2
, (Real
) 0, (Real
) half_pi
<Real
>(), get_convergence_tolerance
<Real
>(), &error
, &L1
);
539 //Q_expected = pow(two_pi<Real>(), 3*half<Real>())/pow(boost::math::tgamma((Real) 1/ (Real) 4), 2);
540 Q_expected
= boost::lexical_cast
<Real
>("1.198140234735592207439922492280323878227212663215651558263674952946405214143915670835885556489793389375907225");
541 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
543 // CRC Section 5.5, integral 585:
544 for (int n
= 0; n
< 3; ++n
) {
545 for (int m
= 0; m
< 3; ++m
) {
546 auto f
= [&](Real x
)->Real
{ return pow(x
, Real(m
))*pow(log(1/x
), Real(n
)); };
547 Q
= integrator
.integrate(f
, (Real
) 0, (Real
) 1, get_convergence_tolerance
<Real
>(), &error
, &L1
);
548 // Calculation of the tgamma function is not exact, giving spurious failures.
549 // Casting to cpp_bin_float_100 beforehand fixes most of them.
550 cpp_bin_float_100 np1
= n
+ 1;
551 cpp_bin_float_100 mp1
= m
+ 1;
552 Q_expected
= boost::lexical_cast
<Real
>((tgamma(np1
)/pow(mp1
, np1
)).str());
553 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
557 // CRC Section 5.5, integral 591
558 // The parameter p allows us to control the strength of the singularity.
559 // Rapid convergence is not guaranteed for this function, as the branch cut makes it non-analytic on a disk.
560 // This converges only when our test type has an extended exponent range as all the area of the integral
561 // occurs so close to 0 (or 1) that we need abscissa values exceptionally small to find it.
562 // "There's a lot of room at the bottom".
563 // We also use a 2 argument functor so that 1-x is evaluated accurately:
564 if (std::numeric_limits
<Real
>::max_exponent
> std::numeric_limits
<double>::max_exponent
)
566 for (Real p
= -0.99; p
< 1; p
+= 0.1) {
567 auto f
= [&](Real x
, Real cx
)->Real
569 //return pow(x, p) / pow(1 - x, p);
570 return cx
< 0 ? exp(p
* (log(x
) - boost::math::log1p(-x
))) : pow(x
, p
) / pow(cx
, p
);
572 Q
= integrator
.integrate(f
, (Real
)0, (Real
)1, get_convergence_tolerance
<Real
>(), &error
, &L1
);
573 Q_expected
= 1 / sinc_pi(p
*pi
<Real
>());
574 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, 10 * tol
);
577 // There is an alternative way to evaluate the above integral: by noticing that all the area of the integral
578 // is near zero for p < 0 and near 1 for p > 0 we can substitute exp(-x) for x and remap the integral to the
579 // domain (0, INF). Internally we need to expand out the terms and evaluate using logs to avoid spurous overflow,
582 for (Real p
= 0.99; p
> 0; p
-= 0.1) {
583 auto f
= [&](Real x
)->Real
585 return exp(-x
* (1 - p
) + p
* log(-boost::math::expm1(-x
)));
587 Q
= integrator
.integrate(f
, 0, boost::math::tools::max_value
<Real
>(), get_convergence_tolerance
<Real
>(), &error
, &L1
);
588 Q_expected
= 1 / sinc_pi(p
*pi
<Real
>());
589 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, 10 * tol
);
592 for (Real p
= -0.99; p
< 0; p
+= 0.1) {
593 auto f
= [&](Real x
)->Real
595 return exp(-p
* log(-boost::math::expm1(-x
)) - (1 + p
) * x
);
597 Q
= integrator
.integrate(f
, 0, boost::math::tools::max_value
<Real
>(), get_convergence_tolerance
<Real
>(), &error
, &L1
);
598 Q_expected
= 1 / sinc_pi(p
*pi
<Real
>());
599 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, 10 * tol
);
602 // CRC Section 5.5, integral 635
603 for (int m
= 0; m
< 10; ++m
) {
604 auto f
= [&](Real x
)->Real
{ return 1/(1 + pow(tan(x
), m
)); };
605 Q
= integrator
.integrate(f
, (Real
) 0, half_pi
<Real
>(), get_convergence_tolerance
<Real
>(), &error
, &L1
);
606 Q_expected
= half_pi
<Real
>()/2;
607 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
610 // CRC Section 5.5, integral 637:
612 // When h gets very close to 1, the strength of the singularity gradually increases until we
613 // no longer have enough exponent range to evaluate the integral....
614 // We also have to use the 2-argument functor version of the integrator to avoid
615 // cancellation error, since the singularity is near PI/2.
617 Real limit
= std::numeric_limits
<Real
>::max_exponent
> std::numeric_limits
<double>::max_exponent
619 for (Real h
= 0.01; h
< limit
; h
+= 0.1) {
620 auto f
= [&](Real x
, Real xc
)->Real
{ return xc
> 0 ? pow(1/tan(xc
), h
) : pow(tan(x
), h
); };
621 Q
= integrator
.integrate(f
, (Real
) 0, half_pi
<Real
>(), get_convergence_tolerance
<Real
>(), &error
, &L1
);
622 Q_expected
= half_pi
<Real
>()/cos(h
*half_pi
<Real
>());
623 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
625 // CRC Section 5.5, integral 637:
627 // Over again, but with argument transformation, we want:
629 // Integral of tan(x)^h over (0, PI/2)
631 // Note that the bulk of the area is next to the singularity at PI/2,
632 // so we'll start by replacing x by PI/2 - x, and that tan(PI/2 - x) == 1/tan(x)
635 // Integral of 1/tan(x)^h over (0, PI/2)
637 // Which is almost the ideal form, except that when h is very close to 1
638 // we run out of exponent range in evaluating the integral arbitrarily close to 0.
639 // So now we substitute exp(-x) for x: this stretches out the range of the
640 // integral to (-log(PI/2), +INF) with the singularity at +INF giving:
642 // Integral of exp(-x)/tan(exp(-x))^h over (-log(PI/2), +INF)
644 // We just need a way to evaluate the function without spurious under/overflow
645 // in the exp terms. Note that for small x: tan(x) ~= x, so making this
646 // substitution and evaluating by logs we have:
648 // exp(-x)/tan(exp(-x))^h ~= exp((h - 1) * x) for x > -log(epsion);
650 // Here's how that looks in code:
652 for (Real i
= 80; i
< 100; ++i
) {
654 auto f
= [&](Real x
)->Real
656 if (x
> -log(boost::math::tools::epsilon
<Real
>()))
657 return exp((h
- 1) * x
);
661 // Need to deal with numeric instability causing et to be greater than PI/2:
662 return et
>= boost::math::constants::half_pi
<Real
>() ? 0 : et
* pow(1 / tan(et
), h
);
665 Q
= integrator
.integrate(f
, -log(half_pi
<Real
>()), boost::math::tools::max_value
<Real
>(), get_convergence_tolerance
<Real
>(), &error
, &L1
);
666 Q_expected
= half_pi
<Real
>() / cos(h
*half_pi
<Real
>());
667 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, 5 * tol
);
670 // CRC Section 5.5, integral 670:
672 auto f3
= [](Real x
)->Real
{ return sqrt(log(1/x
)); };
673 Q
= integrator
.integrate(f3
, (Real
) 0, (Real
) 1, get_convergence_tolerance
<Real
>(), &error
, &L1
);
674 Q_expected
= root_pi
<Real
>()/2;
675 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
679 template <class Real
>
683 // Test some special functions that we already know how to evaluate:
684 std::cout
<< "Testing special functions on type " << boost::typeindex::type_id
<Real
>().pretty_name() << "\n";
685 Real tol
= 10 * boost::math::tools::epsilon
<Real
>();
686 auto integrator
= get_integrator
<Real
>();
689 if (std::numeric_limits
<Real
>::digits10
< 37) // Otherwise too slow
692 auto f
= [&](Real x
)->Real
{ return boost::math::ibeta_derivative(a
, b
, x
); };
693 BOOST_CHECK_CLOSE_FRACTION(integrator
.integrate(f
, 0, Real(0.25)), boost::math::ibeta(100, 15, Real(0.25)), tol
* 10);
694 // Check some really extreme versions:
697 BOOST_CHECK_CLOSE_FRACTION(integrator
.integrate(f
, 0, 1), Real(1), tol
* 10);
699 // This is as extreme as we can get in this domain: otherwise the function has all it's
700 // area so close to zero we never get in there no matter how many levels we go down:
704 BOOST_CHECK_CLOSE_FRACTION(integrator
.integrate(f
, 0, 1), Real(1), tol
* 25);
707 Real x
= 2, y
= 3, z
= 0.5, p
= 0.25;
708 // Elliptic integral RC:
709 Real Q
= integrator
.integrate([&](const Real
& t
) { return 1 / (sqrt(t
+ x
) * (t
+ y
)); }, 0, std::numeric_limits
<Real
>::has_infinity
? std::numeric_limits
<Real
>::infinity() : boost::math::tools::max_value
<Real
>()) / 2;
710 BOOST_CHECK_CLOSE_FRACTION(Q
, boost::math::ellint_rc(x
, y
), tol
);
711 // Elliptic Integral RJ:
712 BOOST_CHECK_CLOSE_FRACTION(Real((Real(3) / 2) * integrator
.integrate([&](Real t
)->Real
{ return 1 / (sqrt((t
+ x
) * (t
+ y
) * (t
+ z
)) * (t
+ p
)); }, 0, std::numeric_limits
<Real
>::has_infinity
? std::numeric_limits
<Real
>::infinity() : boost::math::tools::max_value
<Real
>())), boost::math::ellint_rj(x
, y
, z
, p
), tol
);
715 if (std::numeric_limits
<Real
>::digits10
> 40)
717 else if (!std::numeric_limits
<Real
>::is_specialized
)
719 // tgamma expressed as an integral:
720 BOOST_CHECK_CLOSE_FRACTION(integrator
.integrate([&](Real t
)->Real
{ using std::pow
; using std::exp
; return t
> 10000 ? Real(0) : Real(pow(t
, z
- 1) * exp(-t
)); }, 0, std::numeric_limits
<Real
>::has_infinity
? std::numeric_limits
<Real
>::infinity() : boost::math::tools::max_value
<Real
>()),
721 boost::math::tgamma(z
), tol
);
722 BOOST_CHECK_CLOSE_FRACTION(integrator
.integrate([](const Real
& t
) { using std::exp
; return exp(-t
*t
); }, std::numeric_limits
<Real
>::has_infinity
? -std::numeric_limits
<Real
>::infinity() : -boost::math::tools::max_value
<Real
>(), std::numeric_limits
<Real
>::has_infinity
? std::numeric_limits
<Real
>::infinity() : boost::math::tools::max_value
<Real
>()),
723 boost::math::constants::root_pi
<Real
>(), tol
);
726 template <class Real
>
730 std::cout
<< "Testing 2 argument functors on type " << boost::typeindex::type_id
<Real
>().pretty_name() << "\n";
731 Real tol
= 10 * boost::math::tools::epsilon
<Real
>();
733 auto integrator
= get_integrator
<Real
>();
736 // There are a whole family of integrals of the general form
738 // which have all the interesting behaviour near the 2 singularities
739 // and all converge, see: http://www.wolframalpha.com/input/?i=integrate+(x+*+(1-x))%5E-1%2FN+from+0+to+1
741 Real Q
= integrator
.integrate([&](const Real
& t
, const Real
& tc
)->Real
743 return tc
< 0 ? 1 / sqrt(t
* (1-t
)) : 1 / sqrt(t
* tc
);
745 BOOST_CHECK_CLOSE_FRACTION(Q
, boost::math::constants::pi
<Real
>(), tol
);
746 Q
= integrator
.integrate([&](const Real
& t
, const Real
& tc
)->Real
748 return tc
< 0 ? 1 / boost::math::cbrt(t
* (1-t
)) : 1 / boost::math::cbrt(t
* tc
);
750 BOOST_CHECK_CLOSE_FRACTION(Q
, boost::math::pow
<2>(boost::math::tgamma(Real(2) / 3)) / boost::math::tgamma(Real(4) / 3), tol
* 3);
752 // We can do the same thing with ((1+x)(1-x))^-N ; N < 1
754 Q
= integrator
.integrate([&](const Real
& t
, const Real
& tc
)->Real
756 return t
< 0 ? 1 / sqrt(-tc
* (1-t
)) : 1 / sqrt((t
+ 1) * tc
);
758 BOOST_CHECK_CLOSE_FRACTION(Q
, boost::math::constants::pi
<Real
>(), tol
);
759 Q
= integrator
.integrate([&](const Real
& t
, const Real
& tc
)->Real
761 return t
< 0 ? 1 / sqrt(-tc
* (1-t
)) : 1 / sqrt((t
+ 1) * tc
);
763 BOOST_CHECK_CLOSE_FRACTION(Q
, boost::math::constants::pi
<Real
>(), tol
);
764 Q
= integrator
.integrate([&](const Real
& t
, const Real
& tc
)->Real
766 return t
< 0 ? 1 / boost::math::cbrt(-tc
* (1-t
)) : 1 / boost::math::cbrt((t
+ 1) * tc
);
768 BOOST_CHECK_CLOSE_FRACTION(Q
, sqrt(boost::math::constants::pi
<Real
>()) * boost::math::tgamma(Real(2) / 3) / boost::math::tgamma(Real(7) / 6), tol
* 10);
769 Q
= integrator
.integrate([&](const Real
& t
, const Real
& tc
)->Real
771 return t
< 0 ? 1 / boost::math::cbrt(-tc
* (1-t
)) : 1 / boost::math::cbrt((t
+ 1) * tc
);
773 BOOST_CHECK_CLOSE_FRACTION(Q
, sqrt(boost::math::constants::pi
<Real
>()) * boost::math::tgamma(Real(2) / 3) / boost::math::tgamma(Real(7) / 6), tol
* 10);
775 // These are taken from above, and do not get to full precision via the single arg functor:
777 auto f7
= [](const Real
& t
, const Real
& tc
) { return t
< 1 ? sqrt(tan(t
)) : sqrt(1/tan(tc
)); };
779 Q
= integrator
.integrate(f7
, (Real
)0, (Real
)half_pi
<Real
>(), get_convergence_tolerance
<Real
>(), &error
, &L1
);
780 Real Q_expected
= pi
<Real
>() / root_two
<Real
>();
781 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
782 BOOST_CHECK_CLOSE_FRACTION(L1
, Q_expected
, tol
);
784 auto f8
= [](const Real
& t
, const Real
& tc
) { return t
< 1 ? log(cos(t
)) : log(sin(tc
)); };
785 Q
= integrator
.integrate(f8
, (Real
)0, half_pi
<Real
>(), get_convergence_tolerance
<Real
>(), &error
, &L1
);
786 Q_expected
= -pi
<Real
>()*ln_two
<Real
>()*half
<Real
>();
787 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
788 BOOST_CHECK_CLOSE_FRACTION(L1
, -Q_expected
, tol
);
792 BOOST_AUTO_TEST_CASE(tanh_sinh_quadrature_test
)
794 #ifdef GENERATE_CONSTANTS
796 // Generate pre-computed coefficients:
797 std::cout
<< std::setprecision(35);
798 print_levels(generate_constants
<cpp_bin_float_100
>(10, 8), "L");
804 test_right_limit_infinite
<float>();
805 test_left_limit_infinite
<float>();
806 test_linear
<float>();
807 test_quadratic
<float>();
808 test_singular
<float>();
810 test_three_quadrature_schemes_examples
<float>();
811 test_horrible
<float>();
812 test_integration_over_real_line
<float>();
813 test_nr_examples
<float>();
816 test_early_termination
<float>();
822 test_right_limit_infinite
<double>();
823 test_left_limit_infinite
<double>();
824 test_linear
<double>();
825 test_quadratic
<double>();
826 test_singular
<double>();
828 test_three_quadrature_schemes_examples
<double>();
829 test_horrible
<double>();
830 test_integration_over_real_line
<double>();
831 test_nr_examples
<double>();
832 test_early_termination
<double>();
834 test_2_arg
<double>();
842 test_right_limit_infinite
<long double>();
843 test_left_limit_infinite
<long double>();
844 test_linear
<long double>();
845 test_quadratic
<long double>();
846 test_singular
<long double>();
847 test_ca
<long double>();
848 test_three_quadrature_schemes_examples
<long double>();
849 test_horrible
<long double>();
850 test_integration_over_real_line
<long double>();
851 test_nr_examples
<long double>();
852 test_early_termination
<long double>();
853 test_sf
<long double>();
854 test_2_arg
<long double>();
857 test_crc
<long double>();
863 test_right_limit_infinite
<cpp_bin_float_quad
>();
864 test_left_limit_infinite
<cpp_bin_float_quad
>();
865 test_linear
<cpp_bin_float_quad
>();
866 test_quadratic
<cpp_bin_float_quad
>();
867 test_singular
<cpp_bin_float_quad
>();
868 test_ca
<cpp_bin_float_quad
>();
869 test_three_quadrature_schemes_examples
<cpp_bin_float_quad
>();
870 test_horrible
<cpp_bin_float_quad
>();
871 test_nr_examples
<cpp_bin_float_quad
>();
872 test_early_termination
<cpp_bin_float_quad
>();
873 test_crc
<cpp_bin_float_quad
>();
874 test_sf
<cpp_bin_float_quad
>();
875 test_2_arg
<cpp_bin_float_quad
>();
880 test_sf
<cpp_bin_float_50
>();
881 test_sf
<cpp_bin_float_100
>();
882 test_sf
<boost::multiprecision::number
<boost::multiprecision::cpp_bin_float
<150> > >();
887 test_right_limit_infinite
<boost::math::concepts::real_concept
>();
888 test_left_limit_infinite
<boost::math::concepts::real_concept
>();
889 test_linear
<boost::math::concepts::real_concept
>();
890 test_quadratic
<boost::math::concepts::real_concept
>();
891 test_singular
<boost::math::concepts::real_concept
>();
892 test_ca
<boost::math::concepts::real_concept
>();
893 test_three_quadrature_schemes_examples
<boost::math::concepts::real_concept
>();
894 test_horrible
<boost::math::concepts::real_concept
>();
895 test_integration_over_real_line
<boost::math::concepts::real_concept
>();
896 test_nr_examples
<boost::math::concepts::real_concept
>();
897 test_early_termination
<boost::math::concepts::real_concept
>();
898 test_sf
<boost::math::concepts::real_concept
>();
899 test_2_arg
<boost::math::concepts::real_concept
>();
902 test_crc
<boost::math::concepts::real_concept
>();
906 test_sf
<cpp_dec_float_50
>();
908 #if defined(TEST8) && defined(BOOST_HAS_FLOAT128)
910 test_right_limit_infinite
<boost::multiprecision::float128
>();
911 test_left_limit_infinite
<boost::multiprecision::float128
>();
912 test_linear
<boost::multiprecision::float128
>();
913 test_quadratic
<boost::multiprecision::float128
>();
914 test_singular
<boost::multiprecision::float128
>();
915 test_ca
<boost::multiprecision::float128
>();
916 test_three_quadrature_schemes_examples
<boost::multiprecision::float128
>();
917 test_horrible
<boost::multiprecision::float128
>();
918 test_integration_over_real_line
<boost::multiprecision::float128
>();
919 test_nr_examples
<boost::multiprecision::float128
>();
920 test_early_termination
<boost::multiprecision::float128
>();
921 test_crc
<boost::multiprecision::float128
>();
922 test_sf
<boost::multiprecision::float128
>();
923 test_2_arg
<boost::multiprecision::float128
>();
931 int main() { return 0; }