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 template <typename Integer>
 173 inline BOOST_CXX14_CONSTEXPR Integer gcd_default(Integer a, Integer b) BOOST_GCD_NOEXCEPT(Integer)
 174 {
 175    using boost::math::gcd;
 176    return gcd(a, b);
 177 }
 178
 179
 180 template <class T>
 181 void test_type(const char* name)
 182 {
 183    using namespace boost::math::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>
 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>
 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
 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
 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::math::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::math::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
 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
 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 }