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 sinh_sinh_quadrature_test
8 #include <complex>
9 #include <boost/multiprecision/cpp_complex.hpp>
10 #include <boost/math/concepts/real_concept.hpp>
11 #include <boost/test/included/unit_test.hpp>
12 #include <boost/test/tools/floating_point_comparison.hpp>
13 #include <boost/math/quadrature/sinh_sinh.hpp>
14 #include <boost/math/special_functions/sinc.hpp>
15 #include <boost/multiprecision/cpp_bin_float.hpp>
16
17 #if !BOOST_WORKAROUND(BOOST_MSVC, < 1900)
18 // MSVC-12 has problems if we include 2 different multiprecision types in the same program,
19 // basically random things stop compiling, even though they work fine in isolation :(
20 #include <boost/multiprecision/cpp_dec_float.hpp>
21 #endif
22
23 using std::expm1;
24 using std::exp;
25 using std::sin;
26 using std::cos;
27 using std::atan;
28 using std::tan;
29 using std::log;
30 using std::log1p;
31 using std::asinh;
32 using std::atanh;
33 using std::sqrt;
34 using std::isnormal;
35 using std::abs;
36 using std::sinh;
37 using std::tanh;
38 using std::cosh;
39 using std::pow;
40 using std::string;
41 using boost::multiprecision::cpp_bin_float_quad;
42 using boost::math::quadrature::sinh_sinh;
43 using boost::math::constants::pi;
44 using boost::math::constants::pi_sqr;
45 using boost::math::constants::half_pi;
46 using boost::math::constants::two_div_pi;
47 using boost::math::constants::half;
48 using boost::math::constants::third;
49 using boost::math::constants::half;
50 using boost::math::constants::third;
51 using boost::math::constants::catalan;
52 using boost::math::constants::ln_two;
53 using boost::math::constants::root_two;
54 using boost::math::constants::root_two_pi;
55 using boost::math::constants::root_pi;
56 //
57 // Code for generating the coefficients:
58 //
59 template <class T>
print_levels(const T & v,const char * suffix)60 void print_levels(const T& v, const char* suffix)
61 {
62 std::cout << "{\n";
63 for (unsigned i = 0; i < v.size(); ++i)
64 {
65 std::cout << " { ";
66 for (unsigned j = 0; j < v[i].size(); ++j)
67 {
68 std::cout << v[i][j] << suffix << ", ";
69 }
70 std::cout << "},\n";
71 }
72 std::cout << " };\n";
73 }
74
75 template <class T>
print_levels(const std::pair<T,T> & p,const char * suffix="")76 void print_levels(const std::pair<T, T>& p, const char* suffix = "")
77 {
78 std::cout << " static const std::vector<std::vector<Real> > abscissa = ";
79 print_levels(p.first, suffix);
80 std::cout << " static const std::vector<std::vector<Real> > weights = ";
81 print_levels(p.second, suffix);
82 }
83
84 template <class Real, class TargetType>
generate_constants(unsigned max_rows)85 std::pair<std::vector<std::vector<Real>>, std::vector<std::vector<Real>> > generate_constants(unsigned max_rows)
86 {
87 using boost::math::constants::half_pi;
88 using boost::math::constants::two_div_pi;
89 using boost::math::constants::pi;
90 auto g = [](Real t) { return sinh(half_pi<Real>()*sinh(t)); };
91 auto w = [](Real t) { return cosh(t)*half_pi<Real>()*cosh(half_pi<Real>()*sinh(t)); };
92
93 std::vector<std::vector<Real>> abscissa, weights;
94
95 std::vector<Real> temp;
96
97 Real t_max = log(2 * two_div_pi<Real>()*log(2 * two_div_pi<Real>()*sqrt(boost::math::tools::max_value<TargetType>())));
98
99 std::cout << "m_t_max = " << t_max << ";\n";
100
101 Real h = 1;
102 for (Real t = 1; t < t_max; t += h)
103 {
104 temp.push_back(g(t));
105 }
106 abscissa.push_back(temp);
107 temp.clear();
108
109 for (Real t = 1; t < t_max; t += h)
110 {
111 temp.push_back(w(t * h));
112 }
113 weights.push_back(temp);
114 temp.clear();
115
116 for (unsigned row = 1; row < max_rows; ++row)
117 {
118 h /= 2;
119 for (Real t = h; t < t_max; t += 2 * h)
120 temp.push_back(g(t));
121 abscissa.push_back(temp);
122 temp.clear();
123 }
124 h = 1;
125 for (unsigned row = 1; row < max_rows; ++row)
126 {
127 h /= 2;
128 for (Real t = h; t < t_max; t += 2 * h)
129 temp.push_back(w(t));
130 weights.push_back(temp);
131 temp.clear();
132 }
133
134 return std::make_pair(abscissa, weights);
135 }
136
137 template<class Real>
test_nr_examples()138 void test_nr_examples()
139 {
140 std::cout << "Testing type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
141 Real integration_limit = sqrt(boost::math::tools::epsilon<Real>());
142 Real tol = 10 * boost::math::tools::epsilon<Real>();
143 std::cout << std::setprecision(std::numeric_limits<Real>::digits10);
144 Real Q;
145 Real Q_expected;
146 Real L1;
147 Real error;
148 sinh_sinh<Real> integrator(10);
149
150 auto f0 = [](Real)->Real { return (Real) 0; };
151 Q = integrator.integrate(f0, integration_limit, &error, &L1);
152 Q_expected = 0;
153 BOOST_CHECK_SMALL(Q, tol);
154 BOOST_CHECK_SMALL(L1, tol);
155
156 // In spite of the poles at \pm i, we still get a doubling of the correct digits at each level of refinement.
157 auto f1 = [](const Real& t)->Real { return 1/(1+t*t); };
158 Q = integrator.integrate(f1, integration_limit, &error, &L1);
159 Q_expected = pi<Real>();
160 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
161 BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
162 #if defined(BOOST_MSVC) && (BOOST_MSVC < 1900)
163 auto f2 = [](const Real& x)->Real { return fabs(x) > boost::math::tools::log_max_value<Real>() ? 0 : exp(-x*x); };
164 #else
165 auto f2 = [](const Real& x)->Real { return exp(-x*x); };
166 #endif
167 Q = integrator.integrate(f2, integration_limit, &error, &L1);
168 Q_expected = root_pi<Real>();
169 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
170 BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
171
172 auto f5 = [](const Real& t)->Real { return 1/cosh(t);};
173 Q = integrator.integrate(f5, integration_limit, &error, &L1);
174 Q_expected = pi<Real>();
175 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
176 BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
177
178 // This oscillatory integral has rapid convergence because the oscillations get swamped by the exponential growth of the denominator,
179 // none the less the error is slightly higher than for the other cases:
180 tol *= 10;
181 auto f8 = [](const Real& t)->Real { return cos(t)/cosh(t);};
182 Q = integrator.integrate(f8, integration_limit, &error, &L1);
183 Q_expected = pi<Real>()/cosh(half_pi<Real>());
184 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
185 // Try again with progressively fewer arguments:
186 Q = integrator.integrate(f8, integration_limit);
187 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
188 Q = integrator.integrate(f8);
189 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
190 }
191
192 // Test formulas for in the CRC Handbook of Mathematical functions, 32nd edition.
193 template<class Real>
test_crc()194 void test_crc()
195 {
196 std::cout << "Testing CRC formulas on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
197 Real integration_limit = sqrt(boost::math::tools::epsilon<Real>());
198 Real tol = 10 * boost::math::tools::epsilon<Real>();
199 std::cout << std::setprecision(std::numeric_limits<Real>::digits10);
200 Real Q;
201 Real Q_expected;
202 Real L1;
203 Real error;
204 sinh_sinh<Real> integrator(10);
205
206 // CRC Definite integral 698:
207 auto f0 = [](Real x)->Real {
208 using std::sinh;
209 if(x == 0) {
210 return (Real) 1;
211 }
212 return x/sinh(x);
213 };
214 Q = integrator.integrate(f0, integration_limit, &error, &L1);
215 Q_expected = pi_sqr<Real>()/2;
216 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
217 BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
218
219
220 // CRC Definite integral 695:
221 auto f1 = [](Real x)->Real {
222 using std::sin; using std::sinh;
223 if(x == 0) {
224 return (Real) 1;
225 }
226 return (Real) sin(x)/sinh(x);
227 };
228 Q = integrator.integrate(f1, integration_limit, &error, &L1);
229 Q_expected = pi<Real>()*tanh(half_pi<Real>());
230 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
231 }
232
233 template<class Complex>
test_dirichlet_eta()234 void test_dirichlet_eta()
235 {
236 typedef typename Complex::value_type Real;
237 std::cout << "Testing Dirichlet eta function on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
238 Real tol = 10 * boost::math::tools::epsilon<Real>();
239 Complex Q;
240 sinh_sinh<Real> integrator(10);
241
242 //https://en.wikipedia.org/wiki/Dirichlet_eta_function, integral representations:
243 Complex z = {1,1};
244 auto eta = [&z](Real t)->Complex {
245 using std::exp;
246 using std::pow;
247 using boost::math::constants::pi;
248 Complex i = {0,1};
249 Complex num = pow((Real)1/ (Real)2 + i*t, -z);
250 Real denom = exp(pi<Real>()*t) + exp(-pi<Real>()*t);
251 Complex res = num/denom;
252 return res;
253 };
254 Q = integrator.integrate(eta);
255 // N[DirichletEta[1 + I], 150]
256 Complex Q_expected = {boost::lexical_cast<Real>("0.726559775062463263201495728547241386311129502735725787103568290594808442332084045617744978600192784188182345866652233650512117834307254514480657408096"),
257 boost::lexical_cast<Real>("0.158095863901207324355426285544321998253687969756843115763682522207208309489794631247865357375538028170751576870244296106203144195376645765556607038775")};
258
259 BOOST_CHECK_CLOSE_FRACTION(Q.real(), Q_expected.real(), tol);
260 BOOST_CHECK_CLOSE_FRACTION(Q.imag(), Q_expected.imag(), tol);
261 }
262
263
BOOST_AUTO_TEST_CASE(sinh_sinh_quadrature_test)264 BOOST_AUTO_TEST_CASE(sinh_sinh_quadrature_test)
265 {
266 //
267 // Uncomment the following to print out the coefficients:
268 //
269 /*
270 std::cout << std::scientific << std::setprecision(8);
271 print_levels(generate_constants<cpp_bin_float_100, float>(8), "f");
272 std::cout << std::setprecision(18);
273 print_levels(generate_constants<cpp_bin_float_100, double>(8), "");
274 std::cout << std::setprecision(35);
275 print_levels(generate_constants<cpp_bin_float_100, cpp_bin_float_quad>(8), "L");
276 */
277 test_nr_examples<float>();
278 test_nr_examples<double>();
279 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
280 test_nr_examples<long double>();
281 #endif
282 test_nr_examples<cpp_bin_float_quad>();
283 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
284 test_nr_examples<boost::math::concepts::real_concept>();
285 #endif
286 #if !BOOST_WORKAROUND(BOOST_MSVC, < 1900)
287 test_nr_examples<boost::multiprecision::cpp_dec_float_50>();
288 #endif
289
290 test_crc<float>();
291 test_crc<double>();
292 test_dirichlet_eta<std::complex<double>>();
293 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
294 test_crc<long double>();
295 test_dirichlet_eta<std::complex<long double>>();
296 #endif
297 test_crc<cpp_bin_float_quad>();
298 test_dirichlet_eta<boost::multiprecision::cpp_complex_quad>();
299 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
300 test_crc<boost::math::concepts::real_concept>();
301 #endif
302 #if !BOOST_WORKAROUND(BOOST_MSVC, < 1900)
303 test_crc<boost::multiprecision::cpp_dec_float_50>();
304 #endif
305
306
307 }
308