1 ///////////////////////////////////////////////////////////////////
2 // Copyright Christopher Kormanyos 2019. //
3 // Distributed under the Boost Software License, //
4 // Version 1.0. (See accompanying file LICENSE_1_0.txt //
5 // or copy at http://www.boost.org/LICENSE_1_0.txt) //
6 ///////////////////////////////////////////////////////////////////
7
8 #include <algorithm>
9 #include <cstddef>
10 #include <cstdint>
11 #include <ctime>
12 #include <iomanip>
13 #include <iostream>
14 #include <limits>
15 #include <iterator>
16 #include <vector>
17
18 #include <boost/multiprecision/cpp_int.hpp>
19
20 class random_pcg32_fast_base
21 {
22 protected:
23 using itype = std::uint64_t;
24
25 static constexpr bool is_mcg = false;
26
27 public:
28 virtual ~random_pcg32_fast_base() = default;
29
30 protected:
random_pcg32_fast_base(const itype=itype ())31 explicit random_pcg32_fast_base(const itype = itype()) { }
32
33 random_pcg32_fast_base(const random_pcg32_fast_base&) = default;
34
35 random_pcg32_fast_base& operator=(const random_pcg32_fast_base&) = default;
36
37 template<typename ArithmeticType>
rotr(const ArithmeticType & value_being_shifted,const std::size_t bits_to_shift)38 static ArithmeticType rotr(const ArithmeticType& value_being_shifted,
39 const std::size_t bits_to_shift)
40 {
41 const std::size_t left_shift_amount =
42 std::numeric_limits<ArithmeticType>::digits - bits_to_shift;
43
44 const ArithmeticType part1 = ((bits_to_shift > 0U) ? ArithmeticType(value_being_shifted >> bits_to_shift) : value_being_shifted);
45 const ArithmeticType part2 = ((bits_to_shift > 0U) ? ArithmeticType(value_being_shifted << left_shift_amount) : 0U);
46
47 return ArithmeticType(part1 | part2);
48 }
49
50 template<typename xtype,
51 typename itype>
52 struct xsh_rr_mixin
53 {
outputrandom_pcg32_fast_base::xsh_rr_mixin54 static xtype output(const itype internal_value)
55 {
56 using bitcount_t = std::size_t;
57
58 constexpr bitcount_t bits = bitcount_t(sizeof(itype) * 8U);
59 constexpr bitcount_t xtypebits = bitcount_t(sizeof(xtype) * 8U);
60 constexpr bitcount_t sparebits = bits - xtypebits;
61 constexpr bitcount_t wantedopbits = ((xtypebits >= 128U) ? 7U
62 : ((xtypebits >= 64U) ? 6U
63 : ((xtypebits >= 32U) ? 5U
64 : ((xtypebits >= 16U) ? 4U
65 : 3U))));
66
67 constexpr bitcount_t opbits = ((sparebits >= wantedopbits) ? wantedopbits : sparebits);
68 constexpr bitcount_t amplifier = wantedopbits - opbits;
69 constexpr bitcount_t mask = (1ULL << opbits) - 1U;
70 constexpr bitcount_t topspare = opbits;
71 constexpr bitcount_t bottomspare = sparebits - topspare;
72 constexpr bitcount_t xshift = (topspare + xtypebits) / 2U;
73
74 const bitcount_t rot =
75 ((opbits != 0U) ? (bitcount_t(internal_value >> (bits - opbits)) & mask)
76 : 0U);
77
78 const bitcount_t amprot = (rot << amplifier) & mask;
79
80 const itype internal_value_xor = internal_value ^ itype(internal_value >> xshift);
81
82 const xtype result = rotr(xtype(internal_value_xor >> bottomspare), amprot);
83
84 return result;
85 }
86 };
87 };
88
89 class random_pcg32_fast : public random_pcg32_fast_base
90 {
91 private:
92 static constexpr itype default_multiplier = static_cast<itype>(6364136223846793005ULL);
93 static constexpr itype default_increment = static_cast<itype>(1442695040888963407ULL);
94
95 public:
96 using result_type = std::uint32_t;
97
98 static constexpr itype default_seed = static_cast<itype>(0xCAFEF00DD15EA5E5ULL);
99
random_pcg32_fast(const itype state=default_seed)100 explicit random_pcg32_fast(const itype state = default_seed)
101 : random_pcg32_fast_base(state),
102 my_inc (default_increment),
103 my_state(is_mcg ? state | itype(3U) : bump(state + increment())) { }
104
random_pcg32_fast(const random_pcg32_fast & other)105 random_pcg32_fast(const random_pcg32_fast& other)
106 : random_pcg32_fast_base(other),
107 my_inc (other.my_inc),
108 my_state(other.my_state) { }
109
110 virtual ~random_pcg32_fast() = default;
111
operator =(const random_pcg32_fast & other)112 random_pcg32_fast& operator=(const random_pcg32_fast& other)
113 {
114 static_cast<void>(random_pcg32_fast_base::operator=(other));
115
116 if(this != &other)
117 {
118 my_inc = other.my_inc;
119 my_state = other.my_state;
120 }
121
122 return *this;
123 }
124
seed(const itype state=default_seed)125 void seed(const itype state = default_seed)
126 {
127 my_inc = default_increment;
128
129 my_state = (is_mcg ? state | itype(3U) : bump(state + increment()));
130 }
131
operator ()()132 result_type operator()()
133 {
134 const result_type value =
135 xsh_rr_mixin<result_type, itype>::output(base_generate0());
136
137 return value;
138 }
139
140 private:
141 itype my_inc;
142 itype my_state;
143
multiplier() const144 itype multiplier() const
145 {
146 return default_multiplier;
147 }
148
increment() const149 itype increment() const
150 {
151 return default_increment;
152 }
153
bump(const itype state)154 itype bump(const itype state)
155 {
156 return itype(state * multiplier()) + increment();
157 }
158
base_generate0()159 itype base_generate0()
160 {
161 const itype old_state = my_state;
162
163 my_state = bump(my_state);
164
165 return old_state;
166 }
167 };
168
169
170 template<typename UnsignedIntegralIteratorType,
171 typename RandomEngineType>
get_random_big_uint(RandomEngineType & rng,UnsignedIntegralIteratorType it_out)172 void get_random_big_uint(RandomEngineType& rng, UnsignedIntegralIteratorType it_out)
173 {
174 using local_random_value_type = typename RandomEngineType::result_type;
175
176 using local_uint_type = typename std::iterator_traits<UnsignedIntegralIteratorType>::value_type;
177
178 constexpr std::size_t digits_of_uint___type = static_cast<std::size_t>(std::numeric_limits<local_uint_type>::digits);
179 constexpr std::size_t digits_of_random_type = static_cast<std::size_t>(std::numeric_limits<local_random_value_type>::digits);
180
181 local_random_value_type next_random = rng();
182
183 *it_out = next_random;
184
185 for(std::size_t i = digits_of_random_type; i < digits_of_uint___type; i += digits_of_random_type)
186 {
187 (*it_out) <<= digits_of_random_type;
188
189 next_random = rng();
190
191 (*it_out) |= next_random;
192 }
193 }
194
195 using big_uint_backend_type =
196 boost::multiprecision::cpp_int_backend<8192UL << 1U,
197 8192UL << 1U,
198 boost::multiprecision::unsigned_magnitude>;
199
200 using big_uint_type = boost::multiprecision::number<big_uint_backend_type>;
201
202 namespace local
203 {
204 std::vector<big_uint_type> a(1024U);
205 std::vector<big_uint_type> b(a.size());
206 }
207
main()208 int main()
209 {
210 random_pcg32_fast rng;
211
212 rng.seed(std::clock());
213
214 for(auto i = 0U; i < local::a.size(); ++i)
215 {
216 get_random_big_uint(rng, local::a.begin() + i);
217 get_random_big_uint(rng, local::b.begin() + i);
218 }
219
220 std::size_t count = 0U;
221
222 float total_time = 0.0F;
223
224 const std::clock_t start = std::clock();
225
226 do
227 {
228 const std::size_t index = count % local::a.size();
229
230 local::a[index] * local::b[index];
231
232 ++count;
233 }
234 while((total_time = (float(std::clock() - start) / float(CLOCKS_PER_SEC))) < 10.0F);
235
236 // Boost.Multiprecision 1.71
237 // bits: 16384, kops_per_sec: 4.7
238 // bits: 32768, kops_per_sec: 1.2
239 // bits: 65536, kops_per_sec: 0.3
240 // bits: 131072, kops_per_sec: 0.075
241 // bits: 262144, kops_per_sec: 0.019
242 // bits: 524288, kops_per_sec: 0.0047
243
244 // Boost.Multiprecision + kara mult
245 // bits: 16384, kops_per_sec: 4.8
246 // bits: 32768, kops_per_sec: 1.6
247 // bits: 65536, kops_per_sec: 0.5
248 // bits: 131072, kops_per_sec: 0.15
249 // bits: 262144, kops_per_sec: 0.043
250 // bits: 524288, kops_per_sec: 0.011
251
252 const float kops_per_sec = (float(count) / total_time) / 1000.0F;
253
254 std::cout << "bits: "
255 << std::numeric_limits<big_uint_type>::digits
256 << ", kops_per_sec: "
257 << kops_per_sec
258 << count << std::endl;
259 }
260