1 // Copyright Jeremy Murphy 2016.
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/integer/common_factor_rt.hpp>
11 #include <boost/math/special_functions/prime.hpp>
12 #include <boost/multiprecision/cpp_int.hpp>
13 #include <boost/multiprecision/integer.hpp>
14 #include <boost/random.hpp>
15 #include <boost/array.hpp>
16 #include <iostream>
17 #include <algorithm>
18 #include <numeric>
19 #include <string>
20 #include <tuple>
21 #include <type_traits>
22 #include <vector>
23 #include <functional>
24 #include "fibonacci.hpp"
25 #include "../../test/table_type.hpp"
26 #include "table_helper.hpp"
27 #include "performance.hpp"
28
29
30 using namespace std;
31
32 boost::multiprecision::cpp_int total_sum(0);
33
34 template <typename Func, class Table>
exec_timed_test_foo(Func f,const Table & data,double min_elapsed=0.5)35 double exec_timed_test_foo(Func f, const Table& data, double min_elapsed = 0.5)
36 {
37 double t = 0;
38 unsigned repeats = 1;
39 typename Table::value_type::first_type sum{0};
40 stopwatch<boost::chrono::high_resolution_clock> w;
41 do
42 {
43 for(unsigned count = 0; count < repeats; ++count)
44 {
45 for(typename Table::size_type n = 0; n < data.size(); ++n)
46 sum += f(data[n].first, data[n].second);
47 }
48
49 t = boost::chrono::duration_cast<boost::chrono::duration<double>>(w.elapsed()).count();
50 if(t < min_elapsed)
51 repeats *= 2;
52 }
53 while(t < min_elapsed);
54 total_sum += sum;
55 return t / repeats;
56 }
57
58
59 template <typename T>
60 struct test_function_template
61 {
62 vector<pair<T, T> > const & data;
63 const char* data_name;
64
test_function_templatetest_function_template65 test_function_template(vector<pair<T, T> > const &data, const char* name) : data(data), data_name(name) {}
66
67 template <typename Function>
operator ()test_function_template68 void operator()(pair<Function, string> const &f) const
69 {
70 auto result = exec_timed_test_foo(f.first, data);
71 auto table_name = string("gcd method comparison with ") + compiler_name() + string(" on ") + platform_name();
72
73 report_execution_time(result,
74 table_name,
75 string(data_name),
76 string(f.second) + "\n" + boost_name());
77 }
78 };
79
80 boost::random::mt19937 rng;
81 boost::random::uniform_int_distribution<> d_0_6(0, 6);
82 boost::random::uniform_int_distribution<> d_1_20(1, 20);
83
84 template <class T>
get_prime_products()85 T get_prime_products()
86 {
87 int n_primes = d_0_6(rng);
88 switch(n_primes)
89 {
90 case 0:
91 // Generate a power of 2:
92 return static_cast<T>(1u) << d_1_20(rng);
93 case 1:
94 // prime number:
95 return boost::math::prime(d_1_20(rng) + 3);
96 }
97 T result = 1;
98 for(int i = 0; i < n_primes; ++i)
99 result *= boost::math::prime(d_1_20(rng) + 3) * boost::math::prime(d_1_20(rng) + 3) * boost::math::prime(d_1_20(rng) + 3) * boost::math::prime(d_1_20(rng) + 3) * boost::math::prime(d_1_20(rng) + 3);
100 return result;
101 }
102
103 template <class T>
get_uniform_random()104 T get_uniform_random()
105 {
106 static boost::random::uniform_int_distribution<T> minimax((std::numeric_limits<T>::min)(), (std::numeric_limits<T>::max)());
107 return minimax(rng);
108 }
109
110 template <class T>
even(T const & val)111 inline bool even(T const& val)
112 {
113 return !(val & 1u);
114 }
115
116 template <class Backend, boost::multiprecision::expression_template_option ExpressionTemplates>
even(boost::multiprecision::number<Backend,ExpressionTemplates> const & val)117 inline bool even(boost::multiprecision::number<Backend, ExpressionTemplates> const& val)
118 {
119 return !bit_test(val, 0);
120 }
121
122 template <class T>
euclid_textbook(T a,T b)123 T euclid_textbook(T a, T b)
124 {
125 using std::swap;
126 if(a < b)
127 swap(a, b);
128 while(b)
129 {
130 T t = b;
131 b = a % b;
132 a = t;
133 }
134 return a;
135 }
136
137 template <class T>
binary_textbook(T u,T v)138 T binary_textbook(T u, T v)
139 {
140 if(u && v)
141 {
142 unsigned shifts = (std::min)(boost::multiprecision::lsb(u), boost::multiprecision::lsb(v));
143 if(shifts)
144 {
145 u >>= shifts;
146 v >>= shifts;
147 }
148 while(u)
149 {
150 unsigned bit_index = boost::multiprecision::lsb(u);
151 if(bit_index)
152 {
153 u >>= bit_index;
154 }
155 else if(bit_index = boost::multiprecision::lsb(v))
156 {
157 v >>= bit_index;
158 }
159 else
160 {
161 if(u < v)
162 v = (v - u) >> 1u;
163 else
164 u = (u - v) >> 1u;
165 }
166 }
167 return v << shifts;
168 }
169 return u + v;
170 }
171
172 template <typename Integer>
gcd_default(Integer a,Integer b)173 inline BOOST_CXX14_CONSTEXPR Integer gcd_default(Integer a, Integer b) BOOST_GCD_NOEXCEPT(Integer)
174 {
175 using boost::integer::gcd;
176 return gcd(a, b);
177 }
178
179
180 template <class T>
test_type(const char * name)181 void test_type(const char* name)
182 {
183 using namespace boost::integer::gcd_detail;
184 typedef T int_type;
185 std::vector<pair<int_type, int_type> > data;
186
187 for(unsigned i = 0; i < 1000; ++i)
188 {
189 data.push_back(std::make_pair(get_prime_products<T>(), get_prime_products<T>()));
190 }
191 std::string row_name("gcd<");
192 row_name += name;
193 row_name += "> (random prime number products)";
194
195 typedef pair< function<int_type(int_type, int_type)>, string> f_test;
196 array<f_test, 6> test_functions{ {
197 { gcd_default<int_type>, "gcd" },
198 { Euclid_gcd<int_type>, "Euclid_gcd" },
199 { Stein_gcd<int_type>, "Stein_gcd" } ,
200 { mixed_binary_gcd<int_type>, "mixed_binary_gcd" },
201 { binary_textbook<int_type>, "Stein_gcd_textbook" },
202 { euclid_textbook<int_type>, "gcd_euclid_textbook" },
203 } };
204 for_each(begin(test_functions), end(test_functions), test_function_template<int_type>(data, row_name.c_str()));
205
206 data.clear();
207 for(unsigned i = 0; i < 1000; ++i)
208 {
209 data.push_back(std::make_pair(get_uniform_random<T>(), get_uniform_random<T>()));
210 }
211 row_name.erase();
212 row_name += "gcd<";
213 row_name += name;
214 row_name += "> (uniform random numbers)";
215 for_each(begin(test_functions), end(test_functions), test_function_template<int_type>(data, row_name.c_str()));
216
217 // Fibonacci number tests:
218 row_name.erase();
219 row_name += "gcd<";
220 row_name += name;
221 row_name += "> (adjacent Fibonacci numbers)";
222 for_each(begin(test_functions), end(test_functions), test_function_template<int_type>(fibonacci_numbers_permution_1<T>(), row_name.c_str()));
223
224 row_name.erase();
225 row_name += "gcd<";
226 row_name += name;
227 row_name += "> (permutations of Fibonacci numbers)";
228 for_each(begin(test_functions), end(test_functions), test_function_template<int_type>(fibonacci_numbers_permution_2<T>(), row_name.c_str()));
229
230 row_name.erase();
231 row_name += "gcd<";
232 row_name += name;
233 row_name += "> (Trivial cases)";
234 for_each(begin(test_functions), end(test_functions), test_function_template<int_type>(trivial_gcd_test_cases<T>(), row_name.c_str()));
235 }
236
237 /*******************************************************************************************************************/
238
239 template <class T>
generate_random(unsigned bits_wanted)240 T generate_random(unsigned bits_wanted)
241 {
242 static boost::random::mt19937 gen;
243 typedef boost::random::mt19937::result_type random_type;
244
245 T max_val;
246 unsigned digits;
247 if(std::numeric_limits<T>::is_bounded && (bits_wanted == (unsigned)std::numeric_limits<T>::digits))
248 {
249 max_val = (std::numeric_limits<T>::max)();
250 digits = std::numeric_limits<T>::digits;
251 }
252 else
253 {
254 max_val = T(1) << bits_wanted;
255 digits = bits_wanted;
256 }
257
258 unsigned bits_per_r_val = std::numeric_limits<random_type>::digits - 1;
259 while((random_type(1) << bits_per_r_val) > (gen.max)()) --bits_per_r_val;
260
261 unsigned terms_needed = digits / bits_per_r_val + 1;
262
263 T val = 0;
264 for(unsigned i = 0; i < terms_needed; ++i)
265 {
266 val *= (gen.max)();
267 val += gen();
268 }
269 val %= max_val;
270 return val;
271 }
272
273 template <typename N>
gcd_stein(N m,N n)274 N gcd_stein(N m, N n)
275 {
276 BOOST_ASSERT(m >= static_cast<N>(0));
277 BOOST_ASSERT(n >= static_cast<N>(0));
278 if(m == N(0)) return n;
279 if(n == N(0)) return m;
280 // m > 0 && n > 0
281 unsigned d_m = 0;
282 while(even(m)) { m >>= 1; d_m++; }
283 unsigned d_n = 0;
284 while(even(n)) { n >>= 1; d_n++; }
285 // odd(m) && odd(n)
286 while(m != n) {
287 if(n > m) swap(n, m);
288 m -= n;
289 do m >>= 1; while(even(m));
290 // m == n
291 }
292 return m << (std::min)(d_m, d_n);
293 }
294
295
big_gcd(const boost::multiprecision::cpp_int & a,const boost::multiprecision::cpp_int & b)296 boost::multiprecision::cpp_int big_gcd(const boost::multiprecision::cpp_int& a, const boost::multiprecision::cpp_int& b)
297 {
298 return boost::multiprecision::gcd(a, b);
299 }
300
301 namespace boost { namespace multiprecision { namespace backends {
302
303 template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
304 inline typename enable_if_c<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
eval_gcd_new(cpp_int_backend<MinBits1,MaxBits1,SignType1,Checked1,Allocator1> & result,const cpp_int_backend<MinBits1,MaxBits1,SignType1,Checked1,Allocator1> & a,const cpp_int_backend<MinBits1,MaxBits1,SignType1,Checked1,Allocator1> & b)305 eval_gcd_new(
306 cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
307 const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a,
308 const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& b)
309 {
310 using default_ops::eval_lsb;
311 using default_ops::eval_is_zero;
312 using default_ops::eval_get_sign;
313
314 if(a.size() == 1)
315 {
316 eval_gcd(result, b, *a.limbs());
317 return;
318 }
319 if(b.size() == 1)
320 {
321 eval_gcd(result, a, *b.limbs());
322 return;
323 }
324
325 int shift;
326
327 cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> u(a), v(b), mod;
328
329 int s = eval_get_sign(u);
330
331 /* GCD(0,x) := x */
332 if(s < 0)
333 {
334 u.negate();
335 }
336 else if(s == 0)
337 {
338 result = v;
339 return;
340 }
341 s = eval_get_sign(v);
342 if(s < 0)
343 {
344 v.negate();
345 }
346 else if(s == 0)
347 {
348 result = u;
349 return;
350 }
351
352 /* Let shift := lg K, where K is the greatest power of 2
353 dividing both u and v. */
354
355 unsigned us = eval_lsb(u);
356 unsigned vs = eval_lsb(v);
357 shift = (std::min)(us, vs);
358 eval_right_shift(u, us);
359 eval_right_shift(v, vs);
360
361 // From now on access u and v via pointers, that way we have a trivial swap:
362 cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>* up(&u), *vp(&v), *mp(&mod);
363
364 do
365 {
366 /* Now u and v are both odd, so diff(u, v) is even.
367 Let u = min(u, v), v = diff(u, v)/2. */
368 s = up->compare(*vp);
369 if(s > 0)
370 std::swap(up, vp);
371 if(s == 0)
372 break;
373 if(vp->size() <= 2)
374 {
375 if(vp->size() == 1)
376 *up = boost::integer::gcd_detail::mixed_binary_gcd(*vp->limbs(), *up->limbs());
377 else
378 {
379 double_limb_type i, j;
380 i = vp->limbs()[0] | (static_cast<double_limb_type>(vp->limbs()[1]) << sizeof(limb_type) * CHAR_BIT);
381 j = (up->size() == 1) ? *up->limbs() : up->limbs()[0] | (static_cast<double_limb_type>(up->limbs()[1]) << sizeof(limb_type) * CHAR_BIT);
382 u = boost::integer::gcd_detail::mixed_binary_gcd(i, j);
383 }
384 break;
385 }
386 if(vp->size() > up->size() /*eval_msb(*vp) > eval_msb(*up) + 32*/)
387 {
388 eval_modulus(*mp, *vp, *up);
389 std::swap(vp, mp);
390 eval_subtract(*up, *vp);
391 if(eval_is_zero(*vp) == 0)
392 {
393 vs = eval_lsb(*vp);
394 eval_right_shift(*vp, vs);
395 }
396 else
397 break;
398 if(eval_is_zero(*up) == 0)
399 {
400 vs = eval_lsb(*up);
401 eval_right_shift(*up, vs);
402 }
403 else
404 {
405 std::swap(up, vp);
406 break;
407 }
408 }
409 else
410 {
411 eval_subtract(*vp, *up);
412 vs = eval_lsb(*vp);
413 eval_right_shift(*vp, vs);
414 }
415 }
416 while(true);
417
418 result = *up;
419 eval_left_shift(result, shift);
420 }
421
422 }}}
423
424
big_gcd_new(const boost::multiprecision::cpp_int & a,const boost::multiprecision::cpp_int & b)425 boost::multiprecision::cpp_int big_gcd_new(const boost::multiprecision::cpp_int& a, const boost::multiprecision::cpp_int& b)
426 {
427 boost::multiprecision::cpp_int result;
428 boost::multiprecision::backends::eval_gcd_new(result.backend(), a.backend(), b.backend());
429 return result;
430 }
431
432 #if 0
433 void test_n_bits(unsigned n, std::string data_name, const std::vector<pair<boost::multiprecision::cpp_int, boost::multiprecision::cpp_int> >* p_data = 0)
434 {
435 using namespace boost::math::detail;
436 typedef boost::multiprecision::cpp_int int_type;
437 std::vector<pair<int_type, int_type> > data, data2;
438
439 for(unsigned i = 0; i < 1000; ++i)
440 {
441 data.push_back(std::make_pair(generate_random<int_type>(n), generate_random<int_type>(n)));
442 }
443
444 typedef pair< function<int_type(int_type, int_type)>, string> f_test;
445 array<f_test, 2> test_functions{ { /*{ Stein_gcd<int_type>, "Stein_gcd" } ,{ Euclid_gcd<int_type>, "Euclid_gcd" },{ binary_textbook<int_type>, "Stein_gcd_textbook" },{ euclid_textbook<int_type>, "gcd_euclid_textbook" },{ mixed_binary_gcd<int_type>, "mixed_binary_gcd" },{ gcd_stein<int_type>, "gcd_stein" },*/{ big_gcd, "boost::multiprecision::gcd" },{ big_gcd_new, "big_gcd_new" } } };
446 for_each(begin(test_functions), end(test_functions), test_function_template<int_type>(p_data ? *p_data : data, data_name.c_str(), true));
447 }
448 #endif
449
main()450 int main()
451 {
452 test_type<unsigned short>("unsigned short");
453 test_type<unsigned>("unsigned");
454 test_type<unsigned long>("unsigned long");
455 test_type<unsigned long long>("unsigned long long");
456
457 test_type<boost::multiprecision::uint256_t>("boost::multiprecision::uint256_t");
458 test_type<boost::multiprecision::uint512_t>("boost::multiprecision::uint512_t");
459 test_type<boost::multiprecision::uint1024_t>("boost::multiprecision::uint1024_t");
460
461 /*
462 test_n_bits(16, " 16 bit random values");
463 test_n_bits(32, " 32 bit random values");
464 test_n_bits(64, " 64 bit random values");
465 test_n_bits(125, " 125 bit random values");
466 test_n_bits(250, " 250 bit random values");
467 test_n_bits(500, " 500 bit random values");
468 test_n_bits(1000, " 1000 bit random values");
469 test_n_bits(5000, " 5000 bit random values");
470 test_n_bits(10000, "10000 bit random values");
471 //test_n_bits(100000);
472 //test_n_bits(1000000);
473
474 test_n_bits(0, "consecutive first 1000 fibonacci numbers", &fibonacci_numbers_cpp_int_permution_1());
475 test_n_bits(0, "permutations of first 1000 fibonacci numbers", &fibonacci_numbers_cpp_int_permution_2());
476 */
477 return 0;
478 }
479