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 adaptive_gauss_kronrod_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/gauss_kronrod.hpp>
18 #include <boost/math/special_functions/sinc.hpp>
19 #include <boost/multiprecision/cpp_bin_float.hpp>
20
21 #if !defined(TEST1) && !defined(TEST1A) && !defined(TEST2) && !defined(TEST3)
22 # define TEST1
23 # define TEST1A
24 # define TEST2
25 # define TEST3
26 #endif
27
28 #ifdef _MSC_VER
29 #pragma warning(disable:4127) // Conditional expression is constant
30 #endif
31
32 using std::expm1;
33 using std::atan;
34 using std::tan;
35 using std::log;
36 using std::log1p;
37 using std::asinh;
38 using std::atanh;
39 using std::sqrt;
40 using std::isnormal;
41 using std::abs;
42 using std::sinh;
43 using std::tanh;
44 using std::cosh;
45 using std::pow;
46 using std::exp;
47 using std::sin;
48 using std::cos;
49 using std::string;
50 using boost::math::quadrature::gauss_kronrod;
51 using boost::math::constants::pi;
52 using boost::math::constants::half_pi;
53 using boost::math::constants::two_div_pi;
54 using boost::math::constants::two_pi;
55 using boost::math::constants::half;
56 using boost::math::constants::third;
57 using boost::math::constants::half;
58 using boost::math::constants::third;
59 using boost::math::constants::catalan;
60 using boost::math::constants::ln_two;
61 using boost::math::constants::root_two;
62 using boost::math::constants::root_two_pi;
63 using boost::math::constants::root_pi;
64 using boost::multiprecision::cpp_bin_float_quad;
65
66 template <class Real>
get_termination_condition()67 Real get_termination_condition()
68 {
69 return boost::math::tools::epsilon<Real>() * 1000;
70 }
71
72
73 template<class Real, unsigned Points>
test_linear()74 void test_linear()
75 {
76 std::cout << "Testing linear functions are integrated properly by gauss_kronrod on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
77 Real tol = boost::math::tools::epsilon<Real>() * 10;
78 Real error;
79 auto f = [](const Real& x)
80 {
81 return 5*x + 7;
82 };
83 Real L1;
84 Real Q = gauss_kronrod<Real, Points>::integrate(f, (Real) 0, (Real) 1, 15, get_termination_condition<Real>(), &error, &L1);
85 BOOST_CHECK_CLOSE_FRACTION(Q, 9.5, tol);
86 BOOST_CHECK_CLOSE_FRACTION(L1, 9.5, tol);
87 BOOST_CHECK_LE(fabs(error / Q), get_termination_condition<Real>());
88 BOOST_CHECK_GE(fabs(error), fabs(Q - 9.5));
89 }
90
91 template<class Real, unsigned Points>
test_quadratic()92 void test_quadratic()
93 {
94 std::cout << "Testing quadratic functions are integrated properly by gauss-kronrod on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
95 Real tol = boost::math::tools::epsilon<Real>() * 10;
96 Real error;
97
98 auto f = [](const Real& x) { return 5*x*x + 7*x + 12; };
99 Real L1;
100 Real Q = gauss_kronrod<Real, Points>::integrate(f, 0, 1, 15, get_termination_condition<Real>(), &error, &L1);
101 BOOST_CHECK_CLOSE_FRACTION(Q, (Real) 17 + half<Real>()*third<Real>(), tol);
102 BOOST_CHECK_CLOSE_FRACTION(L1, (Real) 17 + half<Real>()*third<Real>(), tol);
103 BOOST_CHECK_LE(fabs(error / Q), get_termination_condition<Real>());
104 BOOST_CHECK_GE(fabs(error), fabs(Q - ((Real)17 + half<Real>()*third<Real>())));
105 }
106
107 // Examples taken from
108 //http://crd-legacy.lbl.gov/~dhbailey/dhbpapers/quadrature.pdf
109 template<class Real, unsigned Points>
test_ca()110 void test_ca()
111 {
112 std::cout << "Testing integration of C(a) on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
113 Real tol = boost::math::tools::epsilon<Real>() * 10;
114 Real L1;
115 Real error;
116
117 auto f1 = [](const Real& x) { return atan(x)/(x*(x*x + 1)) ; };
118 Real Q = gauss_kronrod<Real, Points>::integrate(f1, 0, 1, 15, get_termination_condition<Real>(), &error, &L1);
119 Real Q_expected = pi<Real>()*ln_two<Real>()/8 + catalan<Real>()*half<Real>();
120 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
121 BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
122 BOOST_CHECK_LE(fabs(error / Q), get_termination_condition<Real>());
123 BOOST_CHECK_GE(fabs(error), fabs(Q - Q_expected));
124
125 auto f2 = [](Real x)->Real { Real t0 = x*x + 1; Real t1 = sqrt(t0); return atan(t1)/(t0*t1); };
126 Q = gauss_kronrod<Real, Points>::integrate(f2, 0 , 1, 15, get_termination_condition<Real>(), &error, &L1);
127 Q_expected = pi<Real>()/4 - pi<Real>()/root_two<Real>() + 3*atan(root_two<Real>())/root_two<Real>();
128 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
129 BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
130 BOOST_CHECK_LE(fabs(error / Q), get_termination_condition<Real>());
131 BOOST_CHECK_GE(fabs(error), fabs(Q - Q_expected));
132
133 auto f5 = [](Real t)->Real { return t*t*log(t)/((t*t - 1)*(t*t*t*t + 1)); };
134 Q = gauss_kronrod<Real, Points>::integrate(f5, 0, 1, 25);
135 Q_expected = pi<Real>()*pi<Real>()*(2 - root_two<Real>())/32;
136 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 100 * tol);
137 }
138
139 template<class Real, unsigned Points>
test_three_quadrature_schemes_examples()140 void test_three_quadrature_schemes_examples()
141 {
142 std::cout << "Testing integral in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
143 Real tol = boost::math::tools::epsilon<Real>() * 10;
144 Real Q;
145 Real Q_expected;
146
147 // Example 1:
148 auto f1 = [](const Real& t) { return t*boost::math::log1p(t); };
149 Q = gauss_kronrod<Real, Points>::integrate(f1, 0 , 1);
150 Q_expected = half<Real>()*half<Real>();
151 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
152
153
154 // Example 2:
155 auto f2 = [](const Real& t) { return t*t*atan(t); };
156 Q = gauss_kronrod<Real, Points>::integrate(f2, 0, 1);
157 Q_expected = (pi<Real>() -2 + 2*ln_two<Real>())/12;
158 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 2 * tol);
159
160 // Example 3:
161 auto f3 = [](const Real& t) { return exp(t)*cos(t); };
162 Q = gauss_kronrod<Real, Points>::integrate(f3, 0, half_pi<Real>());
163 Q_expected = boost::math::expm1(half_pi<Real>())*half<Real>();
164 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
165
166 // Example 4:
167 auto f4 = [](Real x)->Real { Real t0 = sqrt(x*x + 2); return atan(t0)/(t0*(x*x+1)); };
168 Q = gauss_kronrod<Real, Points>::integrate(f4, 0, 1);
169 Q_expected = 5*pi<Real>()*pi<Real>()/96;
170 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
171 }
172
173
174 template<class Real, unsigned Points>
test_integration_over_real_line()175 void test_integration_over_real_line()
176 {
177 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";
178 Real tol = boost::math::tools::epsilon<Real>() * 10;
179 Real Q;
180 Real Q_expected;
181 Real L1;
182 Real error;
183
184 auto f1 = [](const Real& t) { return 1/(1+t*t);};
185 Q = gauss_kronrod<Real, Points>::integrate(f1, -boost::math::tools::max_value<Real>(), boost::math::tools::max_value<Real>(), 15, get_termination_condition<Real>(), &error, &L1);
186 Q_expected = pi<Real>();
187 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
188 BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
189 BOOST_CHECK_LE(fabs(error / Q), get_termination_condition<Real>());
190
191 auto f4 = [](const Real& t) { return 1/cosh(t);};
192 Q = gauss_kronrod<Real, Points>::integrate(f4, -boost::math::tools::max_value<Real>(), boost::math::tools::max_value<Real>(), 15, get_termination_condition<Real>(), &error, &L1);
193 Q_expected = pi<Real>();
194 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
195 BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
196 BOOST_CHECK_LE(fabs(error / Q), get_termination_condition<Real>());
197
198 }
199
200 template<class Real, unsigned Points>
test_right_limit_infinite()201 void test_right_limit_infinite()
202 {
203 std::cout << "Testing right limit infinite for Gauss-Kronrod in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
204 Real tol = boost::math::tools::epsilon<Real>() * 10;
205 Real Q;
206 Real Q_expected;
207 Real L1;
208 Real error;
209
210 // Example 11:
211 auto f1 = [](const Real& t) { return 1/(1+t*t);};
212 Q = gauss_kronrod<Real, Points>::integrate(f1, 0, boost::math::tools::max_value<Real>(), 15, get_termination_condition<Real>(), &error, &L1);
213 Q_expected = half_pi<Real>();
214 BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
215 BOOST_CHECK_LE(fabs(error / Q), get_termination_condition<Real>());
216
217 auto f4 = [](const Real& t) { return 1/(1+t*t); };
218 Q = gauss_kronrod<Real, Points>::integrate(f4, 1, boost::math::tools::max_value<Real>(), 15, get_termination_condition<Real>(), &error, &L1);
219 Q_expected = pi<Real>()/4;
220 BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
221 BOOST_CHECK_LE(fabs(error / Q), get_termination_condition<Real>());
222 }
223
224 template<class Real, unsigned Points>
test_left_limit_infinite()225 void test_left_limit_infinite()
226 {
227 std::cout << "Testing left limit infinite for Gauss-Kronrod in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
228 Real tol = boost::math::tools::epsilon<Real>() * 10;
229 Real Q;
230 Real Q_expected;
231
232 // Example 11:
233 auto f1 = [](const Real& t) { return 1/(1+t*t);};
234 Q = gauss_kronrod<Real, Points>::integrate(f1, -boost::math::tools::max_value<Real>(), 0);
235 Q_expected = half_pi<Real>();
236 BOOST_CHECK_CLOSE(Q, Q_expected, 300*tol);
237 }
238
BOOST_AUTO_TEST_CASE(gauss_quadrature_test)239 BOOST_AUTO_TEST_CASE(gauss_quadrature_test)
240 {
241 #ifdef TEST1
242 std::cout << "Testing with 15 point Gauss-Kronrod rule:\n";
243 test_linear<double, 15>();
244 test_quadratic<double, 15>();
245 test_ca<double, 15>();
246 test_three_quadrature_schemes_examples<double, 15>();
247 test_integration_over_real_line<double, 15>();
248 test_right_limit_infinite<double, 15>();
249 test_left_limit_infinite<double, 15>();
250
251 // test one case where we do not have pre-computed constants:
252 std::cout << "Testing with 17 point Gauss-Kronrod rule:\n";
253 test_linear<double, 17>();
254 test_quadratic<double, 17>();
255 test_ca<double, 17>();
256 test_three_quadrature_schemes_examples<double, 17>();
257 test_integration_over_real_line<double, 17>();
258 test_right_limit_infinite<double, 17>();
259 test_left_limit_infinite<double, 17>();
260 #endif
261 #ifdef TEST1A
262 std::cout << "Testing with 21 point Gauss-Kronrod rule:\n";
263 test_linear<cpp_bin_float_quad, 21>();
264 test_quadratic<cpp_bin_float_quad, 21>();
265 test_ca<cpp_bin_float_quad, 21>();
266 test_three_quadrature_schemes_examples<cpp_bin_float_quad, 21>();
267 test_integration_over_real_line<cpp_bin_float_quad, 21>();
268 test_right_limit_infinite<cpp_bin_float_quad, 21>();
269 test_left_limit_infinite<cpp_bin_float_quad, 21>();
270
271 std::cout << "Testing with 31 point Gauss-Kronrod rule:\n";
272 test_linear<cpp_bin_float_quad, 31>();
273 test_quadratic<cpp_bin_float_quad, 31>();
274 test_ca<cpp_bin_float_quad, 31>();
275 test_three_quadrature_schemes_examples<cpp_bin_float_quad, 31>();
276 test_integration_over_real_line<cpp_bin_float_quad, 31>();
277 test_right_limit_infinite<cpp_bin_float_quad, 31>();
278 test_left_limit_infinite<cpp_bin_float_quad, 31>();
279 #endif
280 #ifdef TEST2
281 std::cout << "Testing with 41 point Gauss-Kronrod rule:\n";
282 test_linear<cpp_bin_float_quad, 41>();
283 test_quadratic<cpp_bin_float_quad, 41>();
284 test_ca<cpp_bin_float_quad, 41>();
285 test_three_quadrature_schemes_examples<cpp_bin_float_quad, 41>();
286 test_integration_over_real_line<cpp_bin_float_quad, 41>();
287 test_right_limit_infinite<cpp_bin_float_quad, 41>();
288 test_left_limit_infinite<cpp_bin_float_quad, 41>();
289
290 std::cout << "Testing with 51 point Gauss-Kronrod rule:\n";
291 test_linear<cpp_bin_float_quad, 51>();
292 test_quadratic<cpp_bin_float_quad, 51>();
293 test_ca<cpp_bin_float_quad, 51>();
294 test_three_quadrature_schemes_examples<cpp_bin_float_quad, 51>();
295 test_integration_over_real_line<cpp_bin_float_quad, 51>();
296 test_right_limit_infinite<cpp_bin_float_quad, 51>();
297 test_left_limit_infinite<cpp_bin_float_quad, 51>();
298 #endif
299 #ifdef TEST3
300 std::cout << "Testing with 61 point Gauss-Kronrod rule:\n";
301 test_linear<cpp_bin_float_quad, 61>();
302 test_quadratic<cpp_bin_float_quad, 61>();
303 test_ca<cpp_bin_float_quad, 61>();
304 test_three_quadrature_schemes_examples<cpp_bin_float_quad, 61>();
305 test_integration_over_real_line<cpp_bin_float_quad, 61>();
306 test_right_limit_infinite<cpp_bin_float_quad, 61>();
307 test_left_limit_infinite<cpp_bin_float_quad, 61>();
308 #endif
309 }
310
311 #else
312
main()313 int main() { return 0; }
314
315 #endif
316