EigenRand  0.3.0
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[1 << bitsize]);
53  alias = std::unique_ptr<_Size[]>(new _Size[1 << bitsize]);
54 
55  std::copy(o.arr.get(), o.arr.get() + (1 << bitsize), arr.get());
56  std::copy(o.alias.get(), o.alias.get() + (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] = 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] = 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] = 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] = 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  while (1)
264  {
265  _Scalar cands = (_Scalar)(rx & bitmask);
266  if (cands <= pdiff) return cands;
267  if (bitcnt + bitsize < 32)
268  {
269  rx >>= bitsize;
270  bitcnt += bitsize;
271  }
272  else
273  {
274  rx = randbits(rng);
275  bitcnt = bitsize;
276  }
277  }
278  }
279  }
280 
281  template<typename Packet, typename Rng>
282  EIGEN_STRONG_INLINE const Packet packetOp(Rng&& rng)
283  {
284  using namespace Eigen::internal;
285  auto rx = randbits.template packetOp<Packet>(rng);
286  auto pbitmask = pset1<Packet>(bitmask);
287  if (pdiff == bitmask)
288  {
289  return padd(pand(rx, pbitmask), pset1<Packet>(pmin));
290  }
291  else
292  {
293  auto& cm = Rand::detail::CompressMask<sizeof(Packet)>::get_inst();
294  auto plen = pset1<Packet>(pdiff + 1);
295  size_t bitcnt = bitsize;
296  while (1)
297  {
298  // accept cands that only < plen
299  auto cands = pand(rx, pbitmask);
300  bool full = false;
301  cache_rest_cnt = cm.compress_append(cands, pcmplt(cands, plen),
302  OptCacheStore::template get<Packet>(), cache_rest_cnt, full);
303  if (full) return padd(cands, pset1<Packet>(pmin));
304 
305  if (bitcnt + bitsize < 32)
306  {
307  rx = psrl(rx, bitsize);
308  bitcnt += bitsize;
309  }
310  else
311  {
312  rx = randbits.template packetOp<Packet>(rng);
313  bitcnt = bitsize;
314  }
315  }
316  }
317  }
318  };
319 
326  template<typename _Scalar, typename Precision = float>
327  class DiscreteGen;
328 
334  template<typename _Scalar>
335  class DiscreteGen<_Scalar, int32_t> : public GenBase<DiscreteGen<_Scalar, int32_t>, _Scalar>
336  {
337  static_assert(std::is_same<_Scalar, int32_t>::value, "discreteDist needs integral types.");
338 #ifdef EIGEN_VECTORIZE_AVX2
339  OptCacheStore cache;
340  bool valid = false;
341 #endif
342  RandbitsGen<int32_t> randbits;
343  std::vector<uint32_t> cdf;
344  AliasMethod<int32_t, _Scalar> alias_table;
345 
346  public:
347  using Scalar = _Scalar;
348 
356  template<typename RealIter>
357  DiscreteGen(RealIter first, RealIter last)
358  {
359  if (std::distance(first, last) < 16)
360  {
361  // use linear or binary search
362  std::vector<double> _cdf;
363  double acc = 0;
364  for (; first != last; ++first)
365  {
366  _cdf.emplace_back(acc += *first);
367  }
368 
369  for (auto& p : _cdf)
370  {
371  cdf.emplace_back((uint32_t)(p / _cdf.back() * 0x80000000));
372  }
373  }
374  else
375  {
376  // use alias table
377  alias_table = AliasMethod<int32_t, _Scalar>{ first, last };
378  }
379  }
380 
388  template<typename Real,
389  typename std::enable_if<std::is_arithmetic<Real>::value, int>::type = 0>
390  DiscreteGen(const std::initializer_list<Real>& il)
391  : DiscreteGen(il.begin(), il.end())
392  {
393  }
394 
395  DiscreteGen()
396  : DiscreteGen({ 1 })
397  {
398  }
399 
400  DiscreteGen(const DiscreteGen&) = default;
401  DiscreteGen(DiscreteGen&&) = default;
402 
403  DiscreteGen& operator=(const DiscreteGen&) = default;
404  DiscreteGen& operator=(DiscreteGen&&) = default;
405 
406  template<typename Rng>
407  EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng)
408  {
409  using namespace Eigen::internal;
410  if (!cdf.empty())
411  {
412  auto rx = randbits(std::forward<Rng>(rng)) & 0x7FFFFFFF;
413  return (_Scalar)(std::lower_bound(cdf.begin(), cdf.end() - 1, rx) - cdf.begin());
414  }
415  else
416  {
417  auto rx = randbits(std::forward<Rng>(rng));
418  auto albit = rx & alias_table.get_bitmask();
419  uint32_t alx = (uint32_t)(rx >> (sizeof(rx) * 8 - 31));
420  if (alx < alias_table.get_prob()[albit]) return albit;
421  return alias_table.get_alias()[albit];
422  }
423  }
424 
425  template<typename Packet, typename Rng>
426  EIGEN_STRONG_INLINE const Packet packetOp(Rng&& rng)
427  {
428  using namespace Eigen::internal;
429 #ifdef EIGEN_VECTORIZE_AVX2
430  if (valid)
431  {
432  valid = false;
433  return cache.template get<Packet>();
434  }
435  using PacketType = Packet8i;
436 #else
437  using PacketType = Packet;
438 #endif
439 
440  PacketType ret;
441  if (!cdf.empty())
442  {
443  ret = pset1<PacketType>(cdf.size() - 1);
444  auto rx = pand(randbits.template packetOp<PacketType>(std::forward<Rng>(rng)), pset1<PacketType>(0x7FFFFFFF));
445  for (size_t i = 0; i < cdf.size() - 1; ++i)
446  {
447  ret = padd(ret, pcmplt(rx, pset1<PacketType>(cdf[i])));
448  }
449  }
450  else
451  {
452  auto rx = randbits.template packetOp<PacketType>(std::forward<Rng>(rng));
453  auto albit = pand(rx, pset1<PacketType>(alias_table.get_bitmask()));
454  auto c = pcmplt(psrl(rx, 1), pgather(alias_table.get_prob(), albit));
455  ret = pblendv(c, albit, pgather(alias_table.get_alias(), albit));
456  }
457 
458 #ifdef EIGEN_VECTORIZE_AVX2
459  valid = true;
460  cache.template get<Packet>() = _mm256_extractf128_si256(ret, 1);
461  return _mm256_extractf128_si256(ret, 0);
462 #else
463  return ret;
464 #endif
465  }
466  };
467 
473  template<typename _Scalar>
474  class DiscreteGen<_Scalar, float> : public GenBase<DiscreteGen<_Scalar, float>, _Scalar>
475  {
476  static_assert(std::is_same<_Scalar, int32_t>::value, "discreteDist needs integral types.");
477  UniformRealGen<float> ur;
478  std::vector<float> cdf;
479  AliasMethod<float, _Scalar> alias_table;
480 
481  public:
482  using Scalar = _Scalar;
483 
491  template<typename RealIter>
492  DiscreteGen(RealIter first, RealIter last)
493  {
494  if (std::distance(first, last) < 16)
495  {
496  // use linear or binary search
497  std::vector<double> _cdf;
498  double acc = 0;
499  for (; first != last; ++first)
500  {
501  _cdf.emplace_back(acc += *first);
502  }
503 
504  for (auto& p : _cdf)
505  {
506  cdf.emplace_back(p / _cdf.back());
507  }
508  }
509  else
510  {
511  // use alias table
512  alias_table = AliasMethod<float, _Scalar>{ first, last };
513  }
514  }
515 
523  template<typename Real,
524  typename std::enable_if<std::is_arithmetic<Real>::value, int>::type = 0>
525  DiscreteGen(const std::initializer_list<Real>& il)
526  : DiscreteGen(il.begin(), il.end())
527  {
528  }
529 
530  DiscreteGen()
531  : DiscreteGen({ 1 })
532  {
533  }
534 
535  DiscreteGen(const DiscreteGen&) = default;
536  DiscreteGen(DiscreteGen&&) = default;
537 
538  DiscreteGen& operator=(const DiscreteGen&) = default;
539  DiscreteGen& operator=(DiscreteGen&&) = default;
540 
541  template<typename Rng>
542  EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng)
543  {
544  using namespace Eigen::internal;
545  if (!cdf.empty())
546  {
547  auto rx = ur(std::forward<Rng>(rng));
548  return (_Scalar)(std::lower_bound(cdf.begin(), cdf.end() - 1, rx) - cdf.begin());
549  }
550  else
551  {
552  auto albit = pfirst(rng()) & alias_table.get_bitmask();
553  auto alx = ur(rng);
554  if (alx < alias_table.get_prob()[albit]) return albit;
555  return alias_table.get_alias()[albit];
556  }
557  }
558 
559  template<typename Packet, typename Rng>
560  EIGEN_STRONG_INLINE const Packet packetOp(Rng&& rng)
561  {
562  using namespace Eigen::internal;
563  using PacketType = decltype(reinterpret_to_float(std::declval<Packet>()));
564 
565  if (!cdf.empty())
566  {
567  auto ret = pset1<Packet>(cdf.size());
568  auto rx = ur.template packetOp<PacketType>(std::forward<Rng>(rng));
569  for (auto& p : cdf)
570  {
571  ret = padd(ret, reinterpret_to_int(pcmplt(rx, pset1<PacketType>(p))));
572  }
573  return ret;
574  }
575  else
576  {
577  using RUtils = RawbitsMaker<Packet, Rng>;
578  auto albit = pand(RUtils{}.rawbits(rng), pset1<Packet>(alias_table.get_bitmask()));
579  auto c = reinterpret_to_int(pcmplt(ur.template packetOp<PacketType>(rng), pgather(alias_table.get_prob(), albit)));
580  return pblendv(c, albit, pgather(alias_table.get_alias(), albit));
581  }
582  }
583  };
584 
590  template<typename _Scalar>
591  class DiscreteGen<_Scalar, double> : public GenBase<DiscreteGen<_Scalar, double>, _Scalar>
592  {
593  static_assert(std::is_same<_Scalar, int32_t>::value, "discreteDist needs integral types.");
594  UniformRealGen<double> ur;
595  std::vector<double> cdf;
596  AliasMethod<double, _Scalar> alias_table;
597 
598  public:
599  using Scalar = _Scalar;
600 
608  template<typename RealIter>
609  DiscreteGen(RealIter first, RealIter last)
610  {
611  if (std::distance(first, last) < 16)
612  {
613  // use linear or binary search
614  std::vector<double> _cdf;
615  double acc = 0;
616  for (; first != last; ++first)
617  {
618  _cdf.emplace_back(acc += *first);
619  }
620 
621  for (auto& p : _cdf)
622  {
623  cdf.emplace_back(p / _cdf.back());
624  }
625  }
626  else
627  {
628  // use alias table
629  alias_table = AliasMethod<double, _Scalar>{ first, last };
630  }
631  }
632 
640  template<typename Real,
641  typename std::enable_if<std::is_arithmetic<Real>::value, int>::type = 0>
642  DiscreteGen(const std::initializer_list<Real>& il)
643  : DiscreteGen(il.begin(), il.end())
644  {
645  }
646 
647  DiscreteGen()
648  : DiscreteGen({ 1 })
649  {
650  }
651 
652  DiscreteGen(const DiscreteGen&) = default;
653  DiscreteGen(DiscreteGen&&) = default;
654 
655  DiscreteGen& operator=(const DiscreteGen&) = default;
656  DiscreteGen& operator=(DiscreteGen&&) = default;
657 
658  template<typename Rng>
659  EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng)
660  {
661  using namespace Eigen::internal;
662  if (!cdf.empty())
663  {
664  auto rx = ur(std::forward<Rng>(rng));
665  return (_Scalar)(std::lower_bound(cdf.begin(), cdf.end() - 1, rx) - cdf.begin());
666  }
667  else
668  {
669  auto albit = pfirst(rng()) & alias_table.get_bitmask();
670  auto alx = ur(rng);
671  if (alx < alias_table.get_prob()[albit]) return albit;
672  return alias_table.get_alias()[albit];
673  }
674  }
675 
676  template<typename Packet, typename Rng>
677  EIGEN_STRONG_INLINE const Packet packetOp(Rng&& rng)
678  {
679  using namespace Eigen::internal;
680  using DPacket = decltype(reinterpret_to_double(std::declval<Packet>()));
681  if (!cdf.empty())
682  {
683  auto ret = pset1<Packet>(cdf.size());
684 #ifdef EIGEN_VECTORIZE_AVX
685  auto rx = ur.template packetOp<Packet4d>(std::forward<Rng>(rng));
686  for (auto& p : cdf)
687  {
688  auto c = reinterpret_to_int(pcmplt(rx, pset1<decltype(rx)>(p)));
689  auto r = combine_low32(c);
690  ret = padd(ret, r);
691  }
692 #else
693  auto rx1 = ur.template packetOp<DPacket>(rng),
694  rx2 = ur.template packetOp<DPacket>(rng);
695  for (auto& p : cdf)
696  {
697  auto pp = pset1<decltype(rx1)>(p);
698  ret = padd(ret, combine_low32(reinterpret_to_int(pcmplt(rx1, pp)), reinterpret_to_int(pcmplt(rx2, pp))));
699  }
700 #endif
701  return ret;
702  }
703  else
704  {
705 #ifdef EIGEN_VECTORIZE_AVX
706  using RUtils = RawbitsMaker<Packet, Rng>;
707  auto albit = pand(RUtils{}.rawbits(rng), pset1<Packet>(alias_table.get_bitmask()));
708  auto c = reinterpret_to_int(pcmplt(ur.template packetOp<Packet4d>(rng), pgather(alias_table.get_prob(), _mm256_castsi128_si256(albit))));
709  return pblendv(combine_low32(c), albit, pgather(alias_table.get_alias(), albit));
710 #else
711  using RUtils = RawbitsMaker<Packet, Rng>;
712  auto albit = pand(RUtils{}.rawbits(rng), pset1<Packet>(alias_table.get_bitmask()));
713  auto c1 = reinterpret_to_int(pcmplt(ur.template packetOp<DPacket>(rng), pgather(alias_table.get_prob(), albit)));
714  auto c2 = reinterpret_to_int(pcmplt(ur.template packetOp<DPacket>(rng), pgather(alias_table.get_prob(), albit, true)));
715  return pblendv(combine_low32(c1, c2), albit, pgather(alias_table.get_alias(), albit));
716 #endif
717  }
718  }
719  };
720 
721  template<typename> class BinomialGen;
722 
728  template<typename _Scalar>
729  class PoissonGen : OptCacheStore, public GenBase<PoissonGen<_Scalar>, _Scalar>
730  {
731  friend BinomialGen<_Scalar>;
732  static_assert(std::is_same<_Scalar, int32_t>::value, "poisson needs integral types.");
733  int cache_rest_cnt = 0;
734  UniformRealGen<float> ur;
735 
736  protected:
737  double mean, ne_mean, sqrt_tmean, log_mean, g1;
738 
739  public:
740  using Scalar = _Scalar;
741 
747  PoissonGen(double _mean = 1)
748  : mean{ _mean }, ne_mean{ std::exp(-_mean) }
749  {
750  sqrt_tmean = std::sqrt(2 * mean);
751  log_mean = std::log(mean);
752  g1 = mean * log_mean - std::lgamma(mean + 1);
753  }
754 
755  PoissonGen(const PoissonGen&) = default;
756  PoissonGen(PoissonGen&&) = default;
757 
758  PoissonGen& operator=(const PoissonGen&) = default;
759  PoissonGen& operator=(PoissonGen&&) = default;
760 
761  template<typename Rng>
762  EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng)
763  {
764  using namespace Eigen::internal;
765  if (mean < 12)
766  {
767  _Scalar res = 0;
768  double val = 1;
769  for (; ; ++res)
770  {
771  val *= ur(rng);
772  if (val <= ne_mean) break;
773  }
774  return res;
775  }
776  else
777  {
778  _Scalar res;
779  double yx;
780  while (1)
781  {
782  yx = std::tan(constant::pi * ur(rng));
783  res = (_Scalar)(sqrt_tmean * yx + mean);
784  if (res >= 0 && ur(rng) <= 0.9 * (1.0 + yx * yx)
785  * std::exp(res * log_mean - std::lgamma(res + 1.0) - g1))
786  {
787  return res;
788  }
789  }
790  }
791  }
792 
793  template<typename Packet, typename Rng>
794  EIGEN_STRONG_INLINE const Packet packetOp(Rng&& rng)
795  {
796  using namespace Eigen::internal;
797  using PacketType = decltype(reinterpret_to_float(std::declval<Packet>()));
798 
799  if (mean < 12)
800  {
801  Packet res = pset1<Packet>(0);
802  PacketType val = pset1<PacketType>(1), pne_mean = pset1<PacketType>(ne_mean);
803  while (1)
804  {
805  val = pmul(val, ur.template packetOp<PacketType>(rng));
806  auto c = reinterpret_to_int(pcmplt(pne_mean, val));
807  if (pmovemask(c) == 0) break;
808  res = padd(res, pnegate(c));
809  }
810  return res;
811  }
812  else
813  {
814  auto& cm = Rand::detail::CompressMask<sizeof(Packet)>::get_inst();
815  const PacketType ppi = pset1<PacketType>(constant::pi),
816  psqrt_tmean = pset1<PacketType>(sqrt_tmean),
817  pmean = pset1<PacketType>(mean),
818  plog_mean = pset1<PacketType>(log_mean),
819  pg1 = pset1<PacketType>(g1);
820  while (1)
821  {
822  PacketType fres, yx, psin, pcos;
823  psincos(pmul(ppi, ur.template packetOp<PacketType>(rng)), psin, pcos);
824  yx = pdiv(psin, pcos);
825  fres = ptruncate(padd(pmul(psqrt_tmean, yx), pmean));
826 
827  auto p1 = pmul(padd(pmul(yx, yx), pset1<PacketType>(1)), pset1<PacketType>(0.9));
828  auto p2 = pexp(psub(psub(pmul(fres, plog_mean), plgamma_approx(padd(fres, pset1<PacketType>(1)))), pg1));
829 
830  auto c1 = pcmple(pset1<PacketType>(0), fres);
831  auto c2 = pcmple(ur.template packetOp<PacketType>(rng), pmul(p1, p2));
832 
833  auto cands = fres;
834  bool full = false;
835  cache_rest_cnt = cm.compress_append(cands, pand(c1, c2),
836  OptCacheStore::template get<PacketType>(), cache_rest_cnt, full);
837  if (full) return pcast<PacketType, Packet>(cands);
838  }
839  }
840  }
841  };
842 
848  template<typename _Scalar>
849  class BinomialGen : public GenBase<BinomialGen<_Scalar>, _Scalar>
850  {
851  static_assert(std::is_same<_Scalar, int32_t>::value, "binomial needs integral types.");
852 
853  PoissonGen<_Scalar> poisson;
854  _Scalar trials;
855  double p, small_p, g1, sqrt_v, log_small_p, log_small_q;
856 
857  public:
858  using Scalar = _Scalar;
859 
866  BinomialGen(_Scalar _trials = 1, double _p = 0.5)
867  : poisson{ _trials * std::min(_p, 1 - _p) },
868  trials{ _trials }, p{ _p }, small_p{ std::min(p, 1 - p) }
869 
870  {
871  if (!(trials < 25 || poisson.mean < 1.0))
872  {
873  g1 = std::lgamma(trials + 1);
874  sqrt_v = std::sqrt(2 * poisson.mean * (1 - small_p));
875  log_small_p = std::log(small_p);
876  log_small_q = std::log(1 - small_p);
877  }
878  }
879 
880  BinomialGen(const BinomialGen&) = default;
881  BinomialGen(BinomialGen&&) = default;
882 
883  BinomialGen& operator=(const BinomialGen&) = default;
884  BinomialGen& operator=(BinomialGen&&) = default;
885 
886  template<typename Rng>
887  EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng)
888  {
889  using namespace Eigen::internal;
890  _Scalar res;
891  if (trials < 25)
892  {
893  res = 0;
894  for (int i = 0; i < trials; ++i)
895  {
896  if (poisson.ur(rng) < p) ++res;
897  }
898  return res;
899  }
900  else if (poisson.mean < 1.0)
901  {
902  res = poisson(rng);
903  }
904  else
905  {
906  while (1)
907  {
908  double ys;
909  ys = std::tan(constant::pi * poisson.ur(rng));
910  res = (_Scalar)(sqrt_v * ys + poisson.mean);
911  if (0 <= res && res <= trials && poisson.ur(rng) <= 1.2 * sqrt_v
912  * (1.0 + ys * ys)
913  * std::exp(g1 - std::lgamma(res + 1)
914  - std::lgamma(trials - res + 1.0)
915  + res * log_small_p
916  + (trials - res) * log_small_q)
917  )
918  {
919  break;
920  }
921  }
922  }
923  return p == small_p ? res : trials - res;
924  }
925 
926  template<typename Packet, typename Rng>
927  EIGEN_STRONG_INLINE const Packet packetOp(Rng&& rng)
928  {
929  using namespace Eigen::internal;
930  using PacketType = decltype(reinterpret_to_float(std::declval<Packet>()));
931  Packet res;
932  if (trials < 25)
933  {
934  PacketType pp = pset1<PacketType>(p);
935  res = pset1<Packet>(trials);
936  for (int i = 0; i < trials; ++i)
937  {
938  auto c = reinterpret_to_int(pcmple(pp, poisson.ur.template packetOp<PacketType>(rng)));
939  res = padd(res, c);
940  }
941  return res;
942  }
943  else if (poisson.mean < 1.0)
944  {
945  res = poisson.template packetOp<Packet>(rng);
946  }
947  else
948  {
949  auto& cm = Rand::detail::CompressMask<sizeof(Packet)>::get_inst();
950  const PacketType ppi = pset1<PacketType>(constant::pi),
951  ptrials = pset1<PacketType>(trials),
952  psqrt_v = pset1<PacketType>(sqrt_v),
953  pmean = pset1<PacketType>(poisson.mean),
954  plog_small_p = pset1<PacketType>(log_small_p),
955  plog_small_q = pset1<PacketType>(log_small_q),
956  pg1 = pset1<PacketType>(g1);
957  while (1)
958  {
959  PacketType fres, ys, psin, pcos;
960  psincos(pmul(ppi, poisson.ur.template packetOp<PacketType>(rng)), psin, pcos);
961  ys = pdiv(psin, pcos);
962  fres = ptruncate(padd(pmul(psqrt_v, ys), pmean));
963 
964  auto p1 = pmul(pmul(pset1<PacketType>(1.2), psqrt_v), padd(pset1<PacketType>(1), pmul(ys, ys)));
965  auto p2 = pexp(
966  padd(padd(psub(
967  psub(pg1, plgamma_approx(padd(fres, pset1<PacketType>(1)))),
968  plgamma_approx(psub(padd(ptrials, pset1<PacketType>(1)), fres))
969  ), pmul(fres, plog_small_p)), pmul(psub(ptrials, fres), plog_small_q))
970  );
971 
972  auto c1 = pand(pcmple(pset1<PacketType>(0), fres), pcmple(fres, ptrials));
973  auto c2 = pcmple(poisson.ur.template packetOp<PacketType>(rng), pmul(p1, p2));
974 
975  auto cands = fres;
976  bool full = false;
977  poisson.cache_rest_cnt = cm.compress_append(cands, pand(c1, c2),
978  poisson.template get<PacketType>(), poisson.cache_rest_cnt, full);
979  if (full)
980  {
981  res = pcast<PacketType, Packet>(cands);
982  break;
983  }
984  }
985  }
986  return p == small_p ? res : psub(pset1<Packet>(trials), res);
987  }
988  };
989 
995  template<typename _Scalar>
996  class GeometricGen : public GenBase<GeometricGen<_Scalar>, _Scalar>
997  {
998  static_assert(std::is_same<_Scalar, int32_t>::value, "geomtric needs integral types.");
999  UniformRealGen<float> ur;
1000  double p, rlog_q;
1001 
1002  public:
1003  using Scalar = _Scalar;
1004 
1010  GeometricGen(double _p = 0.5)
1011  : p{ _p }, rlog_q{ 1 / std::log(1 - p) }
1012  {
1013  }
1014 
1015  GeometricGen(const GeometricGen&) = default;
1016  GeometricGen(GeometricGen&&) = default;
1017 
1018  GeometricGen& operator=(const GeometricGen&) = default;
1019  GeometricGen& operator=(GeometricGen&&) = default;
1020 
1021  template<typename Rng>
1022  EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng)
1023  {
1024  using namespace Eigen::internal;
1025  return (_Scalar)(std::log(1 - ur(std::forward<Rng>(rng))) * rlog_q);
1026  }
1027 
1028  template<typename Packet, typename Rng>
1029  EIGEN_STRONG_INLINE const Packet packetOp(Rng&& rng)
1030  {
1031  using namespace Eigen::internal;
1032  using PacketType = decltype(reinterpret_to_float(std::declval<Packet>()));
1033 
1034  return pcast<PacketType, Packet>(ptruncate(pmul(plog(
1035  psub(pset1<PacketType>(1), ur.template packetOp<PacketType>(std::forward<Rng>(rng)))
1036  ), pset1<PacketType>(rlog_q))));
1037  }
1038  };
1039 
1040 
1041  template<typename Derived, typename Urng>
1042  using UniformIntType = CwiseNullaryOp<internal::scalar_rng_adaptor<UniformIntGen<typename Derived::Scalar>, typename Derived::Scalar, Urng, true>, const Derived>;
1043 
1057  template<typename Derived, typename Urng>
1058  inline const UniformIntType<Derived, Urng>
1059  uniformInt(Index rows, Index cols, Urng&& urng, typename Derived::Scalar min, typename Derived::Scalar max)
1060  {
1061  return {
1062  rows, cols, { std::forward<Urng>(urng), UniformIntGen<typename Derived::Scalar>{ min, max } }
1063  };
1064  }
1065 
1078  template<typename Derived, typename Urng>
1079  inline const UniformIntType<Derived, Urng>
1080  uniformIntLike(Derived& o, Urng&& urng, typename Derived::Scalar min, typename Derived::Scalar max)
1081  {
1082  return {
1083  o.rows(), o.cols(), { std::forward<Urng>(urng), UniformIntGen<typename Derived::Scalar>{ min, max } }
1084  };
1085  }
1086 
1087  template<typename Derived, typename Urng>
1088  using DiscreteFType = CwiseNullaryOp<internal::scalar_rng_adaptor<DiscreteGen<typename Derived::Scalar, float>, typename Derived::Scalar, Urng, true>, const Derived>;
1089 
1104  template<typename Derived, typename Urng, typename RealIter>
1105  inline const DiscreteFType<Derived, Urng>
1106  discreteF(Index rows, Index cols, Urng&& urng, RealIter first, RealIter last)
1107  {
1108  return {
1109  rows, cols, { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, float>{first, last} }
1110  };
1111  }
1112 
1126  template<typename Derived, typename Urng, typename RealIter>
1127  inline const DiscreteFType<Derived, Urng>
1128  discreteFLike(Derived& o, Urng&& urng, RealIter first, RealIter last)
1129  {
1130  return {
1131  o.rows(), o.cols(), { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, float>(first, last) }
1132  };
1133  }
1134 
1149  template<typename Derived, typename Urng, typename Real>
1150  inline const DiscreteFType<Derived, Urng>
1151  discreteF(Index rows, Index cols, Urng&& urng, const std::initializer_list<Real>& il)
1152  {
1153  return {
1154  rows, cols, { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, float>{il.begin(), il.end()} }
1155  };
1156  }
1157 
1171  template<typename Derived, typename Urng, typename Real>
1172  inline const DiscreteFType<Derived, Urng>
1173  discreteFLike(Derived& o, Urng&& urng, const std::initializer_list<Real>& il)
1174  {
1175  return {
1176  o.rows(), o.cols(), { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, float>(il.begin(), il.end()) }
1177  };
1178  }
1179 
1180  template<typename Derived, typename Urng>
1181  using DiscreteDType = CwiseNullaryOp<internal::scalar_rng_adaptor<DiscreteGen<typename Derived::Scalar, double>, typename Derived::Scalar, Urng, true>, const Derived>;
1182 
1197  template<typename Derived, typename Urng, typename RealIter>
1198  inline const DiscreteDType<Derived, Urng>
1199  discreteD(Index rows, Index cols, Urng&& urng, RealIter first, RealIter last)
1200  {
1201  return {
1202  rows, cols, { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, double>{first, last} }
1203  };
1204  }
1205 
1219  template<typename Derived, typename Urng, typename RealIter>
1220  inline const DiscreteDType<Derived, Urng>
1221  discreteDLike(Derived& o, Urng&& urng, RealIter first, RealIter last)
1222  {
1223  return {
1224  o.rows(), o.cols(), { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, double>{first, last} }
1225  };
1226  }
1227 
1242  template<typename Derived, typename Urng, typename Real>
1243  inline const DiscreteDType<Derived, Urng>
1244  discreteD(Index rows, Index cols, Urng&& urng, const std::initializer_list<Real>& il)
1245  {
1246  return {
1247  rows, cols, { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, double>{il.begin(), il.end()} }
1248  };
1249  }
1250 
1264  template<typename Derived, typename Urng, typename Real>
1265  inline const DiscreteDType<Derived, Urng>
1266  discreteDLike(Derived& o, Urng&& urng, const std::initializer_list<Real>& il)
1267  {
1268  return {
1269  o.rows(), o.cols(), { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, double>{il.begin(), il.end()} }
1270  };
1271  }
1272 
1273  template<typename Derived, typename Urng>
1274  using DiscreteType = CwiseNullaryOp<internal::scalar_rng_adaptor<DiscreteGen<typename Derived::Scalar, int32_t>, typename Derived::Scalar, Urng, true>, const Derived>;
1275 
1290  template<typename Derived, typename Urng, typename RealIter>
1291  inline const DiscreteType<Derived, Urng>
1292  discrete(Index rows, Index cols, Urng&& urng, RealIter first, RealIter last)
1293  {
1294  return {
1295  rows, cols, { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, int32_t>{first, last} }
1296  };
1297  }
1298 
1312  template<typename Derived, typename Urng, typename RealIter>
1313  inline const DiscreteType<Derived, Urng>
1314  discreteLike(Derived& o, Urng&& urng, RealIter first, RealIter last)
1315  {
1316  return {
1317  o.rows(), o.cols(), { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, int32_t>{first, last} }
1318  };
1319  }
1320 
1335  template<typename Derived, typename Urng, typename Real>
1336  inline const DiscreteType<Derived, Urng>
1337  discrete(Index rows, Index cols, Urng&& urng, const std::initializer_list<Real>& il)
1338  {
1339  return {
1340  rows, cols, { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, int32_t>{il.begin(), il.end()} }
1341  };
1342  }
1343 
1357  template<typename Derived, typename Urng, typename Real>
1358  inline const DiscreteType<Derived, Urng>
1359  discreteLike(Derived& o, Urng&& urng, const std::initializer_list<Real>& il)
1360  {
1361  return {
1362  o.rows(), o.cols(), { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, int32_t>{il.begin(), il.end()} }
1363  };
1364  }
1365 
1366  template<typename Derived, typename Urng>
1367  using PoissonType = CwiseNullaryOp<internal::scalar_rng_adaptor<PoissonGen<typename Derived::Scalar>, typename Derived::Scalar, Urng, true>, const Derived>;
1368 
1382  template<typename Derived, typename Urng>
1383  inline const PoissonType<Derived, Urng>
1384  poisson(Index rows, Index cols, Urng&& urng, double mean = 1)
1385  {
1386  return {
1387  rows, cols, { std::forward<Urng>(urng), PoissonGen<typename Derived::Scalar>{mean} }
1388  };
1389  }
1390 
1403  template<typename Derived, typename Urng>
1404  inline const PoissonType<Derived, Urng>
1405  poissonLike(Derived& o, Urng&& urng, double mean = 1)
1406  {
1407  return {
1408  o.rows(), o.cols(), { std::forward<Urng>(urng), PoissonGen<typename Derived::Scalar>{mean} }
1409  };
1410  }
1411 
1412  template<typename Derived, typename Urng>
1413  using BinomialType = CwiseNullaryOp<internal::scalar_rng_adaptor<BinomialGen<typename Derived::Scalar>, typename Derived::Scalar, Urng, true>, const Derived>;
1414 
1429  template<typename Derived, typename Urng>
1430  inline const BinomialType<Derived, Urng>
1431  binomial(Index rows, Index cols, Urng&& urng, typename Derived::Scalar trials = 1, double p = 0.5)
1432  {
1433  return {
1434  rows, cols, { std::forward<Urng>(urng), BinomialGen<typename Derived::Scalar>{trials, p} }
1435  };
1436  }
1437 
1451  template<typename Derived, typename Urng>
1452  inline const BinomialType<Derived, Urng>
1453  binomialLike(Derived& o, Urng&& urng, typename Derived::Scalar trials = 1, double p = 0.5)
1454  {
1455  return {
1456  o.rows(), o.cols(), { std::forward<Urng>(urng), BinomialGen<typename Derived::Scalar>{trials, p} }
1457  };
1458  }
1459 
1460  template<typename Derived, typename Urng>
1461  using GeometricType = CwiseNullaryOp<internal::scalar_rng_adaptor<GeometricGen<typename Derived::Scalar>, typename Derived::Scalar, Urng, true>, const Derived>;
1462 
1476  template<typename Derived, typename Urng>
1477  inline const GeometricType<Derived, Urng>
1478  geometric(Index rows, Index cols, Urng&& urng, double p = 0.5)
1479  {
1480  return {
1481  rows, cols, { std::forward<Urng>(urng), GeometricGen<typename Derived::Scalar>{p} }
1482  };
1483  }
1484 
1497  template<typename Derived, typename Urng>
1498  inline const GeometricType<Derived, Urng>
1499  geometricLike(Derived& o, Urng&& urng, double p = 0.5)
1500  {
1501  return {
1502  o.rows(), o.cols(), { std::forward<Urng>(urng), GeometricGen<typename Derived::Scalar>{p} }
1503  };
1504  }
1505  }
1506 }
1507 
1508 #endif
Eigen::Rand::BinomialGen
Generator of integers on a binomial distribution.
Definition: Discrete.h:850
Eigen::Rand::BinomialGen::BinomialGen
BinomialGen(_Scalar _trials=1, double _p=0.5)
Construct a new Binomial Generator.
Definition: Discrete.h:866
Eigen::Rand::DiscreteGen< _Scalar, double >::DiscreteGen
DiscreteGen(RealIter first, RealIter last)
Construct a new Discrete Generator.
Definition: Discrete.h:609
Eigen::Rand::DiscreteGen< _Scalar, double >::DiscreteGen
DiscreteGen(const std::initializer_list< Real > &il)
Construct a new Discrete Generator.
Definition: Discrete.h:642
Eigen::Rand::DiscreteGen< _Scalar, float >::DiscreteGen
DiscreteGen(RealIter first, RealIter last)
Construct a new Discrete Generator.
Definition: Discrete.h:492
Eigen::Rand::DiscreteGen< _Scalar, float >::DiscreteGen
DiscreteGen(const std::initializer_list< Real > &il)
Construct a new Discrete Generator.
Definition: Discrete.h:525
Eigen::Rand::DiscreteGen< _Scalar, int32_t >::DiscreteGen
DiscreteGen(RealIter first, RealIter last)
Construct a new Discrete Generator.
Definition: Discrete.h:357
Eigen::Rand::DiscreteGen< _Scalar, int32_t >::DiscreteGen
DiscreteGen(const std::initializer_list< Real > &il)
Construct a new Discrete Generator.
Definition: Discrete.h:390
Eigen::Rand::DiscreteGen
Generator of integers on the interval [0, n), where the probability of each individual integer i is p...
Definition: Discrete.h:327
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:997
Eigen::Rand::GeometricGen::GeometricGen
GeometricGen(double _p=0.5)
Construct a new Geometric Generator.
Definition: Discrete.h:1010
Eigen::Rand::PoissonGen
Generator of integers on a Poisson distribution.
Definition: Discrete.h:730
Eigen::Rand::PoissonGen::PoissonGen
PoissonGen(double _mean=1)
Construct a new Poisson Generator.
Definition: Discrete.h:747
Eigen::Rand::RandbitsGen
Generator of random bits for integral scalars.
Definition: Basic.h:248
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:1405
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:1106
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:1499
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:1453
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:1059
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:1221
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:1128
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:1080
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:1199
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:1384
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:1292
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:1478
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:1314
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:1431