1 ///////////////////////////////////////////////////////////////
2 // Copyright 2011 John Maddock. Distributed under the Boost
3 // Software License, Version 1.0. (See accompanying file
4 // LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
5
6 #define BOOST_MATH_MAX_ROOT_ITERATION_POLICY 750
7 #define BOOST_MATH_PROMOTE_DOUBLE_POLICY false
8
9 #if !defined(TEST_MPFR) && !defined(TEST_MPREAL) && !defined(TEST_MPF) && !defined(TEST_MPREAL) && !defined(TEST_CPP_DEC_FLOAT) && !defined(TEST_MPFR_CLASS) && !defined(TEST_FLOAT) && !defined(TEST_CPP_BIN_FLOAT)
10 #define TEST_MPFR
11 #define TEST_MPF
12 #define TEST_CPP_DEC_FLOAT
13 #define TEST_CPP_BIN_FLOAT
14 //# define TEST_MPFR_CLASS
15 //# define TEST_MPREAL
16 #define TEST_FLOAT
17 #endif
18
19 #if defined(TEST_MPFR) && !defined(TEST_MPFR_CLASS)
20 #if defined(__has_include)
21 #if __has_include(<gmpfrxx.h>)
22 #define TEST_MPFR_CLASS
23 #endif
24 #endif
25 #endif
26
27 #ifdef TEST_FLOAT
28 #include "arithmetic_backend.hpp"
29 #endif
30 #ifdef TEST_MPFR_CLASS
31 #include <boost/math/bindings/mpfr.hpp>
32 #endif
33 #ifdef TEST_MPFR
34 #include <boost/multiprecision/mpfr.hpp>
35 #endif
36 #ifdef TEST_MPREAL
37 #include <boost/math/bindings/mpreal.hpp>
38 #endif
39 #ifdef TEST_MPF
40 #include <boost/multiprecision/gmp.hpp>
41 #endif
42 #ifdef TEST_CPP_DEC_FLOAT
43 #include <boost/multiprecision/cpp_dec_float.hpp>
44 #endif
45 #ifdef TEST_CPP_BIN_FLOAT
46 #include <boost/multiprecision/cpp_bin_float.hpp>
47 #endif
48 #include <boost/math/special_functions/bessel.hpp>
49 #include <boost/math/tools/rational.hpp>
50 #include <boost/math/distributions/non_central_t.hpp>
51 #include <libs/math/test/table_type.hpp>
52 #include <boost/chrono.hpp>
53 #include <boost/array.hpp>
54 #include <boost/thread.hpp>
55 #include <boost/atomic.hpp>
56
57 template <class Real>
58 Real test_bessel();
59
60 template <class Clock>
61 struct stopwatch
62 {
63 typedef typename Clock::duration duration;
stopwatchstopwatch64 stopwatch()
65 {
66 m_start = Clock::now();
67 }
elapsedstopwatch68 duration elapsed()
69 {
70 return Clock::now() - m_start;
71 }
resetstopwatch72 void reset()
73 {
74 m_start = Clock::now();
75 }
76
77 private:
78 typename Clock::time_point m_start;
79 };
80
81 template <class Real>
test_bessel()82 Real test_bessel()
83 {
84 try
85 {
86 #define T double
87 #define SC_(x) x
88 #include "libs/math/test/bessel_i_int_data.ipp"
89 #include "libs/math/test/bessel_i_data.ipp"
90
91 Real r;
92
93 for (unsigned i = 0; i < bessel_i_int_data.size(); ++i)
94 {
95 r += boost::math::cyl_bessel_i(Real(bessel_i_int_data[i][0]), Real(bessel_i_int_data[i][1]));
96 }
97 for (unsigned i = 0; i < bessel_i_data.size(); ++i)
98 {
99 r += boost::math::cyl_bessel_i(Real(bessel_i_data[i][0]), Real(bessel_i_data[i][1]));
100 }
101
102 #include "libs/math/test/bessel_j_int_data.ipp"
103 for (unsigned i = 0; i < bessel_j_int_data.size(); ++i)
104 {
105 r += boost::math::cyl_bessel_j(Real(bessel_j_int_data[i][0]), Real(bessel_j_int_data[i][1]));
106 }
107
108 #include "libs/math/test/bessel_j_data.ipp"
109 for (unsigned i = 0; i < bessel_j_data.size(); ++i)
110 {
111 r += boost::math::cyl_bessel_j(Real(bessel_j_data[i][0]), Real(bessel_j_data[i][1]));
112 }
113
114 #include "libs/math/test/bessel_j_large_data.ipp"
115 for (unsigned i = 0; i < bessel_j_large_data.size(); ++i)
116 {
117 r += boost::math::cyl_bessel_j(Real(bessel_j_large_data[i][0]), Real(bessel_j_large_data[i][1]));
118 }
119
120 #include "libs/math/test/sph_bessel_data.ipp"
121 for (unsigned i = 0; i < sph_bessel_data.size(); ++i)
122 {
123 r += boost::math::sph_bessel(static_cast<unsigned>(sph_bessel_data[i][0]), Real(sph_bessel_data[i][1]));
124 }
125
126 return r;
127 }
128 catch (const std::exception& e)
129 {
130 std::cout << e.what() << std::endl;
131 }
132 return 0;
133 }
134
135 template <class Real>
test_polynomial()136 Real test_polynomial()
137 {
138 static const unsigned t[] = {
139 2, 3, 4, 5, 6, 7, 8};
140 Real result = 0;
141 for (Real k = 2; k < 1000; ++k)
142 result += boost::math::tools::evaluate_polynomial(t, k);
143
144 return result;
145 }
146
147 template <class Real>
test_nct()148 Real test_nct()
149 {
150 #define T double
151 #include "libs/math/test/nct.ipp"
152
153 Real result = 0;
154 for (unsigned i = 0; i < nct.size(); ++i)
155 {
156 try
157 {
158 result += quantile(boost::math::non_central_t_distribution<Real>(nct[i][0], nct[i][1]), nct[i][3]);
159 result += cdf(boost::math::non_central_t_distribution<Real>(nct[i][0], nct[i][1]), nct[i][2]);
160 }
161 catch (const std::exception&)
162 {}
163 }
164 return result;
165 }
166
167 extern boost::atomic<unsigned> allocation_count;
168 extern std::map<std::string, std::map<std::string, std::pair<double, unsigned> > > result_table;
169
170 template <class Real>
basic_allocation_test(const char * name,Real x)171 void basic_allocation_test(const char* name, Real x)
172 {
173 static const unsigned a[] = {2, 3, 4, 5, 6, 7, 8};
174 allocation_count = 0;
175 Real result = (((((a[6] * x + a[5]) * x + a[4]) * x + a[3]) * x + a[2]) * x + a[1]) * x + a[0];
176 std::cout << "Allocation count for type " << name << " = " << allocation_count << std::endl;
177 }
178
179 template <class Real>
poly_allocation_test(const char * name,Real x)180 void poly_allocation_test(const char* name, Real x)
181 {
182 static const unsigned a[] = {2, 3, 4, 5, 6, 7, 8};
183 allocation_count = 0;
184 Real result = boost::math::tools::evaluate_polynomial(a, x);
185 std::cout << "Allocation count for type " << name << " = " << allocation_count << std::endl;
186 }
187
188 template <class Real>
time_proc(const char * tablename,const char * name,Real (* proc)(),unsigned threads=1)189 void time_proc(const char* tablename, const char* name, Real (*proc)(), unsigned threads = 1)
190 {
191 try
192 {
193 static Real total = 0;
194 allocation_count = 0;
195 boost::chrono::duration<double> time;
196 stopwatch<boost::chrono::high_resolution_clock> c;
197 total += proc();
198 time = c.elapsed();
199 std::cout << "Time for " << name << " = " << time << std::endl;
200 std::cout << "Total allocations for " << name << " = " << allocation_count << std::endl;
201
202 result_table[tablename][name] = std::make_pair(time.count(), (unsigned)allocation_count);
203
204 if (threads > 1)
205 {
206 c.reset();
207 boost::thread_group g;
208 for (unsigned i = 0; i < threads; ++i)
209 g.create_thread(proc);
210 g.join_all();
211 time = c.elapsed();
212 std::cout << "Time for " << name << " (" << (threads) << " theads) = " << time << std::endl;
213 std::cout << "Total allocations for " << name << " = " << allocation_count << std::endl;
214
215 std::ostringstream ss;
216 ss << name << " (" << threads << " concurrent threads)";
217 result_table[tablename][ss.str()] = std::make_pair(time.count(), (unsigned)allocation_count);
218 }
219 }
220 catch (const std::exception& e)
221 {
222 std::cout << e.what() << std::endl;
223 }
224 }
225
226 using namespace boost::multiprecision;
227
228 void basic_tests_1();
229 void basic_tests_2();
230 void basic_tests_3();
231 void basic_tests_4();
232 void basic_tests_5();
233 void basic_tests_6();
234 void basic_tests_7();
235 void basic_tests_8();
236 void basic_tests_9();
237 void bessel_tests();
238 void poly_tests();
239 void nct_tests();
240