EigenRand  0.3.0
RandUtils.h
Go to the documentation of this file.
1 
12 #ifndef EIGENRAND_RAND_UTILS_H
13 #define EIGENRAND_RAND_UTILS_H
14 
15 #include <EigenRand/MorePacketMath.h>
16 #include <EigenRand/PacketFilter.h>
17 #include <EigenRand/PacketRandomEngine.h>
18 
19 namespace Eigen
20 {
21  namespace internal
22  {
23  template<typename Packet, typename Rng,
24  typename RngResult = typename std::remove_reference<Rng>::type::result_type,
25  Rand::RandomEngineType reType = Rand::GetRandomEngineType<
26  typename std::remove_reference<Rng>::type
27  >::value>
28  struct RawbitsMaker;
29 
30  template<typename PacketType, typename Rng>
31  struct UniformRealUtils;
32 
33  template<typename PacketType, typename Rng>
34  struct RandUtils : public UniformRealUtils<PacketType, Rng>
35  {
36  EIGEN_STRONG_INLINE PacketType balanced(Rng& rng)
37  {
38  return psub(pmul(this->zero_to_one(rng), pset1<PacketType>(2)), pset1<PacketType>(1));
39  }
40 
41  EIGEN_STRONG_INLINE PacketType nonzero_uniform_real(Rng& rng)
42  {
43  constexpr auto epsilon = std::numeric_limits<typename unpacket_traits<PacketType>::type>::epsilon() / 8;
44  return padd(this->uniform_real(rng), pset1<PacketType>(epsilon));
45  }
46  };
47 
48  EIGEN_STRONG_INLINE uint32_t collect_upper8bits(uint32_t a, uint32_t b, uint32_t c)
49  {
50  return ((a & 0xFF000000) >> 24) | ((b & 0xFF000000) >> 16) | ((c & 0xFF000000) >> 8);
51  }
52  }
53 }
54 
55 #ifdef EIGEN_VECTORIZE_AVX
56 #include <immintrin.h>
57 
58 namespace Eigen
59 {
60  namespace internal
61  {
62  template<typename Rng>
63  struct RawbitsMaker<Packet4i, Rng, Packet8i, Rand::RandomEngineType::packet>
64  {
65  EIGEN_STRONG_INLINE Packet4i rawbits(Rng& rng)
66  {
67  return rng.half();
68  }
69 
70  EIGEN_STRONG_INLINE Packet4i rawbits_34(Rng& rng)
71  {
72  return rng.half();
73  }
74 
75  EIGEN_STRONG_INLINE Packet4i rawbits_half(Rng& rng)
76  {
77  return rng.half();
78  }
79  };
80 
81  template<typename Rng>
82  struct RawbitsMaker<Packet8i, Rng, Packet4i, Rand::RandomEngineType::packet>
83  {
84  EIGEN_STRONG_INLINE Packet8i rawbits(Rng& rng)
85  {
86  return _mm256_insertf128_si256(_mm256_castsi128_si256(rng()), rng(), 1);
87  }
88 
89  EIGEN_STRONG_INLINE Packet8i rawbits_34(Rng& rng)
90  {
91  return _mm256_insertf128_si256(_mm256_castsi128_si256(rng()), rng(), 1);
92  }
93 
94  EIGEN_STRONG_INLINE Packet4i rawbits_half(Rng& rng)
95  {
96  return rng();
97  }
98  };
99 
100  template<typename Rng, typename RngResult>
101  struct RawbitsMaker<Packet8i, Rng, RngResult, Rand::RandomEngineType::scalar>
102  {
103  EIGEN_STRONG_INLINE Packet8i rawbits(Rng& rng)
104  {
105  if (sizeof(decltype(rng())) == 8)
106  {
107  return _mm256_set_epi64x(rng(), rng(), rng(), rng());
108  }
109  else
110  {
111  return _mm256_set_epi32(rng(), rng(), rng(), rng(),
112  rng(), rng(), rng(), rng());
113  }
114  }
115 
116  EIGEN_STRONG_INLINE Packet8i rawbits_34(Rng& rng)
117  {
118  Packet8i p;
119  if (sizeof(decltype(rng())) == 8)
120  {
121 #ifdef EIGEN_VECTORIZE_AVX2
122  p = _mm256_setr_epi64x(rng(), rng(), rng(), 0);
123  p = _mm256_permutevar8x32_epi32(p, _mm256_setr_epi32(0, 1, 2, 7, 3, 4, 5, 7));
124  p = _mm256_shuffle_epi8(p, _mm256_setr_epi8(
125  0, 1, 2, 3,
126  4, 5, 6, 7,
127  8, 9, 10, 11,
128  3, 7, 11, 11,
129  0, 1, 2, 3,
130  4, 5, 6, 7,
131  8, 9, 10, 11,
132  3, 7, 11, 11
133  ));
134 
135 #else
136  auto v = rng();
137  p = _mm256_setr_epi64x(rng(), v, rng(), v >> 32);
138  Packet4i p1, p2, o = _mm_setr_epi8(
139  0, 1, 2, 3,
140  4, 5, 6, 7,
141  8, 9, 10, 11,
142  3, 7, 11, 11);
143  split_two(p, p1, p2);
144  p = combine_two(_mm_shuffle_epi8(p1, o), _mm_shuffle_epi8(p2, o));
145 #endif
146  }
147  else
148  {
149  p = _mm256_setr_epi32(rng(), rng(), rng(), 0, rng(), rng(), rng(), 0);
150 #ifdef EIGEN_VECTORIZE_AVX2
151  p = _mm256_shuffle_epi8(p, _mm256_setr_epi8(
152  0, 1, 2, 3,
153  4, 5, 6, 7,
154  8, 9, 10, 11,
155  3, 7, 11, 11,
156  0, 1, 2, 3,
157  4, 5, 6, 7,
158  8, 9, 10, 11,
159  3, 7, 11, 11
160  ));
161 #else
162  Packet4i p1, p2, o = _mm_setr_epi8(
163  0, 1, 2, 3,
164  4, 5, 6, 7,
165  8, 9, 10, 11,
166  3, 7, 11, 11);
167  split_two(p, p1, p2);
168  p = combine_two(_mm_shuffle_epi8(p1, o), _mm_shuffle_epi8(p2, o));
169 #endif
170  }
171  return p;
172  }
173 
174  EIGEN_STRONG_INLINE Packet4i rawbits_half(Rng& rng)
175  {
176  if (sizeof(decltype(rng())) == 8)
177  {
178  return _mm_set_epi64x(rng(), rng());
179  }
180  else
181  {
182  return _mm_set_epi32(rng(), rng(), rng(), rng());
183  }
184  }
185  };
186 
187  template<typename Rng>
188  struct RawbitsMaker<Packet8i, Rng, Packet8i, Rand::RandomEngineType::packet>
189  {
190  EIGEN_STRONG_INLINE Packet8i rawbits(Rng& rng)
191  {
192  return rng();
193  }
194 
195  EIGEN_STRONG_INLINE Packet8i rawbits_34(Rng& rng)
196  {
197  return rng();
198  }
199 
200  EIGEN_STRONG_INLINE Packet4i rawbits_half(Rng& rng)
201  {
202  return rng.half();
203  }
204  };
205 
206 #ifndef EIGEN_VECTORIZE_AVX2
207  template<>
208  EIGEN_STRONG_INLINE Packet8f bit_to_ur_float<Packet8i>(const Packet8i& x)
209  {
210  const Packet4i lower = pset1<Packet4i>(0x7FFFFF),
211  upper = pset1<Packet4i>(127 << 23);
212  const Packet8f one = pset1<Packet8f>(1);
213 
214  Packet4i x1, x2;
215  split_two(x, x1, x2);
216 
217  return psub(reinterpret_to_float(
218  combine_two(por(pand(x1, lower), upper), por(pand(x2, lower), upper)
219  )), one);
220  }
221 
222  template<>
223  EIGEN_STRONG_INLINE Packet4d bit_to_ur_double<Packet8i>(const Packet8i& x)
224  {
225  const Packet4i lower = pseti64<Packet4i>(0xFFFFFFFFFFFFFull),
226  upper = pseti64<Packet4i>(1023ull << 52);
227  const Packet4d one = pset1<Packet4d>(1);
228 
229  Packet4i x1, x2;
230  split_two(x, x1, x2);
231 
232  return psub(reinterpret_to_double(
233  combine_two(por(pand(x1, lower), upper), por(pand(x2, lower), upper)
234  )), one);
235  }
236 #endif
237 
238  template<typename Rng>
239  struct UniformRealUtils<Packet8f, Rng> : public RawbitsMaker<Packet8i, Rng>
240  {
241  EIGEN_STRONG_INLINE Packet8f zero_to_one(Rng& rng)
242  {
243  return pdiv(_mm256_cvtepi32_ps(pand(this->rawbits(rng), pset1<Packet8i>(0x7FFFFFFF))),
244  pset1<Packet8f>(0x7FFFFFFF));
245  }
246 
247  EIGEN_STRONG_INLINE Packet8f uniform_real(Rng& rng)
248  {
249  return bit_to_ur_float(this->rawbits_34(rng));
250  }
251  };
252 
253  template<typename Rng>
254  struct UniformRealUtils<Packet4d, Rng> : public RawbitsMaker<Packet8i, Rng>
255  {
256  EIGEN_STRONG_INLINE Packet4d zero_to_one(Rng& rng)
257  {
258  return pdiv(_mm256_cvtepi32_pd(pand(this->rawbits_half(rng), pset1<Packet4i>(0x7FFFFFFF))),
259  pset1<Packet4d>(0x7FFFFFFF));
260  }
261 
262  EIGEN_STRONG_INLINE Packet4d uniform_real(Rng& rng)
263  {
264  return bit_to_ur_double(this->rawbits(rng));
265  }
266  };
267  }
268 }
269 #endif
270 
271 #ifdef EIGEN_VECTORIZE_SSE2
272 #include <xmmintrin.h>
273 
274 namespace Eigen
275 {
276  namespace internal
277  {
278  template<typename Rng, typename RngResult>
279  struct RawbitsMaker<Packet4i, Rng, RngResult, Rand::RandomEngineType::scalar>
280  {
281  EIGEN_STRONG_INLINE Packet4i rawbits(Rng& rng)
282  {
283  if (sizeof(RngResult) == 8)
284  {
285  return _mm_set_epi64x(rng(), rng());
286  }
287  else
288  {
289  return _mm_set_epi32(rng(), rng(), rng(), rng());
290  }
291  }
292 
293  EIGEN_STRONG_INLINE Packet4i rawbits_34(Rng& rng)
294  {
295  if (sizeof(RngResult) == 8)
296  {
297  return _mm_set_epi64x(rng(), rng());
298  }
299  else
300  {
301 #ifdef EIGEN_VECTORIZE_SSSE3
302  Packet4i p = _mm_setr_epi32(rng(), rng(), rng(), 0);
303  return _mm_shuffle_epi8(p, _mm_setr_epi8(
304  0, 1, 2, 3,
305  4, 5, 6, 7,
306  8, 9, 10, 11,
307  3, 7, 11, 11));
308 #else
309  return _mm_set_epi32(rng(), rng(), rng(), rng());
310 #endif
311  }
312  }
313 
314  EIGEN_STRONG_INLINE Packet4i rawbits_half(Rng& rng)
315  {
316  if (sizeof(decltype(rng())) == 8)
317  {
318  return _mm_set_epi64x(0, rng());
319  }
320  else
321  {
322  return _mm_setr_epi32(rng(), rng(), 0, 0);
323  }
324  }
325  };
326 
327  template<typename Rng>
328  struct RawbitsMaker<Packet4i, Rng, Packet4i, Rand::RandomEngineType::packet>
329  {
330  EIGEN_STRONG_INLINE Packet4i rawbits(Rng& rng)
331  {
332  return rng();
333  }
334 
335  EIGEN_STRONG_INLINE Packet4i rawbits_34(Rng& rng)
336  {
337  return rng();
338  }
339 
340  EIGEN_STRONG_INLINE Packet4i rawbits_half(Rng& rng)
341  {
342  return rng();
343  }
344  };
345 
346  template<typename Rng>
347  struct UniformRealUtils<Packet4f, Rng> : public RawbitsMaker<Packet4i, Rng>
348  {
349  EIGEN_STRONG_INLINE Packet4f zero_to_one(Rng& rng)
350  {
351  return pdiv(_mm_cvtepi32_ps(pand(this->rawbits(rng), pset1<Packet4i>(0x7FFFFFFF))),
352  pset1<Packet4f>(0x7FFFFFFF));
353  }
354 
355  EIGEN_STRONG_INLINE Packet4f uniform_real(Rng& rng)
356  {
357  return bit_to_ur_float(this->rawbits_34(rng));
358  }
359  };
360 
361  template<typename Rng>
362  struct UniformRealUtils<Packet2d, Rng> : public RawbitsMaker<Packet4i, Rng>
363  {
364  EIGEN_STRONG_INLINE Packet2d zero_to_one(Rng& rng)
365  {
366  return pdiv(_mm_cvtepi32_pd(pand(this->rawbits_half(rng), pset1<Packet4i>(0x7FFFFFFF))),
367  pset1<Packet2d>(0x7FFFFFFF));
368  }
369 
370  EIGEN_STRONG_INLINE Packet2d uniform_real(Rng& rng)
371  {
372  return bit_to_ur_double(this->rawbits(rng));
373  }
374  };
375  }
376 }
377 #endif
378 
379 namespace Eigen
380 {
381  namespace internal
382  {
383  template<typename Gen, typename _Scalar, typename Rng, bool _mutable = false>
384  struct scalar_rng_adaptor
385  {
386  static_assert(
387  Rand::IsScalarRandomEngine<
388  typename std::remove_reference<Rng>::type
389  >::value ||
390  Rand::IsPacketRandomEngine<
391  typename std::remove_reference<Rng>::type
392  >::value, "Rng must satisfy RandomNumberEngine");
393 
394  Gen gen;
395  Rng rng;
396 
397  scalar_rng_adaptor(const Rng& _rng) : rng{ _rng }
398  {
399  }
400 
401  template<typename _Gen>
402  scalar_rng_adaptor(const Rng& _rng, _Gen&& _gen) : gen{ _gen }, rng{ _rng }
403  {
404  }
405 
406  scalar_rng_adaptor(const scalar_rng_adaptor& o) = default;
407  scalar_rng_adaptor(scalar_rng_adaptor&& o) = default;
408 
409  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const _Scalar operator() () const
410  {
411  return gen(rng);
412  }
413 
414  template<typename Packet>
415  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Packet packetOp() const
416  {
417  return gen.template packetOp<Packet>(rng);
418  }
419  };
420 
421  template<typename Gen, typename _Scalar, typename Rng>
422  struct scalar_rng_adaptor<Gen, _Scalar, Rng, true>
423  {
424  static_assert(
425  Rand::IsScalarRandomEngine<
426  typename std::remove_reference<Rng>::type
427  >::value ||
428  Rand::IsPacketRandomEngine<
429  typename std::remove_reference<Rng>::type
430  >::value, "Rng must satisfy RandomNumberEngine");
431 
432  mutable Gen gen;
433  Rng rng;
434 
435  scalar_rng_adaptor(const Rng& _rng) : rng{ _rng }
436  {
437  }
438 
439  template<typename _Gen>
440  scalar_rng_adaptor(const Rng& _rng, _Gen&& _gen) : gen{ _gen }, rng{ _rng }
441  {
442  }
443 
444  scalar_rng_adaptor(const scalar_rng_adaptor& o) = default;
445  scalar_rng_adaptor(scalar_rng_adaptor&& o) = default;
446 
447  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const _Scalar operator() () const
448  {
449  return gen(rng);
450  }
451 
452  template<typename Packet>
453  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Packet packetOp() const
454  {
455  return gen.template packetOp<Packet>(rng);
456  }
457  };
458 
459  template<typename Gen, typename _Scalar, typename Urng, bool _mutable>
460  struct functor_traits<scalar_rng_adaptor<Gen, _Scalar, Urng, _mutable> >
461  {
462  enum { Cost = HugeCost, PacketAccess = packet_traits<_Scalar>::Vectorizable, IsRepeatable = false };
463  };
464  }
465 }
466 
467 #endif
PacketRandomEngine.h
Eigen::Rand::balanced
const BalancedType< Derived, Urng > balanced(Index rows, Index cols, Urng &&urng)
generates reals in a range [-1, 1]
Definition: Basic.h:397
MorePacketMath.h
PacketFilter.h