1 // Copyright John Maddock 2015.
2 // Use, modification and distribution are subject to the
3 // Boost Software License, Version 1.0. (See accompanying file
4 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
5
6 #ifdef _MSC_VER
7 # pragma warning (disable : 4224)
8 #endif
9
10 #include <boost/array.hpp>
11 #include <boost/lexical_cast.hpp>
12 #include "../../test/table_type.hpp"
13 #include "table_helper.hpp"
14 #include "performance.hpp"
15 #include <iostream>
16
17 #define evaluate_polynomial_c_imp evaluate_polynomial_c_imp_1
18 #undef BOOST_MATH_TOOLS_POLY_EVAL_20_HPP
19 #include <boost/math/tools/detail/polynomial_horner1_20.hpp>
20 #undef evaluate_polynomial_c_imp
21 #undef BOOST_MATH_TOOLS_POLY_EVAL_20_HPP
22 #define evaluate_polynomial_c_imp evaluate_polynomial_c_imp_2
23 #include <boost/math/tools/detail/polynomial_horner2_20.hpp>
24 #undef evaluate_polynomial_c_imp
25 #undef BOOST_MATH_TOOLS_POLY_EVAL_20_HPP
26 #define evaluate_polynomial_c_imp evaluate_polynomial_c_imp_3
27 #include <boost/math/tools/detail/polynomial_horner3_20.hpp>
28 #undef evaluate_polynomial_c_imp
29
30 #undef BOOST_MATH_TOOLS_POLY_RAT_20_HPP
31 #define evaluate_rational_c_imp evaluate_rational_c_imp_1
32 #include <boost/math/tools/detail/rational_horner1_20.hpp>
33 #undef evaluate_rational_c_imp
34 #undef BOOST_MATH_TOOLS_POLY_RAT_20_HPP
35 #define evaluate_rational_c_imp evaluate_rational_c_imp_2
36 #include <boost/math/tools/detail/rational_horner2_20.hpp>
37 #undef evaluate_rational_c_imp
38 #undef BOOST_MATH_TOOLS_RAT_EVAL_20_HPP
39 #define evaluate_rational_c_imp evaluate_rational_c_imp_3
40 #include <boost/math/tools/detail/rational_horner3_20.hpp>
41 #undef evaluate_rational_c_imp
42 #undef BOOST_MATH_TOOLS_POLY_RAT_20_HPP
43
44 static const double num[21] = {
45 static_cast<double>(56906521.91347156388090791033559122686859L),
46 static_cast<double>(103794043.1163445451906271053616070238554L),
47 static_cast<double>(86363131.28813859145546927288977868422342L),
48 static_cast<double>(43338889.32467613834773723740590533316085L),
49 static_cast<double>(14605578.08768506808414169982791359218571L),
50 static_cast<double>(3481712.15498064590882071018964774556468L),
51 static_cast<double>(601859.6171681098786670226533699352302507L),
52 static_cast<double>(75999.29304014542649875303443598909137092L),
53 static_cast<double>(6955.999602515376140356310115515198987526L),
54 static_cast<double>(449.9445569063168119446858607650988409623L),
55 static_cast<double>(19.51992788247617482847860966235652136208L),
56 static_cast<double>(0.5098416655656676188125178644804694509993L),
57 static_cast<double>(0.006061842346248906525783753964555936883222L),
58 0.0
59 };
60 static const double denom[20] = {
61 static_cast<double>(0u),
62 static_cast<double>(39916800u),
63 static_cast<double>(120543840u),
64 static_cast<double>(150917976u),
65 static_cast<double>(105258076u),
66 static_cast<double>(45995730u),
67 static_cast<double>(13339535u),
68 static_cast<double>(2637558u),
69 static_cast<double>(357423u),
70 static_cast<double>(32670u),
71 static_cast<double>(1925u),
72 static_cast<double>(66u),
73 static_cast<double>(1u),
74 0.0
75 };
76 static const boost::uint32_t denom_int[20] = {
77 static_cast<boost::uint32_t>(0u),
78 static_cast<boost::uint32_t>(39916800u),
79 static_cast<boost::uint32_t>(120543840u),
80 static_cast<boost::uint32_t>(150917976u),
81 static_cast<boost::uint32_t>(105258076u),
82 static_cast<boost::uint32_t>(45995730u),
83 static_cast<boost::uint32_t>(13339535u),
84 static_cast<boost::uint32_t>(2637558u),
85 static_cast<boost::uint32_t>(357423u),
86 static_cast<boost::uint32_t>(32670u),
87 static_cast<boost::uint32_t>(1925u),
88 static_cast<boost::uint32_t>(66u),
89 static_cast<boost::uint32_t>(1u),
90 0
91 };
92
make_order_string(int n)93 std::string make_order_string(int n)
94 {
95 std::string result = boost::lexical_cast<std::string>(n);
96 if (result.size() < 2)
97 result.insert(result.begin(), ' ');
98 return result;
99 }
100
test_poly_1(const boost::integral_constant<int,1> &)101 void test_poly_1(const boost::integral_constant<int, 1>&)
102 {
103 }
104
105 template <int N>
test_poly_1(const boost::integral_constant<int,N> &)106 void test_poly_1(const boost::integral_constant<int, N>&)
107 {
108 test_poly_1(boost::integral_constant<int, N - 1>());
109
110 double time = exec_timed_test([](const std::vector<double>& v)
111 {
112 double result = 0;
113 for (unsigned i = 0; i < 10; ++i)
114 result += boost::math::tools::detail::evaluate_polynomial_c_imp_1(denom, v[0] + i, static_cast<boost::integral_constant<int, N>*>(0));
115 return result;
116 });
117 report_execution_time(time, std::string("Polynomial Method Comparison with ") + compiler_name() + std::string(" on ") + platform_name(), "Order " + make_order_string(N), "Method 1[br](Double Coefficients)");
118
119 time = exec_timed_test([](const std::vector<double>& v)
120 {
121 double result = 0;
122 for (unsigned i = 0; i < 10; ++i)
123 result += boost::math::tools::detail::evaluate_polynomial_c_imp_1(denom_int, v[0] + i, static_cast<boost::integral_constant<int, N>*>(0));
124 return result;
125 });
126 report_execution_time(time, std::string("Polynomial Method Comparison with ") + compiler_name() + std::string(" on ") + platform_name(), "Order " + make_order_string(N), "Method 1[br](Integer Coefficients)");
127 }
128
129
test_poly_2(const boost::integral_constant<int,1> &)130 void test_poly_2(const boost::integral_constant<int, 1>&)
131 {
132 }
133
134 template <int N>
test_poly_2(const boost::integral_constant<int,N> &)135 void test_poly_2(const boost::integral_constant<int, N>&)
136 {
137 test_poly_2(boost::integral_constant<int, N - 1>());
138
139 double time = exec_timed_test([](const std::vector<double>& v)
140 {
141 double result = 0;
142 for (unsigned i = 0; i < 10; ++i)
143 result += boost::math::tools::detail::evaluate_polynomial_c_imp_2(denom, v[0] + i, static_cast<boost::integral_constant<int, N>*>(0));
144 return result;
145 });
146 report_execution_time(time, std::string("Polynomial Method Comparison with ") + compiler_name() + std::string(" on ") + platform_name(), "Order " + make_order_string(N), "Method 2[br](Double Coefficients)");
147
148 time = exec_timed_test([](const std::vector<double>& v)
149 {
150 double result = 0;
151 for (unsigned i = 0; i < 10; ++i)
152 result += boost::math::tools::detail::evaluate_polynomial_c_imp_2(denom_int, v[0] + i, static_cast<boost::integral_constant<int, N>*>(0));
153 return result;
154 });
155 report_execution_time(time, std::string("Polynomial Method Comparison with ") + compiler_name() + std::string(" on ") + platform_name(), "Order " + make_order_string(N), "Method 2[br](Integer Coefficients)");
156 }
157
test_poly_3(const boost::integral_constant<int,1> &)158 void test_poly_3(const boost::integral_constant<int, 1>&)
159 {
160 }
161
162 template <int N>
test_poly_3(const boost::integral_constant<int,N> &)163 void test_poly_3(const boost::integral_constant<int, N>&)
164 {
165 test_poly_3(boost::integral_constant<int, N - 1>());
166
167 double time = exec_timed_test([](const std::vector<double>& v) { double result = 0;
168 for (unsigned i = 0; i < 10; ++i)
169 result += boost::math::tools::detail::evaluate_polynomial_c_imp_3(denom, v[0] + i, static_cast<boost::integral_constant<int, N>*>(0));
170 return result;
171 });
172 report_execution_time(time, std::string("Polynomial Method Comparison with ") + compiler_name() + std::string(" on ") + platform_name(), "Order " + make_order_string(N), "Method 3[br](Double Coefficients)");
173
174 time = exec_timed_test([](const std::vector<double>& v) { double result = 0;
175 for (unsigned i = 0; i < 10; ++i)
176 result += boost::math::tools::detail::evaluate_polynomial_c_imp_3(denom_int, v[0] + i, static_cast<boost::integral_constant<int, N>*>(0));
177 return result;
178 });
179 report_execution_time(time, std::string("Polynomial Method Comparison with ") + compiler_name() + std::string(" on ") + platform_name(), "Order " + make_order_string(N), "Method 3[br](Integer Coefficients)");
180 }
181
182 template <class T, class U>
evaluate_polynomial_0(const T * poly,U const & z,std::size_t count)183 U evaluate_polynomial_0(const T* poly, U const& z, std::size_t count)
184 {
185 U sum = static_cast<U>(poly[count - 1]);
186 for (int i = static_cast<int>(count) - 2; i >= 0; --i)
187 {
188 sum *= z;
189 sum += static_cast<U>(poly[i]);
190 }
191 return sum;
192 }
193
test_rat_1(const boost::integral_constant<int,1> &)194 void test_rat_1(const boost::integral_constant<int, 1>&)
195 {
196 }
197
198 template <int N>
test_rat_1(const boost::integral_constant<int,N> &)199 void test_rat_1(const boost::integral_constant<int, N>&)
200 {
201 test_rat_1(boost::integral_constant<int, N - 1>());
202
203 double time = exec_timed_test([](const std::vector<double>& v)
204 {
205 double result = 0;
206 for (unsigned i = 0; i < 10; ++i)
207 result += boost::math::tools::detail::evaluate_rational_c_imp_1(num, denom, v[0] + i, static_cast<boost::integral_constant<int, N>*>(0));
208 return result;
209 });
210 report_execution_time(time, std::string("Rational Method Comparison with ") + compiler_name() + std::string(" on ") + platform_name(), "Order " + make_order_string(N), "Method 1[br](Double Coefficients)");
211
212 time = exec_timed_test([](const std::vector<double>& v)
213 {
214 double result = 0;
215 for (unsigned i = 0; i < 10; ++i)
216 result += boost::math::tools::detail::evaluate_rational_c_imp_1(num, denom_int, v[0] + i, static_cast<boost::integral_constant<int, N>*>(0));
217 return result;
218 });
219 report_execution_time(time, std::string("Rational Method Comparison with ") + compiler_name() + std::string(" on ") + platform_name(), "Order " + make_order_string(N), "Method 1[br](Integer Coefficients)");
220 }
221
test_rat_2(const boost::integral_constant<int,1> &)222 void test_rat_2(const boost::integral_constant<int, 1>&)
223 {
224 }
225
226 template <int N>
test_rat_2(const boost::integral_constant<int,N> &)227 void test_rat_2(const boost::integral_constant<int, N>&)
228 {
229 test_rat_2(boost::integral_constant<int, N - 1>());
230
231 double time = exec_timed_test([](const std::vector<double>& v)
232 {
233 double result = 0;
234 for (unsigned i = 0; i < 10; ++i)
235 result += boost::math::tools::detail::evaluate_rational_c_imp_2(num, denom, v[0] + i, static_cast<boost::integral_constant<int, N>*>(0));
236 return result;
237 });
238 report_execution_time(time, std::string("Rational Method Comparison with ") + compiler_name() + std::string(" on ") + platform_name(), "Order " + make_order_string(N), "Method 2[br](Double Coefficients)");
239
240 time = exec_timed_test([](const std::vector<double>& v)
241 {
242 double result = 0;
243 for (unsigned i = 0; i < 10; ++i)
244 result += boost::math::tools::detail::evaluate_rational_c_imp_2(num, denom_int, v[0] + i, static_cast<boost::integral_constant<int, N>*>(0));
245 return result;
246 });
247 report_execution_time(time, std::string("Rational Method Comparison with ") + compiler_name() + std::string(" on ") + platform_name(), "Order " + make_order_string(N), "Method 2[br](Integer Coefficients)");
248 }
249
test_rat_3(const boost::integral_constant<int,1> &)250 void test_rat_3(const boost::integral_constant<int, 1>&)
251 {
252 }
253
254 template <int N>
test_rat_3(const boost::integral_constant<int,N> &)255 void test_rat_3(const boost::integral_constant<int, N>&)
256 {
257 test_rat_3(boost::integral_constant<int, N - 1>());
258
259 double time = exec_timed_test([](const std::vector<double>& v)
260 {
261 double result = 0;
262 for (unsigned i = 0; i < 10; ++i)
263 result += boost::math::tools::detail::evaluate_rational_c_imp_3(num, denom, v[0] + i, static_cast<boost::integral_constant<int, N>*>(0));
264 return result;
265 });
266 report_execution_time(time, std::string("Rational Method Comparison with ") + compiler_name() + std::string(" on ") + platform_name(), "Order " + make_order_string(N), "Method 3[br](Double Coefficients)");
267
268 time = exec_timed_test([](const std::vector<double>& v)
269 {
270 double result = 0;
271 for (unsigned i = 0; i < 10; ++i)
272 result += boost::math::tools::detail::evaluate_rational_c_imp_3(num, denom_int, v[0] + i, static_cast<boost::integral_constant<int, N>*>(0));
273 return result;
274 });
275 report_execution_time(time, std::string("Rational Method Comparison with ") + compiler_name() + std::string(" on ") + platform_name(), "Order " + make_order_string(N), "Method 3[br](Integer Coefficients)");
276 }
277
278 template <class T, class U, class V>
evaluate_rational_0(const T * num,const U * denom,const V & z_,std::size_t count)279 V evaluate_rational_0(const T* num, const U* denom, const V& z_, std::size_t count)
280 {
281 V z(z_);
282 V s1, s2;
283 if (z <= 1)
284 {
285 s1 = static_cast<V>(num[count - 1]);
286 s2 = static_cast<V>(denom[count - 1]);
287 for (int i = (int)count - 2; i >= 0; --i)
288 {
289 s1 *= z;
290 s2 *= z;
291 s1 += num[i];
292 s2 += denom[i];
293 }
294 }
295 else
296 {
297 z = 1 / z;
298 s1 = static_cast<V>(num[0]);
299 s2 = static_cast<V>(denom[0]);
300 for (unsigned i = 1; i < count; ++i)
301 {
302 s1 *= z;
303 s2 *= z;
304 s1 += num[i];
305 s2 += denom[i];
306 }
307 }
308 return s1 / s2;
309 }
310
311
main()312 int main()
313 {
314 double val = 0.001;
315 while (val < 1)
316 {
317 std::vector<double> v;
318 v.push_back(val);
319 data.push_back(v);
320 val *= 1.1;
321 }
322
323 for (unsigned i = 3; i <= 20; ++i)
324 {
325 double time = exec_timed_test([&](const std::vector<double>& v) {
326 double result = 0;
327 for (unsigned j = 0; j < 10; ++j)
328 result += evaluate_polynomial_0(denom, v[0] + j, i);
329 return result;
330 });
331 report_execution_time(time, std::string("Polynomial Method Comparison with ") + compiler_name() + std::string(" on ") + platform_name(), "Order " + make_order_string(i), "Method 0[br](Double Coefficients)");
332
333 time = exec_timed_test([&](const std::vector<double>& v) {
334 double result = 0;
335 for (unsigned j = 0; j < 10; ++j)
336 result += evaluate_polynomial_0(denom_int, v[0] + j, i);
337 return result;
338 });
339 report_execution_time(time, std::string("Polynomial Method Comparison with ") + compiler_name() + std::string(" on ") + platform_name(), "Order " + make_order_string(i), "Method 0[br](Integer Coefficients)");
340 }
341
342 test_poly_1(boost::integral_constant<int, 20>());
343 test_poly_2(boost::integral_constant<int, 20>());
344 test_poly_3(boost::integral_constant<int, 20>());
345
346 for (unsigned i = 3; i <= 20; ++i)
347 {
348 double time = exec_timed_test([&](const std::vector<double>& v) {
349 double result = 0;
350 for (unsigned j = 0; j < 10; ++j)
351 result += evaluate_rational_0(num, denom, v[0] + j, i);
352 return result;
353 });
354 report_execution_time(time, std::string("Rational Method Comparison with ") + compiler_name() + std::string(" on ") + platform_name(), "Order " + make_order_string(i), "Method 0[br](Double Coefficients)");
355
356 time = exec_timed_test([&](const std::vector<double>& v) {
357 double result = 0;
358 for (unsigned j = 0; j < 10; ++j)
359 result += evaluate_rational_0(num, denom_int, v[0] + j, i);
360 return result;
361 });
362 report_execution_time(time, std::string("Rational Method Comparison with ") + compiler_name() + std::string(" on ") + platform_name(), "Order " + make_order_string(i), "Method 0[br](Integer Coefficients)");
363 }
364
365 test_rat_1(boost::integral_constant<int, 20>());
366 test_rat_2(boost::integral_constant<int, 20>());
367 test_rat_3(boost::integral_constant<int, 20>());
368
369 return 0;
370 }
371
372