ceph/src/boost/libs/math/reporting/performance/test_gcd.cpp

   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/math/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>
  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
  65     test_function_template(vector<pair<T, T> > const &data, const char* name) : data(data), data_name(name) {}
  66
  67     template <typename Function>
  68     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>
  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>
 104 T get_uniform_random()
 105 {
 106    static boost::random::uniform_int_distribution<T> minmax(std::numeric_limits<T>::min(), std::numeric_limits<T>::max());
 107    return minmax(rng);
 108 }
 109
 110 template <class T>
 111 inline bool even(T const& val)
 112 {
 113    return !(val & 1u);
 114 }
 115
 116 template <class Backend, boost::multiprecision::expression_template_option ExpressionTemplates>
 117 inline bool even(boost::multiprecision::number<Backend, ExpressionTemplates> const& val)
 118 {
 119    return !bit_test(val, 0);
 120 }
 121
 122 template <class T>
 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>
 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 //
 173 // The Mixed Binary Euclid Algorithm
 174 // Sidi Mohamed Sedjelmaci
 175 // Electronic Notes in Discrete Mathematics 35 (2009) 169\96176
 176 //
 177 template <class T>
 178 T mixed_binary_gcd(T u, T v)
 179 {
 180    using std::swap;
 181    if(u < v)
 182       swap(u, v);
 183
 184    unsigned shifts = 0;
 185
 186    if(!u)
 187       return v;
 188    if(!v)
 189       return u;
 190
 191    while(even(u) && even(v))
 192    {
 193       u >>= 1u;
 194       v >>= 1u;
 195       ++shifts;
 196    }
 197
 198    while(v > 1)
 199    {
 200       u %= v;
 201       v -= u;
 202       if(!u)
 203          return v << shifts;
 204       if(!v)
 205          return u << shifts;
 206       while(even(u)) u >>= 1u;
 207       while(even(v)) v >>= 1u;
 208       if(u < v)
 209          swap(u, v);
 210    }
 211    return (v == 1 ? v : u) << shifts;
 212 }
 213
 214 template <class T>
 215 void test_type(const char* name)
 216 {
 217    using namespace boost::math::detail;
 218    typedef T int_type;
 219    std::vector<pair<int_type, int_type> > data;
 220
 221    for(unsigned i = 0; i < 1000; ++i)
 222    {
 223       data.push_back(std::make_pair(get_prime_products<T>(), get_prime_products<T>()));
 224    }
 225    std::string row_name("gcd<");
 226    row_name += name;
 227    row_name += "> (random prime number products)";
 228
 229    typedef pair< function<int_type(int_type, int_type)>, string> f_test;
 230    array<f_test, 5> test_functions{ {
 231       { Stein_gcd<int_type>, "Stein_gcd" } ,
 232       { Euclid_gcd<int_type>, "Euclid_gcd" },
 233       { binary_textbook<int_type>, "Stein_gcd_textbook" },
 234       { euclid_textbook<int_type>, "gcd_euclid_textbook" },
 235       { mixed_binary_gcd<int_type>, "mixed_binary_gcd" }
 236    } };
 237    for_each(begin(test_functions), end(test_functions), test_function_template<int_type>(data, row_name.c_str()));
 238
 239    data.clear();
 240    for(unsigned i = 0; i < 1000; ++i)
 241    {
 242       data.push_back(std::make_pair(get_uniform_random<T>(), get_uniform_random<T>()));
 243    }
 244    row_name.erase();
 245    row_name += "gcd<";
 246    row_name += name;
 247    row_name += "> (uniform random numbers)";
 248    for_each(begin(test_functions), end(test_functions), test_function_template<int_type>(data, row_name.c_str()));
 249
 250    // Fibonacci number tests:
 251    row_name.erase();
 252    row_name += "gcd<";
 253    row_name += name;
 254    row_name += "> (adjacent Fibonacci numbers)";
 255    for_each(begin(test_functions), end(test_functions), test_function_template<int_type>(fibonacci_numbers_permution_1<T>(), row_name.c_str()));
 256
 257    row_name.erase();
 258    row_name += "gcd<";
 259    row_name += name;
 260    row_name += "> (permutations of Fibonacci numbers)";
 261    for_each(begin(test_functions), end(test_functions), test_function_template<int_type>(fibonacci_numbers_permution_2<T>(), row_name.c_str()));
 262
 263    row_name.erase();
 264    row_name += "gcd<";
 265    row_name += name;
 266    row_name += "> (Trivial cases)";
 267    for_each(begin(test_functions), end(test_functions), test_function_template<int_type>(trivial_gcd_test_cases<T>(), row_name.c_str()));
 268 }
 269
 270 /*******************************************************************************************************************/
 271
 272 template <class T>
 273 T generate_random(unsigned bits_wanted)
 274 {
 275    static boost::random::mt19937 gen;
 276    typedef boost::random::mt19937::result_type random_type;
 277
 278    T max_val;
 279    unsigned digits;
 280    if(std::numeric_limits<T>::is_bounded && (bits_wanted == (unsigned)std::numeric_limits<T>::digits))
 281    {
 282       max_val = (std::numeric_limits<T>::max)();
 283       digits = std::numeric_limits<T>::digits;
 284    }
 285    else
 286    {
 287       max_val = T(1) << bits_wanted;
 288       digits = bits_wanted;
 289    }
 290
 291    unsigned bits_per_r_val = std::numeric_limits<random_type>::digits - 1;
 292    while((random_type(1) << bits_per_r_val) > (gen.max)()) --bits_per_r_val;
 293
 294    unsigned terms_needed = digits / bits_per_r_val + 1;
 295
 296    T val = 0;
 297    for(unsigned i = 0; i < terms_needed; ++i)
 298    {
 299       val *= (gen.max)();
 300       val += gen();
 301    }
 302    val %= max_val;
 303    return val;
 304 }
 305
 306 template <typename N>
 307 N gcd_stein(N m, N n)
 308 {
 309    BOOST_ASSERT(m >= static_cast<N>(0));
 310    BOOST_ASSERT(n >= static_cast<N>(0));
 311    if(m == N(0)) return n;
 312    if(n == N(0)) return m;
 313            // m > 0 && n > 0
 314    unsigned d_m = 0;
 315    while(even(m)) { m >>= 1; d_m++; }
 316    unsigned d_n = 0;
 317    while(even(n)) { n >>= 1; d_n++; }
 318            // odd(m) && odd(n)
 319       while(m != n) {
 320       if(n > m) swap(n, m);
 321       m -= n;
 322       do m >>= 1; while(even(m));
 323                   // m == n
 324    }
 325    return m << std::min(d_m, d_n);
 326 }
 327
 328 boost::multiprecision::cpp_int big_gcd(const boost::multiprecision::cpp_int& a, const boost::multiprecision::cpp_int& b)
 329 {
 330    return boost::multiprecision::gcd(a, b);
 331 }
 332
 333 namespace boost { namespace multiprecision { namespace backends {
 334
 335 template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
 336 inline typename enable_if_c<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
 337    eval_gcd_new(
 338       cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
 339       const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a,
 340       const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& b)
 341 {
 342    using default_ops::eval_lsb;
 343    using default_ops::eval_is_zero;
 344    using default_ops::eval_get_sign;
 345
 346    if(a.size() == 1)
 347    {
 348       eval_gcd(result, b, *a.limbs());
 349       return;
 350    }
 351    if(b.size() == 1)
 352    {
 353       eval_gcd(result, a, *b.limbs());
 354       return;
 355    }
 356
 357    int shift;
 358
 359    cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> u(a), v(b), mod;
 360
 361    int s = eval_get_sign(u);
 362
 363    /* GCD(0,x) := x */
 364    if(s < 0)
 365    {
 366       u.negate();
 367    }
 368    else if(s == 0)
 369    {
 370       result = v;
 371       return;
 372    }
 373    s = eval_get_sign(v);
 374    if(s < 0)
 375    {
 376       v.negate();
 377    }
 378    else if(s == 0)
 379    {
 380       result = u;
 381       return;
 382    }
 383
 384    /* Let shift := lg K, where K is the greatest power of 2
 385    dividing both u and v. */
 386
 387    unsigned us = eval_lsb(u);
 388    unsigned vs = eval_lsb(v);
 389    shift = (std::min)(us, vs);
 390    eval_right_shift(u, us);
 391    eval_right_shift(v, vs);
 392
 393    // From now on access u and v via pointers, that way we have a trivial swap:
 394    cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>* up(&u), *vp(&v), *mp(&mod);
 395
 396    do
 397    {
 398       /* Now u and v are both odd, so diff(u, v) is even.
 399       Let u = min(u, v), v = diff(u, v)/2. */
 400       s = up->compare(*vp);
 401       if(s > 0)
 402          std::swap(up, vp);
 403       if(s == 0)
 404          break;
 405       if(vp->size() <= 2)
 406       {
 407          if(vp->size() == 1)
 408             *up = mixed_binary_gcd(*vp->limbs(), *up->limbs());
 409          else
 410          {
 411             double_limb_type i, j;
 412             i = vp->limbs()[0] | (static_cast<double_limb_type>(vp->limbs()[1]) << sizeof(limb_type) * CHAR_BIT);
 413             j = (up->size() == 1) ? *up->limbs() : up->limbs()[0] | (static_cast<double_limb_type>(up->limbs()[1]) << sizeof(limb_type) * CHAR_BIT);
 414             u = mixed_binary_gcd(i, j);
 415          }
 416          break;
 417       }
 418       if(vp->size() > up->size() /*eval_msb(*vp) > eval_msb(*up) + 32*/)
 419       {
 420          eval_modulus(*mp, *vp, *up);
 421          std::swap(vp, mp);
 422          eval_subtract(*up, *vp);
 423          if(eval_is_zero(*vp) == 0)
 424          {
 425             vs = eval_lsb(*vp);
 426             eval_right_shift(*vp, vs);
 427          }
 428          else
 429             break;
 430          if(eval_is_zero(*up) == 0)
 431          {
 432             vs = eval_lsb(*up);
 433             eval_right_shift(*up, vs);
 434          }
 435          else
 436          {
 437             std::swap(up, vp);
 438             break;
 439          }
 440       }
 441       else
 442       {
 443          eval_subtract(*vp, *up);
 444          vs = eval_lsb(*vp);
 445          eval_right_shift(*vp, vs);
 446       }
 447    }
 448    while(true);
 449
 450    result = *up;
 451    eval_left_shift(result, shift);
 452 }
 453
 454 }}}
 455
 456
 457 boost::multiprecision::cpp_int big_gcd_new(const boost::multiprecision::cpp_int& a, const boost::multiprecision::cpp_int& b)
 458 {
 459    boost::multiprecision::cpp_int result;
 460    boost::multiprecision::backends::eval_gcd_new(result.backend(), a.backend(), b.backend());
 461    return result;
 462 }
 463
 464 #if 0
 465 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)
 466 {
 467    using namespace boost::math::detail;
 468    typedef boost::multiprecision::cpp_int int_type;
 469    std::vector<pair<int_type, int_type> > data, data2;
 470
 471    for(unsigned i = 0; i < 1000; ++i)
 472    {
 473       data.push_back(std::make_pair(generate_random<int_type>(n), generate_random<int_type>(n)));
 474    }
 475
 476    typedef pair< function<int_type(int_type, int_type)>, string> f_test;
 477    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" } } };
 478    for_each(begin(test_functions), end(test_functions), test_function_template<int_type>(p_data ? *p_data : data, data_name.c_str(), true));
 479 }
 480 #endif
 481
 482 int main()
 483 {
 484     test_type<unsigned short>("unsigned short");
 485     test_type<unsigned>("unsigned");
 486     test_type<unsigned long>("unsigned long");
 487     test_type<unsigned long long>("unsigned long long");
 488
 489     test_type<boost::multiprecision::uint256_t>("boost::multiprecision::uint256_t");
 490     test_type<boost::multiprecision::uint512_t>("boost::multiprecision::uint512_t");
 491     test_type<boost::multiprecision::uint1024_t>("boost::multiprecision::uint1024_t");
 492
 493     /*
 494     test_n_bits(16, "   16 bit random values");
 495     test_n_bits(32, "   32 bit random values");
 496     test_n_bits(64, "   64 bit random values");
 497     test_n_bits(125, "  125 bit random values");
 498     test_n_bits(250, "  250 bit random values");
 499     test_n_bits(500, "  500 bit random values");
 500     test_n_bits(1000, " 1000 bit random values");
 501     test_n_bits(5000, " 5000 bit random values");
 502     test_n_bits(10000, "10000 bit random values");
 503     //test_n_bits(100000);
 504     //test_n_bits(1000000);
 505
 506     test_n_bits(0, "consecutive first 1000 fibonacci numbers", &fibonacci_numbers_cpp_int_permution_1());
 507     test_n_bits(0, "permutations of first 1000 fibonacci numbers", &fibonacci_numbers_cpp_int_permution_2());
 508     */
 509 }