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, _Scalar _max)
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  template<typename Rng>
249  EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng)
250  {
251  using namespace Eigen::internal;
252  auto rx = randbits(rng);
253  if (pdiff == bitmask)
254  {
255  return (_Scalar)(rx & bitmask) + pmin;
256  }
257  else
258  {
259  size_t bitcnt = bitsize;
260  while (1)
261  {
262  _Scalar cands = (_Scalar)(rx & bitmask);
263  if (cands <= pdiff) return cands;
264  if (bitcnt + bitsize < 32)
265  {
266  rx >>= bitsize;
267  bitcnt += bitsize;
268  }
269  else
270  {
271  rx = randbits(rng);
272  bitcnt = bitsize;
273  }
274  }
275  }
276  }
277 
278  template<typename Packet, typename Rng>
279  EIGEN_STRONG_INLINE const Packet packetOp(Rng&& rng)
280  {
281  using namespace Eigen::internal;
282  auto rx = randbits.template packetOp<Packet>(rng);
283  auto pbitmask = pset1<Packet>(bitmask);
284  if (pdiff == bitmask)
285  {
286  return padd(pand(rx, pbitmask), pset1<Packet>(pmin));
287  }
288  else
289  {
290  auto& cm = Rand::detail::CompressMask<sizeof(Packet)>::get_inst();
291  auto plen = pset1<Packet>(pdiff + 1);
292  size_t bitcnt = bitsize;
293  while (1)
294  {
295  // accept cands that only < plen
296  auto cands = pand(rx, pbitmask);
297  bool full = false;
298  cache_rest_cnt = cm.compress_append(cands, pcmplt(cands, plen),
299  OptCacheStore::template get<Packet>(), cache_rest_cnt, full);
300  if (full) return padd(cands, pset1<Packet>(pmin));
301 
302  if (bitcnt + bitsize < 32)
303  {
304  rx = psrl(rx, bitsize);
305  bitcnt += bitsize;
306  }
307  else
308  {
309  rx = randbits.template packetOp<Packet>(rng);
310  bitcnt = bitsize;
311  }
312  }
313  }
314  }
315  };
316 
323  template<typename _Scalar, typename Precision = float>
324  class DiscreteGen;
325 
331  template<typename _Scalar>
332  class DiscreteGen<_Scalar, int32_t> : public GenBase<DiscreteGen<_Scalar, int32_t>, _Scalar>
333  {
334  static_assert(std::is_same<_Scalar, int32_t>::value, "discreteDist needs integral types.");
335 #ifdef EIGEN_VECTORIZE_AVX2
336  OptCacheStore cache;
337  bool valid = false;
338 #endif
339  RandbitsGen<int32_t> randbits;
340  std::vector<uint32_t> cdf;
341  AliasMethod<int32_t, _Scalar> alias_table;
342 
343  public:
344  using Scalar = _Scalar;
345 
353  template<typename RealIter>
354  DiscreteGen(RealIter first, RealIter last)
355  {
356  if (std::distance(first, last) < 16)
357  {
358  // use linear or binary search
359  std::vector<double> _cdf;
360  double acc = 0;
361  for (; first != last; ++first)
362  {
363  _cdf.emplace_back(acc += *first);
364  }
365 
366  for (auto& p : _cdf)
367  {
368  cdf.emplace_back((uint32_t)(p / _cdf.back() * 0x80000000));
369  }
370  }
371  else
372  {
373  // use alias table
374  alias_table = AliasMethod<int32_t, _Scalar>{ first, last };
375  }
376  }
377 
385  template<typename Real,
386  typename std::enable_if<std::is_arithmetic<Real>::value, int>::type = 0>
387  DiscreteGen(const std::initializer_list<Real>& il)
388  : DiscreteGen(il.begin(), il.end())
389  {
390  }
391 
392  template<typename Rng>
393  EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng)
394  {
395  using namespace Eigen::internal;
396  if (!cdf.empty())
397  {
398  auto rx = randbits(std::forward<Rng>(rng)) & 0x7FFFFFFF;
399  return (_Scalar)(std::lower_bound(cdf.begin(), cdf.end() - 1, rx) - cdf.begin());
400  }
401  else
402  {
403  auto rx = randbits(std::forward<Rng>(rng));
404  auto albit = rx & alias_table.get_bitmask();
405  uint32_t alx = (uint32_t)(rx >> (sizeof(rx) * 8 - 31));
406  if (alx < alias_table.get_prob()[albit]) return albit;
407  return alias_table.get_alias()[albit];
408  }
409  }
410 
411  template<typename Packet, typename Rng>
412  EIGEN_STRONG_INLINE const Packet packetOp(Rng&& rng)
413  {
414  using namespace Eigen::internal;
415 #ifdef EIGEN_VECTORIZE_AVX2
416  if (valid)
417  {
418  valid = false;
419  return cache.template get<Packet>();
420  }
421  using PacketType = Packet8i;
422 #else
423  using PacketType = Packet;
424 #endif
425 
426  PacketType ret;
427  if (!cdf.empty())
428  {
429  ret = pset1<PacketType>(cdf.size() - 1);
430  auto rx = pand(randbits.template packetOp<PacketType>(std::forward<Rng>(rng)), pset1<PacketType>(0x7FFFFFFF));
431  for (size_t i = 0; i < cdf.size() - 1; ++i)
432  {
433  ret = padd(ret, pcmplt(rx, pset1<PacketType>(cdf[i])));
434  }
435  }
436  else
437  {
438  auto rx = randbits.template packetOp<PacketType>(std::forward<Rng>(rng));
439  auto albit = pand(rx, pset1<PacketType>(alias_table.get_bitmask()));
440  auto c = pcmplt(psrl(rx, 1), pgather(alias_table.get_prob(), albit));
441  ret = pblendv(c, albit, pgather(alias_table.get_alias(), albit));
442  }
443 
444 #ifdef EIGEN_VECTORIZE_AVX2
445  valid = true;
446  cache.template get<Packet>() = _mm256_extractf128_si256(ret, 1);
447  return _mm256_extractf128_si256(ret, 0);
448 #else
449  return ret;
450 #endif
451  }
452  };
453 
459  template<typename _Scalar>
460  class DiscreteGen<_Scalar, float> : public GenBase<DiscreteGen<_Scalar, float>, _Scalar>
461  {
462  static_assert(std::is_same<_Scalar, int32_t>::value, "discreteDist needs integral types.");
463  UniformRealGen<float> ur;
464  std::vector<float> cdf;
465  AliasMethod<float, _Scalar> alias_table;
466 
467  public:
468  using Scalar = _Scalar;
469 
477  template<typename RealIter>
478  DiscreteGen(RealIter first, RealIter last)
479  {
480  if (std::distance(first, last) < 16)
481  {
482  // use linear or binary search
483  std::vector<double> _cdf;
484  double acc = 0;
485  for (; first != last; ++first)
486  {
487  _cdf.emplace_back(acc += *first);
488  }
489 
490  for (auto& p : _cdf)
491  {
492  cdf.emplace_back(p / _cdf.back());
493  }
494  }
495  else
496  {
497  // use alias table
498  alias_table = AliasMethod<float, _Scalar>{ first, last };
499  }
500  }
501 
509  template<typename Real,
510  typename std::enable_if<std::is_arithmetic<Real>::value, int>::type = 0>
511  DiscreteGen(const std::initializer_list<Real>& il)
512  : DiscreteGen(il.begin(), il.end())
513  {
514  }
515 
516  template<typename Rng>
517  EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng)
518  {
519  using namespace Eigen::internal;
520  if (!cdf.empty())
521  {
522  auto rx = ur(std::forward<Rng>(rng));
523  return (_Scalar)(std::lower_bound(cdf.begin(), cdf.end() - 1, rx) - cdf.begin());
524  }
525  else
526  {
527  auto albit = pfirst(rng()) & alias_table.get_bitmask();
528  auto alx = ur(rng);
529  if (alx < alias_table.get_prob()[albit]) return albit;
530  return alias_table.get_alias()[albit];
531  }
532  }
533 
534  template<typename Packet, typename Rng>
535  EIGEN_STRONG_INLINE const Packet packetOp(Rng&& rng)
536  {
537  using namespace Eigen::internal;
538  using PacketType = decltype(reinterpret_to_float(std::declval<Packet>()));
539 
540  if (!cdf.empty())
541  {
542  auto ret = pset1<Packet>(cdf.size());
543  auto rx = ur.template packetOp<PacketType>(std::forward<Rng>(rng));
544  for (auto& p : cdf)
545  {
546  ret = padd(ret, reinterpret_to_int(pcmplt(rx, pset1<PacketType>(p))));
547  }
548  return ret;
549  }
550  else
551  {
552  using RUtils = RawbitsMaker<Packet, Rng>;
553  auto albit = pand(RUtils{}.rawbits(rng), pset1<Packet>(alias_table.get_bitmask()));
554  auto c = reinterpret_to_int(pcmplt(ur.template packetOp<PacketType>(rng), pgather(alias_table.get_prob(), albit)));
555  return pblendv(c, albit, pgather(alias_table.get_alias(), albit));
556  }
557  }
558  };
559 
565  template<typename _Scalar>
566  class DiscreteGen<_Scalar, double> : public GenBase<DiscreteGen<_Scalar, double>, _Scalar>
567  {
568  static_assert(std::is_same<_Scalar, int32_t>::value, "discreteDist needs integral types.");
569  UniformRealGen<double> ur;
570  std::vector<double> cdf;
571  AliasMethod<double, _Scalar> alias_table;
572 
573  public:
574  using Scalar = _Scalar;
575 
583  template<typename RealIter>
584  DiscreteGen(RealIter first, RealIter last)
585  {
586  if (std::distance(first, last) < 16)
587  {
588  // use linear or binary search
589  std::vector<double> _cdf;
590  double acc = 0;
591  for (; first != last; ++first)
592  {
593  _cdf.emplace_back(acc += *first);
594  }
595 
596  for (auto& p : _cdf)
597  {
598  cdf.emplace_back(p / _cdf.back());
599  }
600  }
601  else
602  {
603  // use alias table
604  alias_table = AliasMethod<double, _Scalar>{ first, last };
605  }
606  }
607 
615  template<typename Real,
616  typename std::enable_if<std::is_arithmetic<Real>::value, int>::type = 0>
617  DiscreteGen(const std::initializer_list<Real>& il)
618  : DiscreteGen(il.begin(), il.end())
619  {
620  }
621 
622  template<typename Rng>
623  EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng)
624  {
625  using namespace Eigen::internal;
626  if (!cdf.empty())
627  {
628  auto rx = ur(std::forward<Rng>(rng));
629  return (_Scalar)(std::lower_bound(cdf.begin(), cdf.end() - 1, rx) - cdf.begin());
630  }
631  else
632  {
633  auto albit = pfirst(rng()) & alias_table.get_bitmask();
634  auto alx = ur(rng);
635  if (alx < alias_table.get_prob()[albit]) return albit;
636  return alias_table.get_alias()[albit];
637  }
638  }
639 
640  template<typename Packet, typename Rng>
641  EIGEN_STRONG_INLINE const Packet packetOp(Rng&& rng)
642  {
643  using namespace Eigen::internal;
644  using DPacket = decltype(reinterpret_to_double(std::declval<Packet>()));
645  if (!cdf.empty())
646  {
647  auto ret = pset1<Packet>(cdf.size());
648 #ifdef EIGEN_VECTORIZE_AVX
649  auto rx = ur.template packetOp<Packet4d>(std::forward<Rng>(rng));
650  for (auto& p : cdf)
651  {
652  auto c = reinterpret_to_int(pcmplt(rx, pset1<decltype(rx)>(p)));
653  auto r = combine_low32(c);
654  ret = padd(ret, r);
655  }
656 #else
657  auto rx1 = ur.template packetOp<DPacket>(rng),
658  rx2 = ur.template packetOp<DPacket>(rng);
659  for (auto& p : cdf)
660  {
661  auto pp = pset1<decltype(rx1)>(p);
662  ret = padd(ret, combine_low32(reinterpret_to_int(pcmplt(rx1, pp)), reinterpret_to_int(pcmplt(rx2, pp))));
663  }
664 #endif
665  return ret;
666  }
667  else
668  {
669 #ifdef EIGEN_VECTORIZE_AVX
670  using RUtils = RawbitsMaker<Packet, Rng>;
671  auto albit = pand(RUtils{}.rawbits(rng), pset1<Packet>(alias_table.get_bitmask()));
672  auto c = reinterpret_to_int(pcmplt(ur.template packetOp<Packet4d>(rng), pgather(alias_table.get_prob(), _mm256_castsi128_si256(albit))));
673  return pblendv(combine_low32(c), albit, pgather(alias_table.get_alias(), albit));
674 #else
675  using RUtils = RawbitsMaker<Packet, Rng>;
676  auto albit = pand(RUtils{}.rawbits(rng), pset1<Packet>(alias_table.get_bitmask()));
677  auto c1 = reinterpret_to_int(pcmplt(ur.template packetOp<DPacket>(rng), pgather(alias_table.get_prob(), albit)));
678  auto c2 = reinterpret_to_int(pcmplt(ur.template packetOp<DPacket>(rng), pgather(alias_table.get_prob(), albit, true)));
679  return pblendv(combine_low32(c1, c2), albit, pgather(alias_table.get_alias(), albit));
680 #endif
681  }
682  }
683  };
684 
685  template<typename> class BinomialGen;
686 
692  template<typename _Scalar>
693  class PoissonGen : OptCacheStore, public GenBase<PoissonGen<_Scalar>, _Scalar>
694  {
695  friend BinomialGen<_Scalar>;
696  static_assert(std::is_same<_Scalar, int32_t>::value, "poisson needs integral types.");
697  int cache_rest_cnt = 0;
698  UniformRealGen<float> ur;
699 
700  protected:
701  double mean, ne_mean, sqrt_tmean, log_mean, g1;
702 
703  public:
704  using Scalar = _Scalar;
705 
711  PoissonGen(double _mean)
712  : mean{ _mean }, ne_mean{ std::exp(-_mean) }
713  {
714  sqrt_tmean = std::sqrt(2 * mean);
715  log_mean = std::log(mean);
716  g1 = mean * log_mean - std::lgamma(mean + 1);
717  }
718 
719  template<typename Rng>
720  EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng)
721  {
722  using namespace Eigen::internal;
723  if (mean < 12)
724  {
725  _Scalar res = 0;
726  double val = 1;
727  for (; ; ++res)
728  {
729  val *= ur(rng);
730  if (val <= ne_mean) break;
731  }
732  return res;
733  }
734  else
735  {
736  _Scalar res;
737  double yx;
738  while (1)
739  {
740  yx = std::tan(constant::pi * ur(rng));
741  res = (_Scalar)(sqrt_tmean * yx + mean);
742  if (res >= 0 && ur(rng) <= 0.9 * (1.0 + yx * yx)
743  * std::exp(res * log_mean - std::lgamma(res + 1.0) - g1))
744  {
745  return res;
746  }
747  }
748  }
749  }
750 
751  template<typename Packet, typename Rng>
752  EIGEN_STRONG_INLINE const Packet packetOp(Rng&& rng)
753  {
754  using namespace Eigen::internal;
755  using PacketType = decltype(reinterpret_to_float(std::declval<Packet>()));
756 
757  if (mean < 12)
758  {
759  Packet res = pset1<Packet>(0);
760  PacketType val = pset1<PacketType>(1), pne_mean = pset1<PacketType>(ne_mean);
761  while (1)
762  {
763  val = pmul(val, ur.template packetOp<PacketType>(rng));
764  auto c = reinterpret_to_int(pcmplt(pne_mean, val));
765  if (pmovemask(c) == 0) break;
766  res = padd(res, pnegate(c));
767  }
768  return res;
769  }
770  else
771  {
772  auto& cm = Rand::detail::CompressMask<sizeof(Packet)>::get_inst();
773  const PacketType ppi = pset1<PacketType>(constant::pi),
774  psqrt_tmean = pset1<PacketType>(sqrt_tmean),
775  pmean = pset1<PacketType>(mean),
776  plog_mean = pset1<PacketType>(log_mean),
777  pg1 = pset1<PacketType>(g1);
778  while (1)
779  {
780  PacketType fres, yx, psin, pcos;
781  psincos(pmul(ppi, ur.template packetOp<PacketType>(rng)), psin, pcos);
782  yx = pdiv(psin, pcos);
783  fres = ptruncate(padd(pmul(psqrt_tmean, yx), pmean));
784 
785  auto p1 = pmul(padd(pmul(yx, yx), pset1<PacketType>(1)), pset1<PacketType>(0.9));
786  auto p2 = pexp(psub(psub(pmul(fres, plog_mean), plgamma(padd(fres, pset1<PacketType>(1)))), pg1));
787 
788  auto c1 = pcmple(pset1<PacketType>(0), fres);
789  auto c2 = pcmple(ur.template packetOp<PacketType>(rng), pmul(p1, p2));
790 
791  auto cands = fres;
792  bool full = false;
793  cache_rest_cnt = cm.compress_append(cands, pand(c1, c2),
794  OptCacheStore::template get<PacketType>(), cache_rest_cnt, full);
795  if (full) return pcast<PacketType, Packet>(cands);
796  }
797  }
798  }
799  };
800 
806  template<typename _Scalar>
807  class BinomialGen : public GenBase<BinomialGen<_Scalar>, _Scalar>
808  {
809  static_assert(std::is_same<_Scalar, int32_t>::value, "binomial needs integral types.");
810 
811  PoissonGen<_Scalar> poisson;
812  _Scalar trials;
813  double p, small_p, g1, sqrt_v, log_small_p, log_small_q;
814 
815  public:
816  using Scalar = _Scalar;
817 
824  BinomialGen(_Scalar _trials = 1, double _p = 0.5)
825  : poisson{ _trials * std::min(_p, 1 - _p) },
826  trials{ _trials }, p{ _p }, small_p{ std::min(p, 1 - p) }
827 
828  {
829  if (!(trials < 25 || poisson.mean < 1.0))
830  {
831  g1 = std::lgamma(trials + 1);
832  sqrt_v = std::sqrt(2 * poisson.mean * (1 - small_p));
833  log_small_p = std::log(small_p);
834  log_small_q = std::log(1 - small_p);
835  }
836  }
837 
838  template<typename Rng>
839  EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng)
840  {
841  using namespace Eigen::internal;
842  _Scalar res;
843  if (trials < 25)
844  {
845  res = 0;
846  for (int i = 0; i < trials; ++i)
847  {
848  if (poisson.ur(rng) < p) ++res;
849  }
850  return res;
851  }
852  else if (poisson.mean < 1.0)
853  {
854  res = poisson(rng);
855  }
856  else
857  {
858  while (1)
859  {
860  double ys;
861  ys = std::tan(constant::pi * poisson.ur(rng));
862  res = (_Scalar)(sqrt_v * ys + poisson.mean);
863  if (0 <= res && res <= trials && poisson.ur(rng) <= 1.2 * sqrt_v
864  * (1.0 + ys * ys)
865  * std::exp(g1 - std::lgamma(res + 1)
866  - std::lgamma(trials - res + 1.0)
867  + res * log_small_p
868  + (trials - res) * log_small_q)
869  )
870  {
871  break;
872  }
873  }
874  }
875  return p == small_p ? res : trials - res;
876  }
877 
878  template<typename Packet, typename Rng>
879  EIGEN_STRONG_INLINE const Packet packetOp(Rng&& rng)
880  {
881  using namespace Eigen::internal;
882  using PacketType = decltype(reinterpret_to_float(std::declval<Packet>()));
883  Packet res;
884  if (trials < 25)
885  {
886  PacketType pp = pset1<PacketType>(p);
887  res = pset1<Packet>(trials);
888  for (int i = 0; i < trials; ++i)
889  {
890  auto c = reinterpret_to_int(pcmple(pp, poisson.ur.template packetOp<PacketType>(rng)));
891  res = padd(res, c);
892  }
893  return res;
894  }
895  else if (poisson.mean < 1.0)
896  {
897  res = poisson.template packetOp<Packet>(rng);
898  }
899  else
900  {
901  auto& cm = Rand::detail::CompressMask<sizeof(Packet)>::get_inst();
902  const PacketType ppi = pset1<PacketType>(constant::pi),
903  ptrials = pset1<PacketType>(trials),
904  psqrt_v = pset1<PacketType>(sqrt_v),
905  pmean = pset1<PacketType>(poisson.mean),
906  plog_small_p = pset1<PacketType>(log_small_p),
907  plog_small_q = pset1<PacketType>(log_small_q),
908  pg1 = pset1<PacketType>(g1);
909  while (1)
910  {
911  PacketType fres, ys, psin, pcos;
912  psincos(pmul(ppi, poisson.ur.template packetOp<PacketType>(rng)), psin, pcos);
913  ys = pdiv(psin, pcos);
914  fres = ptruncate(padd(pmul(psqrt_v, ys), pmean));
915 
916  auto p1 = pmul(pmul(pset1<PacketType>(1.2), psqrt_v), padd(pset1<PacketType>(1), pmul(ys, ys)));
917  auto p2 = pexp(
918  padd(padd(psub(
919  psub(pg1, plgamma(padd(fres, pset1<PacketType>(1)))),
920  plgamma(psub(padd(ptrials, pset1<PacketType>(1)), fres))
921  ), pmul(fres, plog_small_p)), pmul(psub(ptrials, fres), plog_small_q))
922  );
923 
924  auto c1 = pand(pcmple(pset1<PacketType>(0), fres), pcmple(fres, ptrials));
925  auto c2 = pcmple(poisson.ur.template packetOp<PacketType>(rng), pmul(p1, p2));
926 
927  auto cands = fres;
928  bool full = false;
929  poisson.cache_rest_cnt = cm.compress_append(cands, pand(c1, c2),
930  poisson.template get<PacketType>(), poisson.cache_rest_cnt, full);
931  if (full)
932  {
933  res = pcast<PacketType, Packet>(cands);
934  break;
935  }
936  }
937  }
938  return p == small_p ? res : psub(pset1<Packet>(trials), res);
939  }
940  };
941 
947  template<typename _Scalar>
948  class GeometricGen : public GenBase<GeometricGen<_Scalar>, _Scalar>
949  {
950  static_assert(std::is_same<_Scalar, int32_t>::value, "geomtric needs integral types.");
951  UniformRealGen<float> ur;
952  double p, rlog_q;
953 
954  public:
955  using Scalar = _Scalar;
956 
962  GeometricGen(double _p)
963  : p{ _p }, rlog_q{ 1 / std::log(1 - p) }
964  {
965  }
966 
967  template<typename Rng>
968  EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng)
969  {
970  using namespace Eigen::internal;
971  return (_Scalar)(std::log(1 - ur(std::forward<Rng>(rng))) * rlog_q);
972  }
973 
974  template<typename Packet, typename Rng>
975  EIGEN_STRONG_INLINE const Packet packetOp(Rng&& rng)
976  {
977  using namespace Eigen::internal;
978  using PacketType = decltype(reinterpret_to_float(std::declval<Packet>()));
979 
980  return pcast<PacketType, Packet>(ptruncate(pmul(plog(
981  psub(pset1<PacketType>(1), ur.template packetOp<PacketType>(std::forward<Rng>(rng)))
982  ), pset1<PacketType>(rlog_q))));
983  }
984  };
985 
986 
987  template<typename Derived, typename Urng>
988  using UniformIntType = CwiseNullaryOp<internal::scalar_rng_adaptor<UniformIntGen<typename Derived::Scalar>, typename Derived::Scalar, Urng, true>, const Derived>;
989 
1003  template<typename Derived, typename Urng>
1004  inline const UniformIntType<Derived, Urng>
1005  uniformInt(Index rows, Index cols, Urng&& urng, typename Derived::Scalar min, typename Derived::Scalar max)
1006  {
1007  return {
1008  rows, cols, { std::forward<Urng>(urng), UniformIntGen<typename Derived::Scalar>{ min, max } }
1009  };
1010  }
1011 
1024  template<typename Derived, typename Urng>
1025  inline const UniformIntType<Derived, Urng>
1026  uniformIntLike(Derived& o, Urng&& urng, typename Derived::Scalar min, typename Derived::Scalar max)
1027  {
1028  return {
1029  o.rows(), o.cols(), { std::forward<Urng>(urng), UniformIntGen<typename Derived::Scalar>{ min, max } }
1030  };
1031  }
1032 
1033  template<typename Derived, typename Urng>
1034  using DiscreteFType = CwiseNullaryOp<internal::scalar_rng_adaptor<DiscreteGen<typename Derived::Scalar, float>, typename Derived::Scalar, Urng, true>, const Derived>;
1035 
1050  template<typename Derived, typename Urng, typename RealIter>
1051  inline const DiscreteFType<Derived, Urng>
1052  discreteF(Index rows, Index cols, Urng&& urng, RealIter first, RealIter last)
1053  {
1054  return {
1055  rows, cols, { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, float>{first, last} }
1056  };
1057  }
1058 
1072  template<typename Derived, typename Urng, typename RealIter>
1073  inline const DiscreteFType<Derived, Urng>
1074  discreteFLike(Derived& o, Urng&& urng, RealIter first, RealIter last)
1075  {
1076  return {
1077  o.rows(), o.cols(), { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, float>(first, last) }
1078  };
1079  }
1080 
1095  template<typename Derived, typename Urng, typename Real>
1096  inline const DiscreteFType<Derived, Urng>
1097  discreteF(Index rows, Index cols, Urng&& urng, const std::initializer_list<Real>& il)
1098  {
1099  return {
1100  rows, cols, { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, float>{il.begin(), il.end()} }
1101  };
1102  }
1103 
1117  template<typename Derived, typename Urng, typename Real>
1118  inline const DiscreteFType<Derived, Urng>
1119  discreteFLike(Derived& o, Urng&& urng, const std::initializer_list<Real>& il)
1120  {
1121  return {
1122  o.rows(), o.cols(), { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, float>(il.begin(), il.end()) }
1123  };
1124  }
1125 
1126  template<typename Derived, typename Urng>
1127  using DiscreteDType = CwiseNullaryOp<internal::scalar_rng_adaptor<DiscreteGen<typename Derived::Scalar, double>, typename Derived::Scalar, Urng, true>, const Derived>;
1128 
1143  template<typename Derived, typename Urng, typename RealIter>
1144  inline const DiscreteDType<Derived, Urng>
1145  discreteD(Index rows, Index cols, Urng&& urng, RealIter first, RealIter last)
1146  {
1147  return {
1148  rows, cols, { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, double>{first, last} }
1149  };
1150  }
1151 
1165  template<typename Derived, typename Urng, typename RealIter>
1166  inline const DiscreteDType<Derived, Urng>
1167  discreteDLike(Derived& o, Urng&& urng, RealIter first, RealIter last)
1168  {
1169  return {
1170  o.rows(), o.cols(), { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, double>{first, last} }
1171  };
1172  }
1173 
1188  template<typename Derived, typename Urng, typename Real>
1189  inline const DiscreteDType<Derived, Urng>
1190  discreteD(Index rows, Index cols, Urng&& urng, const std::initializer_list<Real>& il)
1191  {
1192  return {
1193  rows, cols, { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, double>{il.begin(), il.end()} }
1194  };
1195  }
1196 
1210  template<typename Derived, typename Urng, typename Real>
1211  inline const DiscreteDType<Derived, Urng>
1212  discreteDLike(Derived& o, Urng&& urng, const std::initializer_list<Real>& il)
1213  {
1214  return {
1215  o.rows(), o.cols(), { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, double>{il.begin(), il.end()} }
1216  };
1217  }
1218 
1219  template<typename Derived, typename Urng>
1220  using DiscreteType = CwiseNullaryOp<internal::scalar_rng_adaptor<DiscreteGen<typename Derived::Scalar, int32_t>, typename Derived::Scalar, Urng, true>, const Derived>;
1221 
1236  template<typename Derived, typename Urng, typename RealIter>
1237  inline const DiscreteType<Derived, Urng>
1238  discrete(Index rows, Index cols, Urng&& urng, RealIter first, RealIter last)
1239  {
1240  return {
1241  rows, cols, { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, int32_t>{first, last} }
1242  };
1243  }
1244 
1258  template<typename Derived, typename Urng, typename RealIter>
1259  inline const DiscreteType<Derived, Urng>
1260  discreteLike(Derived& o, Urng&& urng, RealIter first, RealIter last)
1261  {
1262  return {
1263  o.rows(), o.cols(), { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, int32_t>{first, last} }
1264  };
1265  }
1266 
1281  template<typename Derived, typename Urng, typename Real>
1282  inline const DiscreteType<Derived, Urng>
1283  discrete(Index rows, Index cols, Urng&& urng, const std::initializer_list<Real>& il)
1284  {
1285  return {
1286  rows, cols, { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, int32_t>{il.begin(), il.end()} }
1287  };
1288  }
1289 
1303  template<typename Derived, typename Urng, typename Real>
1304  inline const DiscreteType<Derived, Urng>
1305  discreteLike(Derived& o, Urng&& urng, const std::initializer_list<Real>& il)
1306  {
1307  return {
1308  o.rows(), o.cols(), { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, int32_t>{il.begin(), il.end()} }
1309  };
1310  }
1311 
1312  template<typename Derived, typename Urng>
1313  using PoissonType = CwiseNullaryOp<internal::scalar_rng_adaptor<PoissonGen<typename Derived::Scalar>, typename Derived::Scalar, Urng, true>, const Derived>;
1314 
1328  template<typename Derived, typename Urng>
1329  inline const PoissonType<Derived, Urng>
1330  poisson(Index rows, Index cols, Urng&& urng, double mean = 1)
1331  {
1332  return {
1333  rows, cols, { std::forward<Urng>(urng), PoissonGen<typename Derived::Scalar>{mean} }
1334  };
1335  }
1336 
1349  template<typename Derived, typename Urng>
1350  inline const PoissonType<Derived, Urng>
1351  poissonLike(Derived& o, Urng&& urng, double mean = 1)
1352  {
1353  return {
1354  o.rows(), o.cols(), { std::forward<Urng>(urng), PoissonGen<typename Derived::Scalar>{mean} }
1355  };
1356  }
1357 
1358  template<typename Derived, typename Urng>
1359  using BinomialType = CwiseNullaryOp<internal::scalar_rng_adaptor<BinomialGen<typename Derived::Scalar>, typename Derived::Scalar, Urng, true>, const Derived>;
1360 
1375  template<typename Derived, typename Urng>
1376  inline const BinomialType<Derived, Urng>
1377  binomial(Index rows, Index cols, Urng&& urng, typename Derived::Scalar trials = 1, double p = 0.5)
1378  {
1379  return {
1380  rows, cols, { std::forward<Urng>(urng), BinomialGen<typename Derived::Scalar>{trials, p} }
1381  };
1382  }
1383 
1397  template<typename Derived, typename Urng>
1398  inline const BinomialType<Derived, Urng>
1399  binomialLike(Derived& o, Urng&& urng, typename Derived::Scalar trials = 1, double p = 0.5)
1400  {
1401  return {
1402  o.rows(), o.cols(), { std::forward<Urng>(urng), BinomialGen<typename Derived::Scalar>{trials, p} }
1403  };
1404  }
1405 
1406  template<typename Derived, typename Urng>
1407  using GeometricType = CwiseNullaryOp<internal::scalar_rng_adaptor<GeometricGen<typename Derived::Scalar>, typename Derived::Scalar, Urng, true>, const Derived>;
1408 
1422  template<typename Derived, typename Urng>
1423  inline const GeometricType<Derived, Urng>
1424  geometric(Index rows, Index cols, Urng&& urng, double p = 0.5)
1425  {
1426  return {
1427  rows, cols, { std::forward<Urng>(urng), GeometricGen<typename Derived::Scalar>{p} }
1428  };
1429  }
1430 
1443  template<typename Derived, typename Urng>
1444  inline const GeometricType<Derived, Urng>
1445  geometricLike(Derived& o, Urng&& urng, double p = 0.5)
1446  {
1447  return {
1448  o.rows(), o.cols(), { std::forward<Urng>(urng), GeometricGen<typename Derived::Scalar>{p} }
1449  };
1450  }
1451  }
1452 }
1453 
1454 #endif
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:1424
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:1445
Eigen::Rand::GenBase
Base class of all univariate random generators.
Definition: Basic.h:33
Eigen::Rand::UniformIntGen::UniformIntGen
UniformIntGen(_Scalar _min, _Scalar _max)
Construct a new UniformInt Generator.
Definition: Discrete.h:231
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:1074
Eigen::Rand::GeometricGen::GeometricGen
GeometricGen(double _p)
Construct a new Geometric Generator.
Definition: Discrete.h:962
Eigen::Rand::DiscreteGen< _Scalar, int32_t >::DiscreteGen
DiscreteGen(const std::initializer_list< Real > &il)
Construct a new Discrete Generator.
Definition: Discrete.h:387
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:1330
Eigen::Rand::DiscreteGen< _Scalar, double >::DiscreteGen
DiscreteGen(RealIter first, RealIter last)
Construct a new Discrete Generator.
Definition: Discrete.h:584
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:1377
Eigen::Rand::DiscreteGen< _Scalar, int32_t >::DiscreteGen
DiscreteGen(RealIter first, RealIter last)
Construct a new Discrete Generator.
Definition: Discrete.h:354
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:1145
Eigen::Rand::poissonLike
const PoissonType< Derived, Urng > poissonLike(Derived &o, Urng &&urng, double mean=1)
generates reals on the Poisson distribution.
Definition: Discrete.h:1351
Eigen::Rand::PoissonGen
Generator of integers on a Poisson distribution.
Definition: Discrete.h:694
Eigen::Rand::DiscreteGen< _Scalar, float >::DiscreteGen
DiscreteGen(RealIter first, RealIter last)
Construct a new Discrete Generator.
Definition: Discrete.h:478
Eigen::Rand::GeometricGen
Generator of integers on a geometric distribution.
Definition: Discrete.h:949
Eigen::Rand::DiscreteGen
Generator of integers on the interval [0, n), where the probability of each individual integer i is p...
Definition: Discrete.h:324
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:1238
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:1399
Eigen::Rand::UniformRealGen< float >
Eigen::Rand::PoissonGen::PoissonGen
PoissonGen(double _mean)
Construct a new Poisson Generator.
Definition: Discrete.h:711
Eigen::Rand::BinomialGen
Generator of integers on a binomial distribution.
Definition: Discrete.h:808
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:1260
Eigen::Rand::UniformIntGen
Generator of integers with a given range [min, max]
Definition: Discrete.h:216
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:1005
Eigen::Rand::DiscreteGen< _Scalar, double >::DiscreteGen
DiscreteGen(const std::initializer_list< Real > &il)
Construct a new Discrete Generator.
Definition: Discrete.h:617
Eigen::Rand::BinomialGen::BinomialGen
BinomialGen(_Scalar _trials=1, double _p=0.5)
Construct a new Binomial Generator.
Definition: Discrete.h:824
Eigen::Rand::RandbitsGen
Generator of random bits for integral scalars.
Definition: Basic.h:248
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:1052
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:1026
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:1167
Eigen::Rand::DiscreteGen< _Scalar, float >::DiscreteGen
DiscreteGen(const std::initializer_list< Real > &il)
Construct a new Discrete Generator.
Definition: Discrete.h:511