ceph/src/boost/libs/multiprecision/performance/gcd_bench.cpp

   1 //  Copyright 2020 John Maddock. Distributed under the Boost
   2 //  Software License, Version 1.0. (See accompanying file
   3 //  LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
   4
   5 #include <iostream>
   6 #include <benchmark/benchmark.h>
   7 #include <boost/multiprecision/cpp_int.hpp>
   8 #include <boost/multiprecision/gmp.hpp>
   9 #include <boost/random.hpp>
  10 #include <cmath>
  11
  12 #include <immintrin.h>
  13
  14 using namespace boost::multiprecision;
  15 using namespace boost::random;
  16
  17 namespace boost {
  18    namespace multiprecision {
  19       namespace backends {
  20
  21       template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  22       inline BOOST_MP_CXX14_CONSTEXPR typename enable_if_c<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
  23       eval_gcd_old(
  24           cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>&       result,
  25           const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a,
  26           const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& b)
  27       {
  28          using default_ops::eval_get_sign;
  29          using default_ops::eval_is_zero;
  30          using default_ops::eval_lsb;
  31
  32          if (a.size() == 1)
  33          {
  34             eval_gcd(result, b, *a.limbs());
  35             return;
  36          }
  37          if (b.size() == 1)
  38          {
  39             eval_gcd(result, a, *b.limbs());
  40             return;
  41          }
  42
  43          cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> u(a), v(b);
  44
  45          int s = eval_get_sign(u);
  46
  47          /* GCD(0,x) := x */
  48          if (s < 0)
  49          {
  50             u.negate();
  51          }
  52          else if (s == 0)
  53          {
  54             result = v;
  55             return;
  56          }
  57          s = eval_get_sign(v);
  58          if (s < 0)
  59          {
  60             v.negate();
  61          }
  62          else if (s == 0)
  63          {
  64             result = u;
  65             return;
  66          }
  67
  68          /* Let shift := lg K, where K is the greatest power of 2
  69    dividing both u and v. */
  70
  71          unsigned us    = eval_lsb(u);
  72          unsigned vs    = eval_lsb(v);
  73          int      shift = (std::min)(us, vs);
  74          eval_right_shift(u, us);
  75          eval_right_shift(v, vs);
  76
  77          do
  78          {
  79             /* Now u and v are both odd, so diff(u, v) is even.
  80       Let u = min(u, v), v = diff(u, v)/2. */
  81             s = u.compare(v);
  82             if (s > 0)
  83                u.swap(v);
  84             if (s == 0)
  85                break;
  86
  87             while (((u.size() + 2 < v.size()) && (v.size() * 100 / u.size() > 105)) || ((u.size() <= 2) && (v.size() > 4)))
  88             {
  89                //
  90                // Speical case: if u and v differ considerably in size, then a Euclid step
  91                // is more efficient as we reduce v by several limbs in one go.
  92                // Unfortunately it requires an expensive long division:
  93                //
  94                eval_modulus(v, v, u);
  95                u.swap(v);
  96             }
  97             if (v.size() <= 2)
  98             {
  99                //
 100                // Special case: if v has no more than 2 limbs
 101                // then we can reduce u and v to a pair of integers and perform
 102                // direct integer gcd:
 103                //
 104                if (v.size() == 1)
 105                   u = eval_gcd(*v.limbs(), *u.limbs());
 106                else
 107                {
 108                   double_limb_type i = v.limbs()[0] | (static_cast<double_limb_type>(v.limbs()[1]) << sizeof(limb_type) * CHAR_BIT);
 109                   double_limb_type j = (u.size() == 1) ? *u.limbs() : u.limbs()[0] | (static_cast<double_limb_type>(u.limbs()[1]) << sizeof(limb_type) * CHAR_BIT);
 110                   u                  = eval_gcd(i, j);
 111                }
 112                break;
 113             }
 114             //
 115             // Regular binary gcd case:
 116             //
 117             eval_subtract(v, u);
 118             vs = eval_lsb(v);
 119             eval_right_shift(v, vs);
 120          } while (true);
 121
 122          result = u;
 123          eval_left_shift(result, shift);
 124       }
 125
 126       }
 127    }
 128 }
 129
 130 template <class T>
 131 std::tuple<std::vector<T>, std::vector<T>, std::vector<T> >& get_test_vector(unsigned bits)
 132 {
 133    static std::map<unsigned, std::tuple<std::vector<T>, std::vector<T>, std::vector<T> > > data;
 134
 135    std::tuple<std::vector<T>, std::vector<T>, std::vector<T> >& result = data[bits];
 136
 137    if (std::get<0>(result).size() == 0)
 138    {
 139       mt19937                     mt;
 140       uniform_int_distribution<T> ui(T(1) << (bits - 1), T(1) << bits);
 141
 142       std::vector<T>& a = std::get<0>(result);
 143       std::vector<T>& b = std::get<1>(result);
 144       std::vector<T>& c = std::get<2>(result);
 145
 146       for (unsigned i = 0; i < 1000; ++i)
 147       {
 148          a.push_back(ui(mt));
 149          b.push_back(ui(mt));
 150          if (b.back() > a.back())
 151             b.back().swap(a.back());
 152          c.push_back(0);
 153       }
 154    }
 155    return result;
 156 }
 157
 158 template <class T>
 159 std::vector<T>& get_test_vector_a(unsigned bits)
 160 {
 161    return std::get<0>(get_test_vector<T>(bits));
 162 }
 163 template <class T>
 164 std::vector<T>& get_test_vector_b(unsigned bits)
 165 {
 166    return std::get<1>(get_test_vector<T>(bits));
 167 }
 168 template <class T>
 169 std::vector<T>& get_test_vector_c(unsigned bits)
 170 {
 171    return std::get<2>(get_test_vector<T>(bits));
 172 }
 173
 174
 175 template <typename T>
 176 static void BM_gcd_old(benchmark::State& state)
 177 {
 178    int                         bits = state.range(0);
 179
 180    std::vector<T>& a = get_test_vector_a<T>(bits);
 181    std::vector<T>& b = get_test_vector_b<T>(bits);
 182    std::vector<T>& c = get_test_vector_c<T>(bits);
 183
 184    for (auto _ : state)
 185    {
 186       for (unsigned i = 0; i < a.size(); ++i)
 187          eval_gcd_old(c[i].backend(), a[i].backend(), b[i].backend());
 188    }
 189    state.SetComplexityN(bits);
 190 }
 191
 192 template <typename T>
 193 static void BM_gcd_current(benchmark::State& state)
 194 {
 195    int                         bits = state.range(0);
 196
 197    std::vector<T>& a = get_test_vector_a<T>(bits);
 198    std::vector<T>& b = get_test_vector_b<T>(bits);
 199    std::vector<T>& c = get_test_vector_c<T>(bits);
 200
 201    for (auto _ : state)
 202    {
 203       for (unsigned i = 0; i < a.size(); ++i)
 204          eval_gcd(c[i].backend(), a[i].backend(), b[i].backend());
 205    }
 206    state.SetComplexityN(bits);
 207 }
 208
 209 constexpr unsigned lower_range = 512;
 210 constexpr unsigned upper_range = 1 << 15;
 211
 212 BENCHMARK_TEMPLATE(BM_gcd_old, cpp_int)->RangeMultiplier(2)->Range(lower_range, upper_range)->Unit(benchmark::kMillisecond)->Complexity();
 213 BENCHMARK_TEMPLATE(BM_gcd_current, cpp_int)->RangeMultiplier(2)->Range(lower_range, upper_range)->Unit(benchmark::kMillisecond)->Complexity();
 214 BENCHMARK_TEMPLATE(BM_gcd_old, cpp_int)->RangeMultiplier(2)->Range(lower_range, upper_range)->Unit(benchmark::kMillisecond)->Complexity();
 215 BENCHMARK_TEMPLATE(BM_gcd_current, mpz_int)->RangeMultiplier(2)->Range(lower_range, upper_range)->Unit(benchmark::kMillisecond)->Complexity();
 216 BENCHMARK_TEMPLATE(BM_gcd_current, mpz_int)->RangeMultiplier(2)->Range(lower_range, upper_range)->Unit(benchmark::kMillisecond)->Complexity();
 217
 218 BENCHMARK_MAIN();