EigenRand  0.4.0-alpha
Discrete.h
Go to the documentation of this file.
1 
12 #ifndef EIGENRAND_DISTS_DISCRETE_H
13 #define EIGENRAND_DISTS_DISCRETE_H
14 
15 #include <memory>
16 #include <iterator>
17 #include <limits>
18 
19 namespace Eigen
20 {
21  namespace Rand
22  {
23  template<typename _Precision = uint32_t, typename _Size = uint32_t>
24  class AliasMethod
25  {
26  std::unique_ptr<_Precision[]> arr;
27  std::unique_ptr<_Size[]> alias;
28  size_t msize = 0, bitsize = 0, bitmask = 0;
29 
30  public:
31  AliasMethod()
32  {
33  }
34 
35  AliasMethod(const AliasMethod& o)
36  {
37  operator=(o);
38  }
39 
40  AliasMethod(AliasMethod&& o)
41  {
42  operator=(o);
43  }
44 
45  AliasMethod& operator=(const AliasMethod& o)
46  {
47  msize = o.msize;
48  bitsize = o.bitsize;
49  bitmask = o.bitmask;
50  if (msize)
51  {
52  arr = std::unique_ptr<_Precision[]>(new _Precision[(size_t)1 << bitsize]);
53  alias = std::unique_ptr<_Size[]>(new _Size[(size_t)1 << bitsize]);
54 
55  std::copy(o.arr.get(), o.arr.get() + ((size_t)1 << bitsize), arr.get());
56  std::copy(o.alias.get(), o.alias.get() + ((size_t)1 << bitsize), alias.get());
57  }
58  return *this;
59  }
60 
61  AliasMethod& operator=(AliasMethod&& o)
62  {
63  msize = o.msize;
64  bitsize = o.bitsize;
65  bitmask = o.bitmask;
66  std::swap(arr, o.arr);
67  std::swap(alias, o.alias);
68  return *this;
69  }
70 
71  template<typename _Iter>
72  AliasMethod(_Iter first, _Iter last)
73  {
74  buildTable(first, last);
75  }
76 
77  template<typename _Iter>
78  void buildTable(_Iter first, _Iter last)
79  {
80  size_t psize, nbsize;
81  msize = 0;
82  double sum = 0;
83  for (auto it = first; it != last; ++it, ++msize)
84  {
85  sum += *it;
86  }
87 
88  if (!std::isfinite(sum)) throw std::invalid_argument{ "cannot build NaN value distribution" };
89 
90  // ceil to power of 2
91  nbsize = (size_t)std::ceil(std::log2(msize));
92  psize = (size_t)1 << nbsize;
93 
94  if (nbsize != bitsize)
95  {
96  arr = std::unique_ptr<_Precision[]>(new _Precision[psize]);
97  std::fill(arr.get(), arr.get() + psize, 0);
98  alias = std::unique_ptr<_Size[]>(new _Size[psize]);
99  bitsize = nbsize;
100  bitmask = ((size_t)(-1)) >> (sizeof(size_t) * 8 - bitsize);
101  }
102 
103  sum /= psize;
104 
105  auto f = std::unique_ptr<double[]>(new double[psize]);
106  auto pf = f.get();
107  for (auto it = first; it != last; ++it, ++pf)
108  {
109  *pf = *it / sum;
110  }
111  std::fill(pf, pf + psize - msize, 0);
112 
113  size_t over = 0, under = 0, mm;
114  while (over < psize && f[over] < 1) ++over;
115  while (under < psize && f[under] >= 1) ++under;
116  mm = under + 1;
117 
118  while (over < psize && under < psize)
119  {
120  if (std::is_integral<_Precision>::value)
121  {
122  arr[under] = (_Precision)(f[under] * (std::numeric_limits<_Precision>::max() + 1.0));
123  }
124  else
125  {
126  arr[under] = (_Precision)f[under];
127  }
128  alias[under] = (_Size)over;
129  f[over] += f[under] - 1;
130  if (f[over] >= 1 || mm <= over)
131  {
132  for (under = mm; under < psize && f[under] >= 1; ++under);
133  mm = under + 1;
134  }
135  else
136  {
137  under = over;
138  }
139 
140  while (over < psize && f[over] < 1) ++over;
141  }
142 
143  for (; over < psize; ++over)
144  {
145  if (f[over] >= 1)
146  {
147  if (std::is_integral<_Precision>::value)
148  {
149  arr[over] = std::numeric_limits<_Precision>::max();
150  }
151  else
152  {
153  arr[over] = 1;
154  }
155  alias[over] = (_Size)over;
156  }
157  }
158 
159  if (under < psize)
160  {
161  if (std::is_integral<_Precision>::value)
162  {
163  arr[under] = std::numeric_limits<_Precision>::max();
164  }
165  else
166  {
167  arr[under] = 1;
168  }
169  alias[under] = (_Size)under;
170  for (under = mm; under < msize; ++under)
171  {
172  if (f[under] < 1)
173  {
174  if (std::is_integral<_Precision>::value)
175  {
176  arr[under] = std::numeric_limits<_Precision>::max();
177  }
178  else
179  {
180  arr[under] = 1;
181  }
182  alias[under] = (_Size)under;
183  }
184  }
185  }
186  }
187 
188  size_t get_bitsize() const
189  {
190  return bitsize;
191  }
192 
193  size_t get_bitmask() const
194  {
195  return bitmask;
196  }
197 
198  const _Precision* get_prob() const
199  {
200  return arr.get();
201  }
202 
203  const _Size* get_alias() const
204  {
205  return alias.get();
206  }
207  };
208 
214  template<typename _Scalar>
215  class UniformIntGen : OptCacheStore, public GenBase<UniformIntGen<_Scalar>, _Scalar>
216  {
217  static_assert(std::is_same<_Scalar, int32_t>::value, "uniformInt needs integral types.");
218  int cache_rest_cnt = 0;
219  RandbitsGen<_Scalar> randbits;
220  _Scalar pmin;
221  size_t pdiff, bitsize, bitmask;
222 
223  public:
224  using Scalar = _Scalar;
225 
231  UniformIntGen(_Scalar _min = 0, _Scalar _max = 0)
232  : pmin{ _min }, pdiff{ (size_t)(_max - _min) }
233  {
234  if ((pdiff + 1) > pdiff)
235  {
236  bitsize = (size_t)std::ceil(std::log2(pdiff + 1));
237  }
238  else
239  {
240  bitsize = (size_t)std::ceil(std::log2(pdiff));
241  }
242  bitmask = (_Scalar)(((size_t)-1) >> (sizeof(size_t) * 8 - bitsize));
243  }
244 
245  UniformIntGen(const UniformIntGen&) = default;
246  UniformIntGen(UniformIntGen&&) = default;
247 
248  UniformIntGen& operator=(const UniformIntGen&) = default;
249  UniformIntGen& operator=(UniformIntGen&&) = default;
250 
251  template<typename Rng>
252  EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng)
253  {
254  using namespace Eigen::internal;
255  auto rx = randbits(rng);
256  if (pdiff == bitmask)
257  {
258  return (_Scalar)(rx & bitmask) + pmin;
259  }
260  else
261  {
262  size_t bitcnt = bitsize;
263  for (int _i = 0; ; ++_i)
264  {
265  EIGENRAND_CHECK_INFINITY_LOOP();
266  _Scalar cands = (_Scalar)(rx & bitmask);
267  if (cands <= pdiff) return cands;
268  if (bitcnt + bitsize < 32)
269  {
270  rx >>= bitsize;
271  bitcnt += bitsize;
272  }
273  else
274  {
275  rx = randbits(rng);
276  bitcnt = bitsize;
277  }
278  }
279  }
280  }
281 
282  template<typename Packet, typename Rng>
283  EIGEN_STRONG_INLINE const Packet packetOp(Rng&& rng)
284  {
285  using namespace Eigen::internal;
286  auto rx = randbits.template packetOp<Packet>(rng);
287  auto pbitmask = pset1<Packet>(bitmask);
288  if (pdiff == bitmask)
289  {
290  return padd(pand(rx, pbitmask), pset1<Packet>(pmin));
291  }
292  else
293  {
294  auto& cm = Rand::detail::CompressMask<sizeof(Packet)>::get_inst();
295  auto plen = pset1<Packet>(pdiff + 1);
296  size_t bitcnt = bitsize;
297  for (int _i = 0; ; ++_i)
298  {
299  EIGENRAND_CHECK_INFINITY_LOOP();
300  // accept cands that only < plen
301  auto cands = pand(rx, pbitmask);
302  bool full = false;
303  cache_rest_cnt = cm.compress_append(cands, pcmplt(cands, plen),
304  OptCacheStore::template get<Packet>(), cache_rest_cnt, full);
305  if (full) return padd(cands, pset1<Packet>(pmin));
306 
307  if (bitcnt + bitsize < 32)
308  {
309  rx = psrl<-1>(rx, bitsize);
310  bitcnt += bitsize;
311  }
312  else
313  {
314  rx = randbits.template packetOp<Packet>(rng);
315  bitcnt = bitsize;
316  }
317  }
318  }
319  }
320  };
321 
328  template<typename _Scalar, typename Precision = float>
329  class DiscreteGen;
330 
336  template<typename _Scalar>
337  class DiscreteGen<_Scalar, int32_t> : public GenBase<DiscreteGen<_Scalar, int32_t>, _Scalar>
338  {
339  static_assert(std::is_same<_Scalar, int32_t>::value, "discreteDist needs integral types.");
340 #ifdef EIGEN_VECTORIZE_AVX2
341  OptCacheStore cache;
342  bool valid = false;
343 #endif
344  RandbitsGen<int32_t> randbits;
345  std::vector<uint32_t> cdf;
346  AliasMethod<int32_t, _Scalar> alias_table;
347 
348  public:
349  using Scalar = _Scalar;
350 
358  template<typename RealIter>
359  DiscreteGen(RealIter first, RealIter last)
360  {
361  if (std::distance(first, last) < 16)
362  {
363  // use linear or binary search
364  std::vector<double> _cdf;
365  double acc = 0;
366  for (; first != last; ++first)
367  {
368  _cdf.emplace_back(acc += *first);
369  }
370 
371  for (auto& p : _cdf)
372  {
373  cdf.emplace_back((uint32_t)(p / _cdf.back() * 0x80000000));
374  }
375  }
376  else
377  {
378  // use alias table
379  alias_table = AliasMethod<int32_t, _Scalar>{ first, last };
380  }
381  }
382 
390  template<typename Real,
391  typename std::enable_if<std::is_arithmetic<Real>::value, int>::type = 0>
392  DiscreteGen(const std::initializer_list<Real>& il)
393  : DiscreteGen(il.begin(), il.end())
394  {
395  }
396 
397  DiscreteGen()
398  : DiscreteGen({ 1 })
399  {
400  }
401 
402  DiscreteGen(const DiscreteGen&) = default;
403  DiscreteGen(DiscreteGen&&) = default;
404 
405  DiscreteGen& operator=(const DiscreteGen&) = default;
406  DiscreteGen& operator=(DiscreteGen&&) = default;
407 
408  template<typename Rng>
409  EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng)
410  {
411  using namespace Eigen::internal;
412  if (!cdf.empty())
413  {
414  auto rx = randbits(std::forward<Rng>(rng)) & 0x7FFFFFFF;
415  return (_Scalar)(std::lower_bound(cdf.begin(), cdf.end() - 1, rx) - cdf.begin());
416  }
417  else
418  {
419  auto rx = randbits(std::forward<Rng>(rng));
420  auto albit = rx & alias_table.get_bitmask();
421  uint32_t alx = (uint32_t)(rx >> (sizeof(rx) * 8 - 31));
422  if (alx < alias_table.get_prob()[albit]) return (_Scalar)albit;
423  return alias_table.get_alias()[albit];
424  }
425  }
426 
427  template<typename Packet, typename Rng>
428  EIGEN_STRONG_INLINE const Packet packetOp(Rng&& rng)
429  {
430  using namespace Eigen::internal;
431 #ifdef EIGEN_VECTORIZE_AVX2
432  if (valid)
433  {
434  valid = false;
435  return cache.template get<Packet>();
436  }
437  using PacketType = Packet8i;
438 #else
439  using PacketType = Packet;
440 #endif
441 
442  PacketType ret;
443  if (!cdf.empty())
444  {
445  ret = pset1<PacketType>(cdf.size() - 1);
446  auto rx = pand(randbits.template packetOp<PacketType>(std::forward<Rng>(rng)), pset1<PacketType>(0x7FFFFFFF));
447  for (size_t i = 0; i < cdf.size() - 1; ++i)
448  {
449  ret = padd(ret, pcmplt(rx, pset1<PacketType>(cdf[i])));
450  }
451  }
452  else
453  {
454  auto rx = randbits.template packetOp<PacketType>(std::forward<Rng>(rng));
455  auto albit = pand(rx, pset1<PacketType>(alias_table.get_bitmask()));
456  auto c = pcmplt(psrl<1>(rx), pgather(alias_table.get_prob(), albit));
457  ret = pblendv(c, albit, pgather(alias_table.get_alias(), albit));
458  }
459 
460 #ifdef EIGEN_VECTORIZE_AVX2
461  valid = true;
462  cache.template get<Packet>() = _mm256_extractf128_si256(ret, 1);
463  return _mm256_extractf128_si256(ret, 0);
464 #else
465  return ret;
466 #endif
467  }
468  };
469 
475  template<typename _Scalar>
476  class DiscreteGen<_Scalar, float> : public GenBase<DiscreteGen<_Scalar, float>, _Scalar>
477  {
478  static_assert(std::is_same<_Scalar, int32_t>::value, "discreteDist needs integral types.");
479  UniformRealGen<float> ur;
480  std::vector<float> cdf;
481  AliasMethod<float, _Scalar> alias_table;
482 
483  public:
484  using Scalar = _Scalar;
485 
493  template<typename RealIter>
494  DiscreteGen(RealIter first, RealIter last)
495  {
496  if (std::distance(first, last) < 16)
497  {
498  // use linear or binary search
499  std::vector<double> _cdf;
500  double acc = 0;
501  for (; first != last; ++first)
502  {
503  _cdf.emplace_back(acc += *first);
504  }
505 
506  for (auto& p : _cdf)
507  {
508  cdf.emplace_back(p / _cdf.back());
509  }
510  }
511  else
512  {
513  // use alias table
514  alias_table = AliasMethod<float, _Scalar>{ first, last };
515  }
516  }
517 
525  template<typename Real,
526  typename std::enable_if<std::is_arithmetic<Real>::value, int>::type = 0>
527  DiscreteGen(const std::initializer_list<Real>& il)
528  : DiscreteGen(il.begin(), il.end())
529  {
530  }
531 
532  DiscreteGen()
533  : DiscreteGen({ 1 })
534  {
535  }
536 
537  DiscreteGen(const DiscreteGen&) = default;
538  DiscreteGen(DiscreteGen&&) = default;
539 
540  DiscreteGen& operator=(const DiscreteGen&) = default;
541  DiscreteGen& operator=(DiscreteGen&&) = default;
542 
543  template<typename Rng>
544  EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng)
545  {
546  using namespace Eigen::internal;
547  if (!cdf.empty())
548  {
549  auto rx = ur(std::forward<Rng>(rng));
550  return (_Scalar)(std::lower_bound(cdf.begin(), cdf.end() - 1, rx) - cdf.begin());
551  }
552  else
553  {
554  auto albit = pfirst(rng()) & alias_table.get_bitmask();
555  auto alx = ur(rng);
556  if (alx < alias_table.get_prob()[albit]) return albit;
557  return alias_table.get_alias()[albit];
558  }
559  }
560 
561  template<typename Packet, typename Rng>
562  EIGEN_STRONG_INLINE const Packet packetOp(Rng&& rng)
563  {
564  using namespace Eigen::internal;
565  using PacketType = decltype(reinterpret_to_float(std::declval<Packet>()));
566 
567  if (!cdf.empty())
568  {
569  auto ret = pset1<Packet>(cdf.size());
570  auto rx = ur.template packetOp<PacketType>(std::forward<Rng>(rng));
571  for (auto& p : cdf)
572  {
573  ret = padd(ret, reinterpret_to_int(pcmplt(rx, pset1<PacketType>(p))));
574  }
575  return ret;
576  }
577  else
578  {
579  using RUtils = RawbitsMaker<Packet, Rng>;
580  auto albit = pand(RUtils{}.rawbits(rng), pset1<Packet>(alias_table.get_bitmask()));
581  auto c = reinterpret_to_int(pcmplt(ur.template packetOp<PacketType>(rng), pgather(alias_table.get_prob(), albit)));
582  return pblendv(c, albit, pgather(alias_table.get_alias(), albit));
583  }
584  }
585  };
586 
592  template<typename _Scalar>
593  class DiscreteGen<_Scalar, double> : public GenBase<DiscreteGen<_Scalar, double>, _Scalar>
594  {
595  static_assert(std::is_same<_Scalar, int32_t>::value, "discreteDist needs integral types.");
596  UniformRealGen<double> ur;
597  std::vector<double> cdf;
598  AliasMethod<double, _Scalar> alias_table;
599 
600  public:
601  using Scalar = _Scalar;
602 
610  template<typename RealIter>
611  DiscreteGen(RealIter first, RealIter last)
612  {
613  if (std::distance(first, last) < 16)
614  {
615  // use linear or binary search
616  std::vector<double> _cdf;
617  double acc = 0;
618  for (; first != last; ++first)
619  {
620  _cdf.emplace_back(acc += *first);
621  }
622 
623  for (auto& p : _cdf)
624  {
625  cdf.emplace_back(p / _cdf.back());
626  }
627  }
628  else
629  {
630  // use alias table
631  alias_table = AliasMethod<double, _Scalar>{ first, last };
632  }
633  }
634 
642  template<typename Real,
643  typename std::enable_if<std::is_arithmetic<Real>::value, int>::type = 0>
644  DiscreteGen(const std::initializer_list<Real>& il)
645  : DiscreteGen(il.begin(), il.end())
646  {
647  }
648 
649  DiscreteGen()
650  : DiscreteGen({ 1 })
651  {
652  }
653 
654  DiscreteGen(const DiscreteGen&) = default;
655  DiscreteGen(DiscreteGen&&) = default;
656 
657  DiscreteGen& operator=(const DiscreteGen&) = default;
658  DiscreteGen& operator=(DiscreteGen&&) = default;
659 
660  template<typename Rng>
661  EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng)
662  {
663  using namespace Eigen::internal;
664  if (!cdf.empty())
665  {
666  auto rx = ur(std::forward<Rng>(rng));
667  return (_Scalar)(std::lower_bound(cdf.begin(), cdf.end() - 1, rx) - cdf.begin());
668  }
669  else
670  {
671  auto albit = pfirst(rng()) & alias_table.get_bitmask();
672  auto alx = ur(rng);
673  if (alx < alias_table.get_prob()[albit]) return albit;
674  return alias_table.get_alias()[albit];
675  }
676  }
677 
678 #ifdef EIGEN_VECTORIZE_NEON
679 #else
680  template<typename Packet, typename Rng>
681  EIGEN_STRONG_INLINE const Packet packetOp(Rng&& rng)
682  {
683  using namespace Eigen::internal;
684  using DPacket = decltype(reinterpret_to_double(std::declval<Packet>()));
685  if (!cdf.empty())
686  {
687  auto ret = pset1<Packet>(cdf.size());
688  #ifdef EIGEN_VECTORIZE_AVX
689  auto rx = ur.template packetOp<Packet4d>(std::forward<Rng>(rng));
690  for (auto& p : cdf)
691  {
692  auto c = reinterpret_to_int(pcmplt(rx, pset1<decltype(rx)>(p)));
693  auto r = combine_low32(c);
694  ret = padd(ret, r);
695  }
696  #else
697  auto rx1 = ur.template packetOp<DPacket>(rng),
698  rx2 = ur.template packetOp<DPacket>(rng);
699  for (auto& p : cdf)
700  {
701  auto pp = pset1<decltype(rx1)>(p);
702  ret = padd(ret, combine_low32(reinterpret_to_int(pcmplt(rx1, pp)), reinterpret_to_int(pcmplt(rx2, pp))));
703  }
704  #endif
705  return ret;
706  }
707  else
708  {
709  #ifdef EIGEN_VECTORIZE_AVX
710  using RUtils = RawbitsMaker<Packet, Rng>;
711  auto albit = pand(RUtils{}.rawbits(rng), pset1<Packet>(alias_table.get_bitmask()));
712  auto c = reinterpret_to_int(pcmplt(ur.template packetOp<Packet4d>(rng), pgather(alias_table.get_prob(), _mm256_castsi128_si256(albit))));
713  return pblendv(combine_low32(c), albit, pgather(alias_table.get_alias(), albit));
714  #else
715  using RUtils = RawbitsMaker<Packet, Rng>;
716  auto albit = pand(RUtils{}.rawbits(rng), pset1<Packet>(alias_table.get_bitmask()));
717  auto c1 = reinterpret_to_int(pcmplt(ur.template packetOp<DPacket>(rng), pgather(alias_table.get_prob(), albit)));
718  auto c2 = reinterpret_to_int(pcmplt(ur.template packetOp<DPacket>(rng), pgather(alias_table.get_prob(), albit, true)));
719  return pblendv(combine_low32(c1, c2), albit, pgather(alias_table.get_alias(), albit));
720  #endif
721  }
722  }
723 #endif
724  };
725 
726  template<typename> class BinomialGen;
727 
733  template<typename _Scalar>
734  class PoissonGen : OptCacheStore, public GenBase<PoissonGen<_Scalar>, _Scalar>
735  {
736  friend BinomialGen<_Scalar>;
737  static_assert(std::is_same<_Scalar, int32_t>::value, "poisson needs integral types.");
738  int cache_rest_cnt = 0;
739  UniformRealGen<float> ur;
740 
741  protected:
742  double mean, ne_mean, sqrt_tmean, log_mean, g1;
743 
744  public:
745  using Scalar = _Scalar;
746 
752  PoissonGen(double _mean = 1)
753  : mean{ _mean }, ne_mean{ std::exp(-_mean) }
754  {
755  sqrt_tmean = std::sqrt(2 * mean);
756  log_mean = std::log(mean);
757  g1 = mean * log_mean - std::lgamma(mean + 1);
758  }
759 
760  PoissonGen(const PoissonGen&) = default;
761  PoissonGen(PoissonGen&&) = default;
762 
763  PoissonGen& operator=(const PoissonGen&) = default;
764  PoissonGen& operator=(PoissonGen&&) = default;
765 
766  template<typename Rng>
767  EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng)
768  {
769  using namespace Eigen::internal;
770  if (mean < 12)
771  {
772  _Scalar res = 0;
773  double val = 1;
774  for (; ; ++res)
775  {
776  val *= ur(rng);
777  if (val <= ne_mean) break;
778  }
779  return res;
780  }
781  else
782  {
783  _Scalar res;
784  double yx;
785  for (int _i = 0; ; ++_i)
786  {
787  EIGENRAND_CHECK_INFINITY_LOOP();
788  yx = std::tan(constant::pi * ur(rng));
789  res = (_Scalar)(sqrt_tmean * yx + mean);
790  if (res >= 0 && ur(rng) <= 0.9 * (1.0 + yx * yx)
791  * std::exp(res * log_mean - std::lgamma(res + 1.0) - g1))
792  {
793  return res;
794  }
795  }
796  }
797  }
798 
799  template<typename Packet, typename Rng>
800  EIGEN_STRONG_INLINE const Packet packetOp(Rng&& rng)
801  {
802  using namespace Eigen::internal;
803  using PacketType = decltype(reinterpret_to_float(std::declval<Packet>()));
804 
805  if (mean < 12)
806  {
807  Packet res = pset1<Packet>(0);
808  PacketType val = pset1<PacketType>(1), pne_mean = pset1<PacketType>(ne_mean);
809  for (int _i = 0; ; ++_i)
810  {
811  EIGENRAND_CHECK_INFINITY_LOOP();
812  val = pmul(val, ur.template packetOp<PacketType>(rng));
813  auto c = reinterpret_to_int(pcmplt(pne_mean, val));
814  if (pmovemask(c) == 0) break;
815  res = padd(res, pnegate(c));
816  }
817  return res;
818  }
819  else
820  {
821  auto& cm = Rand::detail::CompressMask<sizeof(Packet)>::get_inst();
822  const PacketType ppi = pset1<PacketType>(constant::pi),
823  psqrt_tmean = pset1<PacketType>(sqrt_tmean),
824  pmean = pset1<PacketType>(mean),
825  plog_mean = pset1<PacketType>(log_mean),
826  pg1 = pset1<PacketType>(g1);
827  for (int _i = 0; ; ++_i)
828  {
829  EIGENRAND_CHECK_INFINITY_LOOP();
830  PacketType fres, yx, psin, pcos;
831  psincos(pmul(ppi, ur.template packetOp<PacketType>(rng)), psin, pcos);
832  yx = pdiv(psin, pcos);
833  fres = ptruncate(padd(pmul(psqrt_tmean, yx), pmean));
834 
835  auto p1 = pmul(padd(pmul(yx, yx), pset1<PacketType>(1)), pset1<PacketType>(0.9));
836  auto p2 = pexp(psub(psub(pmul(fres, plog_mean), plgamma_approx(padd(fres, pset1<PacketType>(1)))), pg1));
837 
838  auto c1 = pcmple(pset1<PacketType>(0), fres);
839  auto c2 = pcmple(ur.template packetOp<PacketType>(rng), pmul(p1, p2));
840 
841  auto cands = fres;
842  bool full = false;
843  cache_rest_cnt = cm.compress_append(cands, pand(c1, c2),
844  OptCacheStore::template get<PacketType>(), cache_rest_cnt, full);
845  if (full) return pcast<PacketType, Packet>(cands);
846  }
847  }
848  }
849  };
850 
856  template<typename _Scalar>
857  class BinomialGen : public GenBase<BinomialGen<_Scalar>, _Scalar>
858  {
859  static_assert(std::is_same<_Scalar, int32_t>::value, "binomial needs integral types.");
860 
861  PoissonGen<_Scalar> poisson;
862  _Scalar trials;
863  double p = 0, small_p = 0, g1 = 0, sqrt_v = 0, log_small_p = 0, log_small_q = 0;
864 
865  public:
866  using Scalar = _Scalar;
867 
874  BinomialGen(_Scalar _trials = 1, double _p = 0.5)
875  : poisson{ _trials * std::min(_p, 1 - _p) },
876  trials{ _trials }, p{ _p }, small_p{ std::min(p, 1 - p) }
877 
878  {
879  if (!(trials < 25 || poisson.mean < 1.0))
880  {
881  g1 = std::lgamma(trials + 1);
882  sqrt_v = std::sqrt(2 * poisson.mean * (1 - small_p));
883  log_small_p = std::log(small_p);
884  log_small_q = std::log(1 - small_p);
885  }
886  }
887 
888  BinomialGen(const BinomialGen&) = default;
889  BinomialGen(BinomialGen&&) = default;
890 
891  BinomialGen& operator=(const BinomialGen&) = default;
892  BinomialGen& operator=(BinomialGen&&) = default;
893 
894  template<typename Rng>
895  EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng)
896  {
897  using namespace Eigen::internal;
898  _Scalar res;
899  if (trials < 25)
900  {
901  res = 0;
902  for (int i = 0; i < trials; ++i)
903  {
904  if (poisson.ur(rng) < p) ++res;
905  }
906  return res;
907  }
908  else if (poisson.mean < 1.0)
909  {
910  res = poisson(rng);
911  }
912  else
913  {
914  for (int _i = 0; ; ++_i)
915  {
916  EIGENRAND_CHECK_INFINITY_LOOP();
917  double ys;
918  ys = std::tan(constant::pi * poisson.ur(rng));
919  res = (_Scalar)(sqrt_v * ys + poisson.mean);
920  if (0 <= res && res <= trials && poisson.ur(rng) <= 1.2 * sqrt_v
921  * (1.0 + ys * ys)
922  * std::exp(g1 - std::lgamma(res + 1)
923  - std::lgamma(trials - res + 1.0)
924  + res * log_small_p
925  + (trials - res) * log_small_q)
926  )
927  {
928  break;
929  }
930  }
931  }
932  return p == small_p ? res : trials - res;
933  }
934 
935  template<typename Packet, typename Rng>
936  EIGEN_STRONG_INLINE const Packet packetOp(Rng&& rng)
937  {
938  using namespace Eigen::internal;
939  using PacketType = decltype(reinterpret_to_float(std::declval<Packet>()));
940  Packet res;
941  if (trials < 25)
942  {
943  PacketType pp = pset1<PacketType>(p);
944  res = pset1<Packet>(trials);
945  for (int i = 0; i < trials; ++i)
946  {
947  auto c = reinterpret_to_int(pcmple(pp, poisson.ur.template packetOp<PacketType>(rng)));
948  res = padd(res, c);
949  }
950  return res;
951  }
952  else if (poisson.mean < 1.0)
953  {
954  res = poisson.template packetOp<Packet>(rng);
955  }
956  else
957  {
958  auto& cm = Rand::detail::CompressMask<sizeof(Packet)>::get_inst();
959  const PacketType ppi = pset1<PacketType>(constant::pi),
960  ptrials = pset1<PacketType>(trials),
961  psqrt_v = pset1<PacketType>(sqrt_v),
962  pmean = pset1<PacketType>(poisson.mean),
963  plog_small_p = pset1<PacketType>(log_small_p),
964  plog_small_q = pset1<PacketType>(log_small_q),
965  pg1 = pset1<PacketType>(g1);
966  for (int _i = 0; ; ++_i)
967  {
968  EIGENRAND_CHECK_INFINITY_LOOP();
969  PacketType fres, ys, psin, pcos;
970  psincos(pmul(ppi, poisson.ur.template packetOp<PacketType>(rng)), psin, pcos);
971  ys = pdiv(psin, pcos);
972  fres = ptruncate(padd(pmul(psqrt_v, ys), pmean));
973 
974  auto p1 = pmul(pmul(pset1<PacketType>(1.2), psqrt_v), padd(pset1<PacketType>(1), pmul(ys, ys)));
975  auto p2 = pexp(
976  padd(padd(psub(
977  psub(pg1, plgamma_approx(padd(fres, pset1<PacketType>(1)))),
978  plgamma_approx(psub(padd(ptrials, pset1<PacketType>(1)), fres))
979  ), pmul(fres, plog_small_p)), pmul(psub(ptrials, fres), plog_small_q))
980  );
981 
982  auto c1 = pand(pcmple(pset1<PacketType>(0), fres), pcmple(fres, ptrials));
983  auto c2 = pcmple(poisson.ur.template packetOp<PacketType>(rng), pmul(p1, p2));
984 
985  auto cands = fres;
986  bool full = false;
987  poisson.cache_rest_cnt = cm.compress_append(cands, pand(c1, c2),
988  poisson.template get<PacketType>(), poisson.cache_rest_cnt, full);
989  if (full)
990  {
991  res = pcast<PacketType, Packet>(cands);
992  break;
993  }
994  }
995  }
996  return p == small_p ? res : psub(pset1<Packet>(trials), res);
997  }
998  };
999 
1005  template<typename _Scalar>
1006  class GeometricGen : public GenBase<GeometricGen<_Scalar>, _Scalar>
1007  {
1008  static_assert(std::is_same<_Scalar, int32_t>::value, "geomtric needs integral types.");
1009  UniformRealGen<float> ur;
1010  double p, rlog_q;
1011 
1012  public:
1013  using Scalar = _Scalar;
1014 
1020  GeometricGen(double _p = 0.5)
1021  : p{ _p }, rlog_q{ 1 / std::log(1 - p) }
1022  {
1023  }
1024 
1025  GeometricGen(const GeometricGen&) = default;
1026  GeometricGen(GeometricGen&&) = default;
1027 
1028  GeometricGen& operator=(const GeometricGen&) = default;
1029  GeometricGen& operator=(GeometricGen&&) = default;
1030 
1031  template<typename Rng>
1032  EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng)
1033  {
1034  using namespace Eigen::internal;
1035  return (_Scalar)(std::log(1 - ur(std::forward<Rng>(rng))) * rlog_q);
1036  }
1037 
1038  template<typename Packet, typename Rng>
1039  EIGEN_STRONG_INLINE const Packet packetOp(Rng&& rng)
1040  {
1041  using namespace Eigen::internal;
1042  using PacketType = decltype(reinterpret_to_float(std::declval<Packet>()));
1043 
1044  return pcast<PacketType, Packet>(ptruncate(pmul(plog(
1045  psub(pset1<PacketType>(1), ur.template packetOp<PacketType>(std::forward<Rng>(rng)))
1046  ), pset1<PacketType>(rlog_q))));
1047  }
1048  };
1049 
1050 
1051  template<typename Derived, typename Urng>
1052  using UniformIntType = CwiseNullaryOp<internal::scalar_rng_adaptor<UniformIntGen<typename Derived::Scalar>, typename Derived::Scalar, Urng, true>, const Derived>;
1053 
1067  template<typename Derived, typename Urng>
1068  inline const UniformIntType<Derived, Urng>
1069  uniformInt(Index rows, Index cols, Urng&& urng, typename Derived::Scalar min, typename Derived::Scalar max)
1070  {
1071  return {
1072  rows, cols, { std::forward<Urng>(urng), UniformIntGen<typename Derived::Scalar>{ min, max } }
1073  };
1074  }
1075 
1088  template<typename Derived, typename Urng>
1089  inline const UniformIntType<Derived, Urng>
1090  uniformIntLike(Derived& o, Urng&& urng, typename Derived::Scalar min, typename Derived::Scalar max)
1091  {
1092  return {
1093  o.rows(), o.cols(), { std::forward<Urng>(urng), UniformIntGen<typename Derived::Scalar>{ min, max } }
1094  };
1095  }
1096 
1097  template<typename Derived, typename Urng>
1098  using DiscreteFType = CwiseNullaryOp<internal::scalar_rng_adaptor<DiscreteGen<typename Derived::Scalar, float>, typename Derived::Scalar, Urng, true>, const Derived>;
1099 
1114  template<typename Derived, typename Urng, typename RealIter>
1115  inline const DiscreteFType<Derived, Urng>
1116  discreteF(Index rows, Index cols, Urng&& urng, RealIter first, RealIter last)
1117  {
1118  return {
1119  rows, cols, { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, float>{first, last} }
1120  };
1121  }
1122 
1136  template<typename Derived, typename Urng, typename RealIter>
1137  inline const DiscreteFType<Derived, Urng>
1138  discreteFLike(Derived& o, Urng&& urng, RealIter first, RealIter last)
1139  {
1140  return {
1141  o.rows(), o.cols(), { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, float>(first, last) }
1142  };
1143  }
1144 
1159  template<typename Derived, typename Urng, typename Real>
1160  inline const DiscreteFType<Derived, Urng>
1161  discreteF(Index rows, Index cols, Urng&& urng, const std::initializer_list<Real>& il)
1162  {
1163  return {
1164  rows, cols, { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, float>{il.begin(), il.end()} }
1165  };
1166  }
1167 
1181  template<typename Derived, typename Urng, typename Real>
1182  inline const DiscreteFType<Derived, Urng>
1183  discreteFLike(Derived& o, Urng&& urng, const std::initializer_list<Real>& il)
1184  {
1185  return {
1186  o.rows(), o.cols(), { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, float>(il.begin(), il.end()) }
1187  };
1188  }
1189 
1190  template<typename Derived, typename Urng>
1191  using DiscreteDType = CwiseNullaryOp<internal::scalar_rng_adaptor<DiscreteGen<typename Derived::Scalar, double>, typename Derived::Scalar, Urng, true>, const Derived>;
1192 
1207  template<typename Derived, typename Urng, typename RealIter>
1208  inline const DiscreteDType<Derived, Urng>
1209  discreteD(Index rows, Index cols, Urng&& urng, RealIter first, RealIter last)
1210  {
1211  return {
1212  rows, cols, { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, double>{first, last} }
1213  };
1214  }
1215 
1229  template<typename Derived, typename Urng, typename RealIter>
1230  inline const DiscreteDType<Derived, Urng>
1231  discreteDLike(Derived& o, Urng&& urng, RealIter first, RealIter last)
1232  {
1233  return {
1234  o.rows(), o.cols(), { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, double>{first, last} }
1235  };
1236  }
1237 
1252  template<typename Derived, typename Urng, typename Real>
1253  inline const DiscreteDType<Derived, Urng>
1254  discreteD(Index rows, Index cols, Urng&& urng, const std::initializer_list<Real>& il)
1255  {
1256  return {
1257  rows, cols, { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, double>{il.begin(), il.end()} }
1258  };
1259  }
1260 
1274  template<typename Derived, typename Urng, typename Real>
1275  inline const DiscreteDType<Derived, Urng>
1276  discreteDLike(Derived& o, Urng&& urng, const std::initializer_list<Real>& il)
1277  {
1278  return {
1279  o.rows(), o.cols(), { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, double>{il.begin(), il.end()} }
1280  };
1281  }
1282 
1283  template<typename Derived, typename Urng>
1284  using DiscreteType = CwiseNullaryOp<internal::scalar_rng_adaptor<DiscreteGen<typename Derived::Scalar, int32_t>, typename Derived::Scalar, Urng, true>, const Derived>;
1285 
1300  template<typename Derived, typename Urng, typename RealIter>
1301  inline const DiscreteType<Derived, Urng>
1302  discrete(Index rows, Index cols, Urng&& urng, RealIter first, RealIter last)
1303  {
1304  return {
1305  rows, cols, { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, int32_t>{first, last} }
1306  };
1307  }
1308 
1322  template<typename Derived, typename Urng, typename RealIter>
1323  inline const DiscreteType<Derived, Urng>
1324  discreteLike(Derived& o, Urng&& urng, RealIter first, RealIter last)
1325  {
1326  return {
1327  o.rows(), o.cols(), { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, int32_t>{first, last} }
1328  };
1329  }
1330 
1345  template<typename Derived, typename Urng, typename Real>
1346  inline const DiscreteType<Derived, Urng>
1347  discrete(Index rows, Index cols, Urng&& urng, const std::initializer_list<Real>& il)
1348  {
1349  return {
1350  rows, cols, { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, int32_t>{il.begin(), il.end()} }
1351  };
1352  }
1353 
1367  template<typename Derived, typename Urng, typename Real>
1368  inline const DiscreteType<Derived, Urng>
1369  discreteLike(Derived& o, Urng&& urng, const std::initializer_list<Real>& il)
1370  {
1371  return {
1372  o.rows(), o.cols(), { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, int32_t>{il.begin(), il.end()} }
1373  };
1374  }
1375 
1376  template<typename Derived, typename Urng>
1377  using PoissonType = CwiseNullaryOp<internal::scalar_rng_adaptor<PoissonGen<typename Derived::Scalar>, typename Derived::Scalar, Urng, true>, const Derived>;
1378 
1392  template<typename Derived, typename Urng>
1393  inline const PoissonType<Derived, Urng>
1394  poisson(Index rows, Index cols, Urng&& urng, double mean = 1)
1395  {
1396  return {
1397  rows, cols, { std::forward<Urng>(urng), PoissonGen<typename Derived::Scalar>{mean} }
1398  };
1399  }
1400 
1413  template<typename Derived, typename Urng>
1414  inline const PoissonType<Derived, Urng>
1415  poissonLike(Derived& o, Urng&& urng, double mean = 1)
1416  {
1417  return {
1418  o.rows(), o.cols(), { std::forward<Urng>(urng), PoissonGen<typename Derived::Scalar>{mean} }
1419  };
1420  }
1421 
1422  template<typename Derived, typename Urng>
1423  using BinomialType = CwiseNullaryOp<internal::scalar_rng_adaptor<BinomialGen<typename Derived::Scalar>, typename Derived::Scalar, Urng, true>, const Derived>;
1424 
1439  template<typename Derived, typename Urng>
1440  inline const BinomialType<Derived, Urng>
1441  binomial(Index rows, Index cols, Urng&& urng, typename Derived::Scalar trials = 1, double p = 0.5)
1442  {
1443  return {
1444  rows, cols, { std::forward<Urng>(urng), BinomialGen<typename Derived::Scalar>{trials, p} }
1445  };
1446  }
1447 
1461  template<typename Derived, typename Urng>
1462  inline const BinomialType<Derived, Urng>
1463  binomialLike(Derived& o, Urng&& urng, typename Derived::Scalar trials = 1, double p = 0.5)
1464  {
1465  return {
1466  o.rows(), o.cols(), { std::forward<Urng>(urng), BinomialGen<typename Derived::Scalar>{trials, p} }
1467  };
1468  }
1469 
1470  template<typename Derived, typename Urng>
1471  using GeometricType = CwiseNullaryOp<internal::scalar_rng_adaptor<GeometricGen<typename Derived::Scalar>, typename Derived::Scalar, Urng, true>, const Derived>;
1472 
1486  template<typename Derived, typename Urng>
1487  inline const GeometricType<Derived, Urng>
1488  geometric(Index rows, Index cols, Urng&& urng, double p = 0.5)
1489  {
1490  return {
1491  rows, cols, { std::forward<Urng>(urng), GeometricGen<typename Derived::Scalar>{p} }
1492  };
1493  }
1494 
1507  template<typename Derived, typename Urng>
1508  inline const GeometricType<Derived, Urng>
1509  geometricLike(Derived& o, Urng&& urng, double p = 0.5)
1510  {
1511  return {
1512  o.rows(), o.cols(), { std::forward<Urng>(urng), GeometricGen<typename Derived::Scalar>{p} }
1513  };
1514  }
1515  }
1516 
1517 #ifdef EIGEN_VECTORIZE_NEON
1518  namespace internal
1519  {
1520  template<typename _Scalar, typename Urng, bool _mutable>
1521  struct functor_traits<scalar_rng_adaptor<Rand::DiscreteGen<_Scalar, double>, _Scalar, Urng, _mutable> >
1522  {
1523  enum { Cost = HugeCost, PacketAccess = 0, IsRepeatable = false };
1524  };
1525  }
1526 #endif
1527 }
1528 
1529 #endif
Eigen::Rand::BinomialGen
Generator of integers on a binomial distribution.
Definition: Discrete.h:858
Eigen::Rand::BinomialGen::BinomialGen
BinomialGen(_Scalar _trials=1, double _p=0.5)
Construct a new Binomial Generator.
Definition: Discrete.h:874
Eigen::Rand::DiscreteGen< _Scalar, double >::DiscreteGen
DiscreteGen(RealIter first, RealIter last)
Construct a new Discrete Generator.
Definition: Discrete.h:611
Eigen::Rand::DiscreteGen< _Scalar, double >::DiscreteGen
DiscreteGen(const std::initializer_list< Real > &il)
Construct a new Discrete Generator.
Definition: Discrete.h:644
Eigen::Rand::DiscreteGen< _Scalar, float >::DiscreteGen
DiscreteGen(RealIter first, RealIter last)
Construct a new Discrete Generator.
Definition: Discrete.h:494
Eigen::Rand::DiscreteGen< _Scalar, float >::DiscreteGen
DiscreteGen(const std::initializer_list< Real > &il)
Construct a new Discrete Generator.
Definition: Discrete.h:527
Eigen::Rand::DiscreteGen< _Scalar, int32_t >::DiscreteGen
DiscreteGen(RealIter first, RealIter last)
Construct a new Discrete Generator.
Definition: Discrete.h:359
Eigen::Rand::DiscreteGen< _Scalar, int32_t >::DiscreteGen
DiscreteGen(const std::initializer_list< Real > &il)
Construct a new Discrete Generator.
Definition: Discrete.h:392
Eigen::Rand::DiscreteGen
Generator of integers on the interval [0, n), where the probability of each individual integer i is p...
Definition: Discrete.h:329
Eigen::Rand::GenBase
Base class of all univariate random generators.
Definition: Basic.h:33
Eigen::Rand::GeometricGen
Generator of integers on a geometric distribution.
Definition: Discrete.h:1007
Eigen::Rand::GeometricGen::GeometricGen
GeometricGen(double _p=0.5)
Construct a new Geometric Generator.
Definition: Discrete.h:1020
Eigen::Rand::PoissonGen
Generator of integers on a Poisson distribution.
Definition: Discrete.h:735
Eigen::Rand::PoissonGen::PoissonGen
PoissonGen(double _mean=1)
Construct a new Poisson Generator.
Definition: Discrete.h:752
Eigen::Rand::RandbitsGen
Generator of random bits for integral scalars.
Definition: Basic.h:273
Eigen::Rand::UniformIntGen
Generator of integers with a given range [min, max]
Definition: Discrete.h:216
Eigen::Rand::UniformIntGen::UniformIntGen
UniformIntGen(_Scalar _min=0, _Scalar _max=0)
Construct a new UniformInt Generator.
Definition: Discrete.h:231
Eigen::Rand::poissonLike
const PoissonType< Derived, Urng > poissonLike(Derived &o, Urng &&urng, double mean=1)
generates reals on the Poisson distribution.
Definition: Discrete.h:1415
Eigen::Rand::discreteF
const DiscreteFType< Derived, Urng > discreteF(Index rows, Index cols, Urng &&urng, RealIter first, RealIter last)
generates random integers on the interval [0, n), where the probability of each individual integer i ...
Definition: Discrete.h:1116
Eigen::Rand::geometricLike
const GeometricType< Derived, Urng > geometricLike(Derived &o, Urng &&urng, double p=0.5)
generates reals on the geometric distribution.
Definition: Discrete.h:1509
Eigen::Rand::binomialLike
const BinomialType< Derived, Urng > binomialLike(Derived &o, Urng &&urng, typename Derived::Scalar trials=1, double p=0.5)
generates reals on the binomial distribution.
Definition: Discrete.h:1463
Eigen::Rand::uniformInt
const UniformIntType< Derived, Urng > uniformInt(Index rows, Index cols, Urng &&urng, typename Derived::Scalar min, typename Derived::Scalar max)
generates integers with a given range [min, max]
Definition: Discrete.h:1069
Eigen::Rand::discreteDLike
const DiscreteDType< Derived, Urng > discreteDLike(Derived &o, Urng &&urng, RealIter first, RealIter last)
generates random integers on the interval [0, n), where the probability of each individual integer i ...
Definition: Discrete.h:1231
Eigen::Rand::discreteFLike
const DiscreteFType< Derived, Urng > discreteFLike(Derived &o, Urng &&urng, RealIter first, RealIter last)
generates random integers on the interval [0, n), where the probability of each individual integer i ...
Definition: Discrete.h:1138
Eigen::Rand::uniformIntLike
const UniformIntType< Derived, Urng > uniformIntLike(Derived &o, Urng &&urng, typename Derived::Scalar min, typename Derived::Scalar max)
generates integers with a given range [min, max]
Definition: Discrete.h:1090
Eigen::Rand::discreteD
const DiscreteDType< Derived, Urng > discreteD(Index rows, Index cols, Urng &&urng, RealIter first, RealIter last)
generates random integers on the interval [0, n), where the probability of each individual integer i ...
Definition: Discrete.h:1209
Eigen::Rand::poisson
const PoissonType< Derived, Urng > poisson(Index rows, Index cols, Urng &&urng, double mean=1)
generates reals on the Poisson distribution.
Definition: Discrete.h:1394
Eigen::Rand::discrete
const DiscreteType< Derived, Urng > discrete(Index rows, Index cols, Urng &&urng, RealIter first, RealIter last)
generates random integers on the interval [0, n), where the probability of each individual integer i ...
Definition: Discrete.h:1302
Eigen::Rand::geometric
const GeometricType< Derived, Urng > geometric(Index rows, Index cols, Urng &&urng, double p=0.5)
generates reals on the geometric distribution.
Definition: Discrete.h:1488
Eigen::Rand::discreteLike
const DiscreteType< Derived, Urng > discreteLike(Derived &o, Urng &&urng, RealIter first, RealIter last)
generates random integers on the interval [0, n), where the probability of each individual integer i ...
Definition: Discrete.h:1324
Eigen::Rand::binomial
const BinomialType< Derived, Urng > binomial(Index rows, Index cols, Urng &&urng, typename Derived::Scalar trials=1, double p=0.5)
generates reals on the binomial distribution.
Definition: Discrete.h:1441