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 tanh_sinh_quadrature_test
8
9 #include <boost/config.hpp>
10 #include <boost/detail/workaround.hpp>
11
12 #if !defined(BOOST_NO_CXX11_DECLTYPE) && !defined(BOOST_NO_CXX11_TRAILING_RESULT_TYPES) && !defined(BOOST_NO_SFINAE_EXPR)
13
14 #include <boost/math/concepts/real_concept.hpp>
15 #include <boost/test/included/unit_test.hpp>
16 #include <boost/test/tools/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>
26
27 #ifdef BOOST_HAS_FLOAT128
28 #include <boost/multiprecision/float128.hpp>
29 #endif
30
31 #ifdef _MSC_VER
32 #pragma warning(disable:4127) // Conditional expression is constant
33 #endif
34
35 #if !defined(TEST1) && !defined(TEST2) && !defined(TEST3) && !defined(TEST4) && !defined(TEST5) && !defined(TEST6) && !defined(TEST7) && !defined(TEST8)\
36 && !defined(TEST1A) && !defined(TEST1B) && !defined(TEST2A) && !defined(TEST3A) && !defined(TEST6A) && !defined(TEST9)
37 # define TEST1
38 # define TEST2
39 # define TEST3
40 # define TEST4
41 # define TEST5
42 # define TEST6
43 # define TEST7
44 # define TEST8
45 # define TEST1A
46 # define TEST1B
47 # define TEST2A
48 # define TEST3A
49 # define TEST6A
50 # define TEST9
51 #endif
52
53 using std::expm1;
54 using std::atan;
55 using std::tan;
56 using std::log;
57 using std::log1p;
58 using std::asinh;
59 using std::atanh;
60 using std::sqrt;
61 using std::isnormal;
62 using std::abs;
63 using std::sinh;
64 using std::tanh;
65 using std::cosh;
66 using std::pow;
67 using std::exp;
68 using std::sin;
69 using std::cos;
70 using std::string;
71 using boost::multiprecision::cpp_bin_float_50;
72 using boost::multiprecision::cpp_bin_float_100;
73 using boost::multiprecision::cpp_dec_float_50;
74 using boost::multiprecision::cpp_dec_float_100;
75 using boost::multiprecision::cpp_bin_float_quad;
76 using boost::math::sinc_pi;
77 using boost::math::quadrature::tanh_sinh;
78 using boost::math::quadrature::detail::tanh_sinh_detail;
79 using boost::math::constants::pi;
80 using boost::math::constants::half_pi;
81 using boost::math::constants::two_div_pi;
82 using boost::math::constants::two_pi;
83 using boost::math::constants::half;
84 using boost::math::constants::third;
85 using boost::math::constants::half;
86 using boost::math::constants::third;
87 using boost::math::constants::catalan;
88 using boost::math::constants::ln_two;
89 using boost::math::constants::root_two;
90 using boost::math::constants::root_two_pi;
91 using boost::math::constants::root_pi;
92
93 template <class T>
print_levels(const T & v,const char * suffix)94 void print_levels(const T& v, const char* suffix)
95 {
96 std::cout << "{\n";
97 for (unsigned i = 0; i < v.size(); ++i)
98 {
99 std::cout << " { ";
100 for (unsigned j = 0; j < v[i].size(); ++j)
101 {
102 std::cout << v[i][j] << suffix << ", ";
103 }
104 std::cout << "},\n";
105 }
106 std::cout << " };\n";
107 }
108
109 template <class T>
print_complement_indexes(const T & v)110 void print_complement_indexes(const T& v)
111 {
112 std::cout << "\n {";
113 for (unsigned i = 0; i < v.size(); ++i)
114 {
115 unsigned index = 0;
116 while (v[i][index] >= 0)
117 ++index;
118 std::cout << index << ", ";
119 }
120 std::cout << "};\n";
121 }
122
123 template <class T>
print_levels(const std::pair<T,T> & p,const char * suffix="")124 void print_levels(const std::pair<T, T>& p, const char* suffix = "")
125 {
126 std::cout << " static const std::vector<std::vector<Real> > abscissa = ";
127 print_levels(p.first, suffix);
128 std::cout << " static const std::vector<std::vector<Real> > weights = ";
129 print_levels(p.second, suffix);
130 std::cout << " static const std::vector<unsigned > indexes = ";
131 print_complement_indexes(p.first);
132 }
133
134 template <class Real>
generate_constants(unsigned max_index,unsigned max_rows)135 std::pair<std::vector<std::vector<Real>>, std::vector<std::vector<Real>> > generate_constants(unsigned max_index, unsigned max_rows)
136 {
137 using boost::math::constants::half_pi;
138 using boost::math::constants::two_div_pi;
139 using boost::math::constants::pi;
140 auto g = [](Real t) { return tanh(half_pi<Real>()*sinh(t)); };
141 auto w = [](Real t) { Real cs = cosh(half_pi<Real>() * sinh(t)); return half_pi<Real>() * cosh(t) / (cs * cs); };
142 auto gc = [](Real t) { Real u2 = half_pi<Real>() * sinh(t); return 1 / (exp(u2) *cosh(u2)); };
143 auto g_inv = [](float x)->float { return asinh(two_div_pi<float>()*atanh(x)); };
144 auto gc_inv = [](float x)
145 {
146 float l = log(sqrt((2 - x) / x));
147 return log((sqrt(4 * l * l + pi<float>() * pi<float>()) + 2 * l) / pi<float>());
148 };
149
150 std::vector<std::vector<Real>> abscissa, weights;
151
152 std::vector<Real> temp;
153
154 float t_crossover = gc_inv(0.5f);
155
156 Real h = 1;
157 for (unsigned i = 0; i < max_index; ++i)
158 {
159 temp.push_back((i < t_crossover) ? g(i * h) : -gc(i * h));
160 }
161 abscissa.push_back(temp);
162 temp.clear();
163
164 for (unsigned i = 0; i < max_index; ++i)
165 {
166 temp.push_back(w(i * h));
167 }
168 weights.push_back(temp);
169 temp.clear();
170
171 for (unsigned row = 1; row < max_rows; ++row)
172 {
173 h /= 2;
174 for (Real t = h; t < max_index - 1; t += 2 * h)
175 temp.push_back((t < t_crossover) ? g(t) : -gc(t));
176 abscissa.push_back(temp);
177 temp.clear();
178 }
179 h = 1;
180 for (unsigned row = 1; row < max_rows; ++row)
181 {
182 h /= 2;
183 for (Real t = h; t < max_index - 1; t += 2 * h)
184 temp.push_back(w(t));
185 weights.push_back(temp);
186 temp.clear();
187 }
188
189 return std::make_pair(abscissa, weights);
190 }
191
192 template <class Real>
get_integrator()193 const tanh_sinh<Real>& get_integrator()
194 {
195 static const tanh_sinh<Real> integrator(15);
196 return integrator;
197 }
198
199 template <class Real>
get_convergence_tolerance()200 Real get_convergence_tolerance()
201 {
202 return boost::math::tools::root_epsilon<Real>();
203 }
204
205
206 template<class Real>
test_linear()207 void test_linear()
208 {
209 std::cout << "Testing linear functions are integrated properly by tanh_sinh on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
210 Real tol = 10*boost::math::tools::epsilon<Real>();
211 auto integrator = get_integrator<Real>();
212 auto f = [](const Real& x)
213 {
214 return 5*x + 7;
215 };
216 Real error;
217 Real L1;
218 Real Q = integrator.integrate(f, (Real) 0, (Real) 1, get_convergence_tolerance<Real>(), &error, &L1);
219 BOOST_CHECK_CLOSE_FRACTION(Q, 9.5, tol);
220 BOOST_CHECK_CLOSE_FRACTION(L1, 9.5, tol);
221 Q = integrator.integrate(f, (Real) 1, (Real) 0, get_convergence_tolerance<Real>(), &error, &L1);
222 BOOST_CHECK_CLOSE_FRACTION(Q, -9.5, tol);
223 BOOST_CHECK_CLOSE_FRACTION(L1, 9.5, tol);
224 Q = integrator.integrate(f, (Real) 1, (Real) 1, get_convergence_tolerance<Real>(), &error, &L1);
225 BOOST_CHECK_EQUAL(Q, Real(0));
226 }
227
228
229 template<class Real>
test_quadratic()230 void test_quadratic()
231 {
232 std::cout << "Testing quadratic functions are integrated properly by tanh_sinh on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
233 Real tol = 10*boost::math::tools::epsilon<Real>();
234 auto integrator = get_integrator<Real>();
235 auto f = [](const Real& x) { return 5*x*x + 7*x + 12; };
236 Real error;
237 Real L1;
238 Real Q = integrator.integrate(f, (Real) 0, (Real) 1, get_convergence_tolerance<Real>(), &error, &L1);
239 BOOST_CHECK_CLOSE_FRACTION(Q, (Real) 17 + half<Real>()*third<Real>(), tol);
240 BOOST_CHECK_CLOSE_FRACTION(L1, (Real) 17 + half<Real>()*third<Real>(), tol);
241 }
242
243
244 template<class Real>
test_singular()245 void test_singular()
246 {
247 std::cout << "Testing integration of endpoint singularities on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
248 Real tol = 10*boost::math::tools::epsilon<Real>();
249 Real error;
250 Real L1;
251 auto integrator = get_integrator<Real>();
252 auto f = [](const Real& x) { return log(x)*log(1-x); };
253 Real Q = integrator.integrate(f, (Real) 0, (Real) 1, get_convergence_tolerance<Real>(), &error, &L1);
254 Real Q_expected = 2 - pi<Real>()*pi<Real>()*half<Real>()*third<Real>();
255
256 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
257 BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
258 }
259
260
261 // Examples taken from
262 //http://crd-legacy.lbl.gov/~dhbailey/dhbpapers/quadrature.pdf
263 template<class Real>
test_ca()264 void test_ca()
265 {
266 std::cout << "Testing integration of C(a) on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
267 Real tol = 10 * boost::math::tools::epsilon<Real>();
268 Real error;
269 Real L1;
270
271 auto integrator = get_integrator<Real>();
272 auto f1 = [](const Real& x) { return atan(x)/(x*(x*x + 1)) ; };
273 Real Q = integrator.integrate(f1, (Real) 0, (Real) 1, get_convergence_tolerance<Real>(), &error, &L1);
274 Real Q_expected = pi<Real>()*ln_two<Real>()/8 + catalan<Real>()*half<Real>();
275 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
276 BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
277
278 auto f2 = [](Real x)->Real { Real t0 = x*x + 1; Real t1 = sqrt(t0); return atan(t1)/(t0*t1); };
279 Q = integrator.integrate(f2, (Real) 0 , (Real) 1, get_convergence_tolerance<Real>(), &error, &L1);
280 Q_expected = pi<Real>()/4 - pi<Real>()/root_two<Real>() + 3*atan(root_two<Real>())/root_two<Real>();
281 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
282 BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
283
284 auto f5 = [](Real t)->Real { return t*t*log(t)/((t*t - 1)*(t*t*t*t + 1)); };
285 Q = integrator.integrate(f5, (Real) 0 , (Real) 1);
286 Q_expected = pi<Real>()*pi<Real>()*(2 - root_two<Real>())/32;
287 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
288
289
290 // Oh it suffers on this one.
291 auto f6 = [](Real t)->Real { Real x = log(t); return x*x; };
292 Q = integrator.integrate(f6, (Real) 0 , (Real) 1, get_convergence_tolerance<Real>(), &error, &L1);
293 Q_expected = 2;
294
295 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 50*tol);
296 BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, 50*tol);
297
298
299 // Although it doesn't get to the requested tolerance on this integral, the error bounds can be queried and are reasonable:
300 tol = sqrt(boost::math::tools::epsilon<Real>());
301 auto f7 = [](const Real& t) { return sqrt(tan(t)); };
302 Q = integrator.integrate(f7, (Real) 0 , (Real) half_pi<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
303 Q_expected = pi<Real>()/root_two<Real>();
304 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
305 BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
306
307 auto f8 = [](const Real& t) { return log(cos(t)); };
308 Q = integrator.integrate(f8, (Real) 0 , half_pi<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
309 Q_expected = -pi<Real>()*ln_two<Real>()*half<Real>();
310 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
311 BOOST_CHECK_CLOSE_FRACTION(L1, -Q_expected, tol);
312 }
313
314
315 template<class Real>
test_three_quadrature_schemes_examples()316 void test_three_quadrature_schemes_examples()
317 {
318 std::cout << "Testing integral in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
319 Real tol = 10 * boost::math::tools::epsilon<Real>();
320 Real Q;
321 Real Q_expected;
322
323 auto integrator = get_integrator<Real>();
324 // Example 1:
325 auto f1 = [](const Real& t) { return t*boost::math::log1p(t); };
326 Q = integrator.integrate(f1, (Real) 0 , (Real) 1);
327 Q_expected = half<Real>()*half<Real>();
328 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
329
330
331 // Example 2:
332 auto f2 = [](const Real& t) { return t*t*atan(t); };
333 Q = integrator.integrate(f2, (Real) 0 , (Real) 1);
334 Q_expected = (pi<Real>() -2 + 2*ln_two<Real>())/12;
335 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 2 * tol);
336
337 // Example 3:
338 auto f3 = [](const Real& t) { return exp(t)*cos(t); };
339 Q = integrator.integrate(f3, (Real) 0, half_pi<Real>());
340 Q_expected = boost::math::expm1(half_pi<Real>())*half<Real>();
341 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
342
343 // Example 4:
344 auto f4 = [](Real x)->Real { Real t0 = sqrt(x*x + 2); return atan(t0)/(t0*(x*x+1)); };
345 Q = integrator.integrate(f4, (Real) 0 , (Real) 1);
346 Q_expected = 5*pi<Real>()*pi<Real>()/96;
347 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
348
349 // Example 5:
350 auto f5 = [](const Real& t) { return sqrt(t)*log(t); };
351 Q = integrator.integrate(f5, (Real) 0 , (Real) 1);
352 Q_expected = -4/ (Real) 9;
353 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
354
355 // Example 6:
356 auto f6 = [](const Real& t) { return sqrt(1 - t*t); };
357 Q = integrator.integrate(f6, (Real) 0 , (Real) 1);
358 Q_expected = pi<Real>()/4;
359 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
360 }
361
362
363 template<class Real>
test_integration_over_real_line()364 void test_integration_over_real_line()
365 {
366 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";
367 Real tol = 10 * boost::math::tools::epsilon<Real>();
368 Real Q;
369 Real Q_expected;
370 Real error;
371 Real L1;
372 auto integrator = get_integrator<Real>();
373
374 auto f1 = [](const Real& t) { return 1/(1+t*t);};
375 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);
376 Q_expected = pi<Real>();
377 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
378 BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
379
380 auto f2 = [](const Real& t) { return exp(-t*t*half<Real>()); };
381 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);
382 Q_expected = root_two_pi<Real>();
383 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol * 2);
384 BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol * 2);
385
386 // This test shows how oscillatory integrals are approximated very poorly by this method:
387 //std::cout << "Testing sinc integral: \n";
388 //Q = integrator.integrate(boost::math::sinc_pi<Real>, -std::numeric_limits<Real>::infinity(), std::numeric_limits<Real>::infinity(), &error, &L1);
389 //std::cout << "Error estimate of sinc integral: " << error << std::endl;
390 //std::cout << "L1 norm of sinc integral " << L1 << std::endl;
391 //Q_expected = pi<Real>();
392 //BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
393
394 auto f4 = [](const Real& t) { return 1/cosh(t);};
395 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);
396 Q_expected = pi<Real>();
397 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
398 BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
399
400 }
401
402 template<class Real>
test_right_limit_infinite()403 void test_right_limit_infinite()
404 {
405 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";
406 Real tol = 10 * boost::math::tools::epsilon<Real>();
407 Real Q;
408 Real Q_expected;
409 Real error;
410 Real L1;
411 auto integrator = get_integrator<Real>();
412
413 // Example 11:
414 auto f1 = [](const Real& t) { return 1/(1+t*t);};
415 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);
416 Q_expected = half_pi<Real>();
417 BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
418
419 // Example 12
420 auto f2 = [](const Real& t) { return exp(-t)/sqrt(t); };
421 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);
422 Q_expected = root_pi<Real>();
423 BOOST_CHECK_CLOSE(Q, Q_expected, 1000*tol);
424
425 auto f3 = [](const Real& t) { return exp(-t)*cos(t); };
426 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);
427 Q_expected = half<Real>();
428 BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
429
430 auto f4 = [](const Real& t) { return 1/(1+t*t); };
431 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);
432 Q_expected = pi<Real>()/4;
433 BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
434 }
435
436 template<class Real>
test_left_limit_infinite()437 void test_left_limit_infinite()
438 {
439 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";
440 Real tol = 10 * boost::math::tools::epsilon<Real>();
441 Real Q;
442 Real Q_expected;
443 auto integrator = get_integrator<Real>();
444
445 // Example 11:
446 auto f1 = [](const Real& t) { return 1/(1+t*t);};
447 Q = integrator.integrate(f1, std::numeric_limits<Real>::has_infinity ? -std::numeric_limits<Real>::infinity() : -boost::math::tools::max_value<Real>(), Real(0));
448 Q_expected = half_pi<Real>();
449 BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
450 }
451
452
453 // A horrible function taken from
454 // http://www.chebfun.org/examples/quad/GaussClenCurt.html
455 template<class Real>
test_horrible()456 void test_horrible()
457 {
458 std::cout << "Testing Trefenthen's horrible integral on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
459 // We only know the integral to double precision, so requesting a higher tolerance doesn't make sense.
460 Real tol = 10 * std::numeric_limits<float>::epsilon();
461 Real Q;
462 Real Q_expected;
463 Real error;
464 Real L1;
465 auto integrator = get_integrator<Real>();
466
467 auto f = [](Real x)->Real { return x*sin(2*exp(2*sin(2*exp(2*x) ) ) ); };
468 Q = integrator.integrate(f, (Real) -1, (Real) 1, get_convergence_tolerance<Real>(), &error, &L1);
469 // NIntegrate[x*Sin[2*Exp[2*Sin[2*Exp[2*x]]]], {x, -1, 1}, WorkingPrecision -> 130, MaxRecursion -> 100]
470 Q_expected = boost::lexical_cast<Real>("0.33673283478172753598559003181355241139806404130031017259552729882281");
471 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
472 // Over again without specifying the bounds:
473 Q = integrator.integrate(f, get_convergence_tolerance<Real>(), &error, &L1);
474 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
475 }
476
477 // Some examples of tough integrals from NR, section 4.5.4:
478 template<class Real>
test_nr_examples()479 void test_nr_examples()
480 {
481 std::cout << "Testing singular integrals from NR 4.5.4 on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
482 Real tol = 10 * boost::math::tools::epsilon<Real>();
483 Real Q;
484 Real Q_expected;
485 Real error;
486 Real L1;
487 auto integrator = get_integrator<Real>();
488
489 auto f1 = [](Real x)->Real
490 {
491 return (sin(x * half<Real>()) * exp(-x) / x) / sqrt(x);
492 };
493 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);
494 Q_expected = sqrt(pi<Real>()*(sqrt((Real) 5) - 2));
495 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 25*tol);
496
497 auto f2 = [](Real x)->Real { return pow(x, -(Real) 2/(Real) 7)*exp(-x*x); };
498 Q = integrator.integrate(f2, 0, std::numeric_limits<Real>::has_infinity ? std::numeric_limits<Real>::infinity() : boost::math::tools::max_value<Real>());
499 Q_expected = half<Real>()*boost::math::tgamma((Real) 5/ (Real) 14);
500 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol * 6);
501
502 }
503
504 // Test integrand known to fool some termination schemes:
505 template<class Real>
test_early_termination()506 void test_early_termination()
507 {
508 std::cout << "Testing Clenshaw & Curtis's example of integrand which fools termination schemes on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
509 Real tol = 10 * boost::math::tools::epsilon<Real>();
510 Real Q;
511 Real Q_expected;
512 Real error;
513 Real L1;
514 auto integrator = get_integrator<Real>();
515
516 auto f1 = [](Real x)->Real { return 23*cosh(x)/ (Real) 25 - cos(x) ; };
517 Q = integrator.integrate(f1, (Real) -1, (Real) 1, get_convergence_tolerance<Real>(), &error, &L1);
518 Q_expected = 46*sinh((Real) 1)/(Real) 25 - 2*sin((Real) 1);
519 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
520 // Over again with no bounds:
521 Q = integrator.integrate(f1);
522 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
523 }
524
525 // Test some definite integrals from the CRC handbook, 32nd edition:
526 template<class Real>
test_crc()527 void test_crc()
528 {
529 std::cout << "Testing CRC formulas on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
530 Real tol = 10 * boost::math::tools::epsilon<Real>();
531 Real Q;
532 Real Q_expected;
533 Real error;
534 Real L1;
535 auto integrator = get_integrator<Real>();
536
537 // CRC Definite integral 585
538 auto f1 = [](Real x)->Real { Real t = log(1/x); return x*x*t*t*t; };
539 Q = integrator.integrate(f1, (Real) 0, (Real) 1, get_convergence_tolerance<Real>(), &error, &L1);
540 Q_expected = (Real) 2/ (Real) 27;
541 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
542
543 // CRC 636:
544 auto f2 = [](Real x)->Real { return sqrt(cos(x)); };
545 Q = integrator.integrate(f2, (Real) 0, (Real) half_pi<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
546 //Q_expected = pow(two_pi<Real>(), 3*half<Real>())/pow(boost::math::tgamma((Real) 1/ (Real) 4), 2);
547 Q_expected = boost::lexical_cast<Real>("1.198140234735592207439922492280323878227212663215651558263674952946405214143915670835885556489793389375907225");
548 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
549
550 // CRC Section 5.5, integral 585:
551 for (int n = 0; n < 3; ++n) {
552 for (int m = 0; m < 3; ++m) {
553 auto f = [&](Real x)->Real { return pow(x, Real(m))*pow(log(1/x), Real(n)); };
554 Q = integrator.integrate(f, (Real) 0, (Real) 1, get_convergence_tolerance<Real>(), &error, &L1);
555 // Calculation of the tgamma function is not exact, giving spurious failures.
556 // Casting to cpp_bin_float_100 beforehand fixes most of them.
557 cpp_bin_float_100 np1 = n + 1;
558 cpp_bin_float_100 mp1 = m + 1;
559 Q_expected = boost::lexical_cast<Real>((tgamma(np1)/pow(mp1, np1)).str());
560 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
561 }
562 }
563
564 // CRC Section 5.5, integral 591
565 // The parameter p allows us to control the strength of the singularity.
566 // Rapid convergence is not guaranteed for this function, as the branch cut makes it non-analytic on a disk.
567 // This converges only when our test type has an extended exponent range as all the area of the integral
568 // occurs so close to 0 (or 1) that we need abscissa values exceptionally small to find it.
569 // "There's a lot of room at the bottom".
570 // We also use a 2 argument functor so that 1-x is evaluated accurately:
571 if (std::numeric_limits<Real>::max_exponent > std::numeric_limits<double>::max_exponent)
572 {
573 for (Real p = Real (-0.99); p < 1; p += Real(0.1)) {
574 auto f = [&](Real x, Real cx)->Real
575 {
576 //return pow(x, p) / pow(1 - x, p);
577 return cx < 0 ? exp(p * (log(x) - boost::math::log1p(-x))) : pow(x, p) / pow(cx, p);
578 };
579 Q = integrator.integrate(f, (Real)0, (Real)1, get_convergence_tolerance<Real>(), &error, &L1);
580 Q_expected = 1 / sinc_pi(p*pi<Real>());
581 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 10 * tol);
582 }
583 }
584 // There is an alternative way to evaluate the above integral: by noticing that all the area of the integral
585 // 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
586 // domain (0, INF). Internally we need to expand out the terms and evaluate using logs to avoid spurious overflow,
587 // this gives us
588 // for p > 0:
589 for (Real p = Real(0.99); p > 0; p -= Real(0.1)) {
590 auto f = [&](Real x)->Real
591 {
592 return exp(-x * (1 - p) + p * log(-boost::math::expm1(-x)));
593 };
594 Q = integrator.integrate(f, 0, boost::math::tools::max_value<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
595 Q_expected = 1 / sinc_pi(p*pi<Real>());
596 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 10 * tol);
597 }
598 // and for p < 1:
599 for (Real p = Real (-0.99); p < 0; p += Real(0.1)) {
600 auto f = [&](Real x)->Real
601 {
602 return exp(-p * log(-boost::math::expm1(-x)) - (1 + p) * x);
603 };
604 Q = integrator.integrate(f, 0, boost::math::tools::max_value<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
605 Q_expected = 1 / sinc_pi(p*pi<Real>());
606 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 10 * tol);
607 }
608
609 // CRC Section 5.5, integral 635
610 for (int m = 0; m < 10; ++m) {
611 auto f = [&](Real x)->Real { return 1/(1 + pow(tan(x), m)); };
612 Q = integrator.integrate(f, (Real) 0, half_pi<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
613 Q_expected = half_pi<Real>()/2;
614 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
615 }
616
617 // CRC Section 5.5, integral 637:
618 //
619 // When h gets very close to 1, the strength of the singularity gradually increases until we
620 // no longer have enough exponent range to evaluate the integral....
621 // We also have to use the 2-argument functor version of the integrator to avoid
622 // cancellation error, since the singularity is near PI/2.
623 //
624 Real limit = std::numeric_limits<Real>::max_exponent > std::numeric_limits<double>::max_exponent
625 ? .95f : .85f;
626 for (Real h = Real(0.01); h < limit; h += Real(0.1)) {
627 auto f = [&](Real x, Real xc)->Real { return xc > 0 ? pow(1/tan(xc), h) : pow(tan(x), h); };
628 Q = integrator.integrate(f, (Real) 0, half_pi<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
629 Q_expected = half_pi<Real>()/cos(h*half_pi<Real>());
630 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
631 }
632 // CRC Section 5.5, integral 637:
633 //
634 // Over again, but with argument transformation, we want:
635 //
636 // Integral of tan(x)^h over (0, PI/2)
637 //
638 // Note that the bulk of the area is next to the singularity at PI/2,
639 // so we'll start by replacing x by PI/2 - x, and that tan(PI/2 - x) == 1/tan(x)
640 // so we now have:
641 //
642 // Integral of 1/tan(x)^h over (0, PI/2)
643 //
644 // Which is almost the ideal form, except that when h is very close to 1
645 // we run out of exponent range in evaluating the integral arbitrarily close to 0.
646 // So now we substitute exp(-x) for x: this stretches out the range of the
647 // integral to (-log(PI/2), +INF) with the singularity at +INF giving:
648 //
649 // Integral of exp(-x)/tan(exp(-x))^h over (-log(PI/2), +INF)
650 //
651 // We just need a way to evaluate the function without spurious under/overflow
652 // in the exp terms. Note that for small x: tan(x) ~= x, so making this
653 // substitution and evaluating by logs we have:
654 //
655 // exp(-x)/tan(exp(-x))^h ~= exp((h - 1) * x) for x > -log(epsilon);
656 //
657 // Here's how that looks in code:
658 //
659 for (Real i = 80; i < 100; ++i) {
660 Real h = i / 100;
661 auto f = [&](Real x)->Real
662 {
663 if (x > -log(boost::math::tools::epsilon<Real>()))
664 return exp((h - 1) * x);
665 else
666 {
667 Real et = exp(-x);
668 // Need to deal with numeric instability causing et to be greater than PI/2:
669 return et >= boost::math::constants::half_pi<Real>() ? 0 : et * pow(1 / tan(et), h);
670 }
671 };
672 Q = integrator.integrate(f, -log(half_pi<Real>()), boost::math::tools::max_value<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
673 Q_expected = half_pi<Real>() / cos(h*half_pi<Real>());
674 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 5 * tol);
675 }
676
677 // CRC Section 5.5, integral 670:
678
679 auto f3 = [](Real x)->Real { return sqrt(log(1/x)); };
680 Q = integrator.integrate(f3, (Real) 0, (Real) 1, get_convergence_tolerance<Real>(), &error, &L1);
681 Q_expected = root_pi<Real>()/2;
682 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
683
684 }
685
686 template <class Real>
test_sf()687 void test_sf()
688 {
689 using std::sqrt;
690 // Test some special functions that we already know how to evaluate:
691 std::cout << "Testing special functions on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
692 Real tol = 10 * boost::math::tools::epsilon<Real>();
693 auto integrator = get_integrator<Real>();
694
695 // incomplete beta:
696 if (std::numeric_limits<Real>::digits10 < 37) // Otherwise too slow
697 {
698 Real a(100), b(15);
699 auto f = [&](Real x)->Real { return boost::math::ibeta_derivative(a, b, x); };
700 BOOST_CHECK_CLOSE_FRACTION(integrator.integrate(f, 0, Real(0.25)), boost::math::ibeta(100, 15, Real(0.25)), tol * 10);
701 // Check some really extreme versions:
702 a = 1000;
703 b = 500;
704 BOOST_CHECK_CLOSE_FRACTION(integrator.integrate(f, 0, 1), Real(1), tol * 15);
705 //
706 // This is as extreme as we can get in this domain: otherwise the function has all it's
707 // area so close to zero we never get in there no matter how many levels we go down:
708 //
709 a = Real(1) / 15;
710 b = 30;
711 BOOST_CHECK_CLOSE_FRACTION(integrator.integrate(f, 0, 1), Real(1), tol * 25);
712 }
713
714 Real x = 2, y = 3, z = 0.5, p = 0.25;
715 // Elliptic integral RC:
716 Real Q = integrator.integrate([&](const Real& t)->Real { 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;
717 BOOST_CHECK_CLOSE_FRACTION(Q, boost::math::ellint_rc(x, y), tol);
718 // Elliptic Integral RJ:
719 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);
720
721 z = 5.5;
722 if (std::numeric_limits<Real>::digits10 > 40)
723 tol *= 200;
724 else if (!std::numeric_limits<Real>::is_specialized)
725 tol *= 3;
726 // tgamma expressed as an integral:
727 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>()),
728 boost::math::tgamma(z), tol);
729 BOOST_CHECK_CLOSE_FRACTION(integrator.integrate([](const Real& t)->Real { 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>()),
730 boost::math::constants::root_pi<Real>(), tol);
731 }
732
733 template <class Real>
test_2_arg()734 void test_2_arg()
735 {
736 BOOST_MATH_STD_USING
737 std::cout << "Testing 2 argument functors on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
738 Real tol = 10 * boost::math::tools::epsilon<Real>();
739
740 auto integrator = get_integrator<Real>();
741
742 //
743 // There are a whole family of integrals of the general form
744 // x(1-x)^-N ; N < 1
745 // which have all the interesting behaviour near the 2 singularities
746 // and all converge, see: http://www.wolframalpha.com/input/?i=integrate+(x+*+(1-x))%5E-1%2FN+from+0+to+1
747 //
748 Real Q = integrator.integrate([&](const Real& t, const Real & tc)->Real
749 {
750 return tc < 0 ? 1 / sqrt(t * (1-t)) : 1 / sqrt(t * tc);
751 }, 0, 1);
752 BOOST_CHECK_CLOSE_FRACTION(Q, boost::math::constants::pi<Real>(), tol);
753 Q = integrator.integrate([&](const Real& t, const Real & tc)->Real
754 {
755 return tc < 0 ? 1 / boost::math::cbrt(t * (1-t)) : 1 / boost::math::cbrt(t * tc);
756 }, 0, 1);
757 BOOST_CHECK_CLOSE_FRACTION(Q, boost::math::pow<2>(boost::math::tgamma(Real(2) / 3)) / boost::math::tgamma(Real(4) / 3), tol * 3);
758 //
759 // We can do the same thing with ((1+x)(1-x))^-N ; N < 1
760 //
761 Q = integrator.integrate([&](const Real& t, const Real & tc)->Real
762 {
763 return t < 0 ? 1 / sqrt(-tc * (1-t)) : 1 / sqrt((t + 1) * tc);
764 }, -1, 1);
765 BOOST_CHECK_CLOSE_FRACTION(Q, boost::math::constants::pi<Real>(), tol);
766 Q = integrator.integrate([&](const Real& t, const Real & tc)->Real
767 {
768 return t < 0 ? 1 / sqrt(-tc * (1-t)) : 1 / sqrt((t + 1) * tc);
769 });
770 BOOST_CHECK_CLOSE_FRACTION(Q, boost::math::constants::pi<Real>(), tol);
771 Q = integrator.integrate([&](const Real& t, const Real & tc)->Real
772 {
773 return t < 0 ? 1 / boost::math::cbrt(-tc * (1-t)) : 1 / boost::math::cbrt((t + 1) * tc);
774 }, -1, 1);
775 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);
776 Q = integrator.integrate([&](const Real& t, const Real & tc)->Real
777 {
778 return t < 0 ? 1 / boost::math::cbrt(-tc * (1-t)) : 1 / boost::math::cbrt((t + 1) * tc);
779 });
780 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);
781 //
782 // These are taken from above, and do not get to full precision via the single arg functor:
783 //
784 auto f7 = [](const Real& t, const Real& tc) { return t < 1 ? sqrt(tan(t)) : sqrt(1/tan(tc)); };
785 Real error, L1;
786 Q = integrator.integrate(f7, (Real)0, (Real)half_pi<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
787 Real Q_expected = pi<Real>() / root_two<Real>();
788 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
789 BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
790
791 auto f8 = [](const Real& t, const Real& tc) { return t < 1 ? log(cos(t)) : log(sin(tc)); };
792 Q = integrator.integrate(f8, (Real)0, half_pi<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
793 Q_expected = -pi<Real>()*ln_two<Real>()*half<Real>();
794 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
795 BOOST_CHECK_CLOSE_FRACTION(L1, -Q_expected, tol);
796 }
797
798 template <class Complex>
test_complex()799 void test_complex()
800 {
801 typedef typename Complex::value_type value_type;
802 //
803 // Integral version of the confluent hypergeometric function:
804 // https://dlmf.nist.gov/13.4#E1
805 //
806 value_type tol = 10 * boost::math::tools::epsilon<value_type>();
807 Complex a(2, 3), b(3, 4), z(0.5, -2);
808
809 auto f = [&](value_type t)
810 {
811 return exp(z * t) * pow(t, a - value_type(1)) * pow(value_type(1) - t, b - a - value_type(1));
812 };
813
814 auto integrator = get_integrator<value_type>();
815 auto Q = integrator.integrate(f, value_type(0), value_type(1), get_convergence_tolerance<value_type>());
816 //
817 // Expected result computed from http://www.wolframalpha.com/input/?i=1F1%5B(2%2B3i),+(3%2B4i);+(0.5-2i)%5D+*+gamma(2%2B3i)+*+gamma(1%2Bi)+%2F+gamma(3%2B4i)
818 //
819 Complex expected(boost::lexical_cast<value_type>("-0.2911081612888249710582867318081776512805281815037891183828405999609246645054069649838607112484426042883371996"),
820 boost::lexical_cast<value_type>("0.4507983563969959578849120188097153649211346293694903758252662015991543519595834937475296809912196906074655385"));
821
822 value_type error = abs(expected - Q);
823 BOOST_CHECK_LE(error, tol);
824
825 //
826 // Sin Integral https://dlmf.nist.gov/6.2#E9
827 //
828 auto f2 = [z](value_type t)
829 {
830 return -exp(-z * cos(t)) * cos(z * sin(t));
831 };
832 Q = integrator.integrate(f2, value_type(0), boost::math::constants::half_pi<value_type>(), get_convergence_tolerance<value_type>());
833
834 expected = Complex(boost::lexical_cast<value_type>("0.8893822921008980697856313681734926564752476188106405688951257340480164694708337246829840859633322683740376134733"),
835 -boost::lexical_cast<value_type>("2.381380802906111364088958767973164614925936185337231718483495612539455538280372745733208000514737758457795502168"));
836 expected -= boost::math::constants::half_pi<value_type>();
837
838 error = abs(expected - Q);
839 BOOST_CHECK_LE(error, tol);
840 }
841
842
BOOST_AUTO_TEST_CASE(tanh_sinh_quadrature_test)843 BOOST_AUTO_TEST_CASE(tanh_sinh_quadrature_test)
844 {
845 #ifdef GENERATE_CONSTANTS
846 //
847 // Generate pre-computed coefficients:
848 std::cout << std::setprecision(35);
849 print_levels(generate_constants<cpp_bin_float_100>(10, 8), "L");
850
851 #else
852
853 #ifdef TEST1
854
855 test_right_limit_infinite<float>();
856 test_left_limit_infinite<float>();
857 test_linear<float>();
858 test_quadratic<float>();
859 test_singular<float>();
860 test_ca<float>();
861 test_three_quadrature_schemes_examples<float>();
862 test_horrible<float>();
863 test_integration_over_real_line<float>();
864 test_nr_examples<float>();
865 #endif
866 #ifdef TEST1A
867 test_early_termination<float>();
868 test_2_arg<float>();
869 #endif
870 #ifdef TEST1B
871 test_crc<float>();
872 #endif
873 #ifdef TEST2
874 test_right_limit_infinite<double>();
875 test_left_limit_infinite<double>();
876 test_linear<double>();
877 test_quadratic<double>();
878 test_singular<double>();
879 test_ca<double>();
880 test_three_quadrature_schemes_examples<double>();
881 test_horrible<double>();
882 test_integration_over_real_line<double>();
883 test_nr_examples<double>();
884 test_early_termination<double>();
885 test_sf<double>();
886 test_2_arg<double>();
887 #endif
888 #ifdef TEST2A
889 test_crc<double>();
890 #endif
891
892 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
893
894 #ifdef TEST3
895 test_right_limit_infinite<long double>();
896 test_left_limit_infinite<long double>();
897 test_linear<long double>();
898 test_quadratic<long double>();
899 test_singular<long double>();
900 test_ca<long double>();
901 test_three_quadrature_schemes_examples<long double>();
902 test_horrible<long double>();
903 test_integration_over_real_line<long double>();
904 test_nr_examples<long double>();
905 test_early_termination<long double>();
906 test_sf<long double>();
907 test_2_arg<long double>();
908 #endif
909 #ifdef TEST3A
910 test_crc<long double>();
911
912 #endif
913 #endif
914
915 #ifdef TEST4
916 test_right_limit_infinite<cpp_bin_float_quad>();
917 test_left_limit_infinite<cpp_bin_float_quad>();
918 test_linear<cpp_bin_float_quad>();
919 test_quadratic<cpp_bin_float_quad>();
920 test_singular<cpp_bin_float_quad>();
921 test_ca<cpp_bin_float_quad>();
922 test_three_quadrature_schemes_examples<cpp_bin_float_quad>();
923 test_horrible<cpp_bin_float_quad>();
924 test_nr_examples<cpp_bin_float_quad>();
925 test_early_termination<cpp_bin_float_quad>();
926 test_crc<cpp_bin_float_quad>();
927 test_sf<cpp_bin_float_quad>();
928 test_2_arg<cpp_bin_float_quad>();
929
930 #endif
931 #ifdef TEST5
932
933 test_sf<cpp_bin_float_50>();
934 test_sf<cpp_bin_float_100>();
935 test_sf<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<150> > >();
936
937 #endif
938 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
939 #ifdef TEST6
940
941 test_right_limit_infinite<boost::math::concepts::real_concept>();
942 test_left_limit_infinite<boost::math::concepts::real_concept>();
943 test_linear<boost::math::concepts::real_concept>();
944 test_quadratic<boost::math::concepts::real_concept>();
945 test_singular<boost::math::concepts::real_concept>();
946 test_ca<boost::math::concepts::real_concept>();
947 test_three_quadrature_schemes_examples<boost::math::concepts::real_concept>();
948 test_horrible<boost::math::concepts::real_concept>();
949 test_integration_over_real_line<boost::math::concepts::real_concept>();
950 test_nr_examples<boost::math::concepts::real_concept>();
951 test_early_termination<boost::math::concepts::real_concept>();
952 test_sf<boost::math::concepts::real_concept>();
953 test_2_arg<boost::math::concepts::real_concept>();
954 #endif
955 #ifdef TEST6A
956 test_crc<boost::math::concepts::real_concept>();
957
958 #endif
959 #endif
960 #ifdef TEST7
961 test_sf<cpp_dec_float_50>();
962 #endif
963 #if defined(TEST8) && defined(BOOST_HAS_FLOAT128)
964
965 test_right_limit_infinite<boost::multiprecision::float128>();
966 test_left_limit_infinite<boost::multiprecision::float128>();
967 test_linear<boost::multiprecision::float128>();
968 test_quadratic<boost::multiprecision::float128>();
969 test_singular<boost::multiprecision::float128>();
970 test_ca<boost::multiprecision::float128>();
971 test_three_quadrature_schemes_examples<boost::multiprecision::float128>();
972 test_horrible<boost::multiprecision::float128>();
973 test_integration_over_real_line<boost::multiprecision::float128>();
974 test_nr_examples<boost::multiprecision::float128>();
975 test_early_termination<boost::multiprecision::float128>();
976 test_crc<boost::multiprecision::float128>();
977 test_sf<boost::multiprecision::float128>();
978 test_2_arg<boost::multiprecision::float128>();
979
980 #endif
981 #ifdef TEST9
982 test_complex<std::complex<double> >();
983 test_complex<std::complex<float> >();
984 #endif
985
986
987 #endif
988 }
989
990 #else
991
main()992 int main() { return 0; }
993
994 #endif
995