EigenRand  0.5.0
 
Loading...
Searching...
No Matches
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
19namespace 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 + pmin;
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 StdUniformRealGen<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 StdUniformRealGen<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 StdUniformRealGen<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
1000 template<typename _Scalar>
1001 class BinomialVGen
1002 {
1003 static_assert(std::is_same<_Scalar, int32_t>::value, "binomial needs integral types.");
1004 };
1005
1006 template<>
1007 class BinomialVGen<int32_t> : public BinaryGenBase<BinomialVGen<int32_t>, int32_t, int32_t, float>
1008 {
1009 StdUniformRealGen<float> ur;
1010
1011 public:
1012 using Scalar = int32_t;
1013
1014 template<typename Rng>
1015 EIGEN_STRONG_INLINE const Scalar operator() (Rng&& rng, Scalar trials, float p)
1016 {
1017 using namespace Eigen::internal;
1018
1019 float small_p = std::min(p, 1 - p);
1020 float mean = trials * small_p;
1021 float g1 = std::lgamma(trials + 1);
1022 float sqrt_v = std::sqrt(2 * mean * (1 - small_p));
1023 float log_small_p = std::log(small_p);
1024 float log_small_q = std::log(1 - small_p);
1025
1026 Scalar res = 0;
1027 if (trials < 25)
1028 {
1029 res = 0;
1030 for (int i = 0; i < trials; ++i)
1031 {
1032 if (ur(rng) < p) ++res;
1033 }
1034 return res;
1035 }
1036 else if (mean < 1.)
1037 {
1038 float ne_mean = std::exp(-mean);
1039 float val = 1;
1040 for (; ; ++res)
1041 {
1042 val *= ur(rng);
1043 if (val <= ne_mean) break;
1044 }
1045 }
1046 else
1047 {
1048 for (int _i = 0; ; ++_i)
1049 {
1050 EIGENRAND_CHECK_INFINITY_LOOP();
1051 double ys;
1052 ys = std::tan(constant::pi * ur(rng));
1053 res = (Scalar)(sqrt_v * ys + mean);
1054 if (0 <= res && res <= trials && ur(rng) <= 1.2 * sqrt_v
1055 * (1.0 + ys * ys)
1056 * std::exp(g1 - std::lgamma(res + 1)
1057 - std::lgamma(trials - res + 1.0)
1058 + res * log_small_p
1059 + (trials - res) * log_small_q)
1060 )
1061 {
1062 break;
1063 }
1064 }
1065 }
1066 return p == small_p ? res : trials - res;
1067 }
1068
1069 template<typename Packet, typename FPacket, typename Rng>
1070 EIGEN_STRONG_INLINE const Packet packetOp(Rng&& rng, const Packet& trials, const FPacket& p)
1071 {
1072 using namespace Eigen::internal;
1073 Packet valid = pset1<Packet>(0);
1074 FPacket valid_res;
1075
1076 const auto psmall_p = pmin(p, psub(pset1<FPacket>(1), p));
1077 const auto ptrials = pcast<Packet, FPacket>(trials);
1078 const auto pmean = pmul(ptrials, psmall_p);
1079 const auto pg1 = plgamma_approx(padd(ptrials, pset1<FPacket>(1)));
1080 const auto psqrt_v = psqrt(pmul(pmul(pset1<FPacket>(2), pmean), psub(pset1<FPacket>(1), psmall_p)));
1081 const auto plog_small_p = plog(psmall_p);
1082 const auto plog_small_q = plog(psub(pset1<FPacket>(1), psmall_p));
1083 const auto ppi = pset1<FPacket>(constant::pi);
1084 valid = reinterpret_to_int(pcmplt(pmean, pset1<FPacket>(1)));
1085 if (predux_any(valid))
1086 {
1087 Packet res = pset1<Packet>(0);
1088 FPacket val = pset1<FPacket>(1), pne_mean = pexp(pnegate(pmean));
1089 for (int _i = 0; ; ++_i)
1090 {
1091 EIGENRAND_CHECK_INFINITY_LOOP();
1092 val = pmul(val, ur.template packetOp<FPacket>(rng));
1093 auto c = pand(reinterpret_to_int(pcmplt(pne_mean, val)), valid);
1094 if (!predux_any(c)) break;
1095 res = padd(res, pnegate(c));
1096 }
1097 valid_res = pcast<Packet, FPacket>(res);
1098 }
1099
1100 if (!predux_all(valid))
1101 {
1102 for (int _i = 0; ; ++_i)
1103 {
1104 EIGENRAND_CHECK_INFINITY_LOOP();
1105 FPacket fres, ys, psin, pcos;
1106 psincos(pmul(ppi, ur.template packetOp<FPacket>(rng)), psin, pcos);
1107 ys = pdiv(psin, pcos);
1108 fres = ptruncate(padd(pmul(psqrt_v, ys), pmean));
1109
1110 auto p1 = pmul(pmul(pset1<FPacket>(1.2), psqrt_v), padd(pset1<FPacket>(1), pmul(ys, ys)));
1111 auto p2 = pexp(
1112 padd(padd(psub(
1113 psub(pg1, plgamma_approx(padd(fres, pset1<FPacket>(1)))),
1114 plgamma_approx(psub(padd(ptrials, pset1<FPacket>(1)), fres))
1115 ), pmul(fres, plog_small_p)), pmul(psub(ptrials, fres), plog_small_q))
1116 );
1117
1118 auto c1 = pand(pcmple(pset1<FPacket>(0), fres), pcmple(fres, ptrials));
1119 auto c2 = pcmple(ur.template packetOp<FPacket>(rng), pmul(p1, p2));
1120
1121 auto fvalid = pand(c1, c2);
1122 valid_res = pblendv(pnew_andnot(fvalid, reinterpret_to_float(valid)), fres, valid_res);
1123 valid = por(reinterpret_to_int(fvalid), valid);
1124
1125 if (predux_all(valid))
1126 {
1127 break;
1128 }
1129 }
1130 }
1131
1132 return pblendv(
1133 reinterpret_to_int(pcmpeq(p, psmall_p)),
1134 pcast<FPacket, Packet>(valid_res),
1135 psub(trials, pcast<FPacket, Packet>(valid_res))
1136 );
1137 }
1138 };
1139
1140 template<>
1141 class BinomialVGen<float> : public BinaryGenBase<BinomialVGen<float>, float, float, float>
1142 {
1143 StdUniformRealGen<float> ur;
1144
1145 public:
1146 using Scalar = float;
1147
1148 template<typename Rng>
1149 EIGEN_STRONG_INLINE const Scalar operator() (Rng&& rng, Scalar _trials, float p)
1150 {
1151 using namespace Eigen::internal;
1152
1153 auto trials = reinterpret_cast<int32_t&>(_trials);
1154
1155 float small_p = std::min(p, 1 - p);
1156 float mean = trials * small_p;
1157 float g1 = std::lgamma(trials + 1);
1158 float sqrt_v = std::sqrt(2 * mean * (1 - small_p));
1159 float log_small_p = std::log(small_p);
1160 float log_small_q = std::log(1 - small_p);
1161
1162 int32_t res = 0;
1163 if (trials < 25)
1164 {
1165 res = 0;
1166 for (int i = 0; i < trials; ++i)
1167 {
1168 if (ur(rng) < p) ++res;
1169 }
1170 return res;
1171 }
1172 else if (mean < 1.)
1173 {
1174 float ne_mean = std::exp(-mean);
1175 float val = 1;
1176 for (; ; ++res)
1177 {
1178 val *= ur(rng);
1179 if (val <= ne_mean) break;
1180 }
1181 }
1182 else
1183 {
1184 for (int _i = 0; ; ++_i)
1185 {
1186 EIGENRAND_CHECK_INFINITY_LOOP();
1187 double ys;
1188 ys = std::tan(constant::pi * ur(rng));
1189 res = (int32_t)(sqrt_v * ys + mean);
1190 if (0 <= res && res <= trials && ur(rng) <= 1.2 * sqrt_v
1191 * (1.0 + ys * ys)
1192 * std::exp(g1 - std::lgamma(res + 1)
1193 - std::lgamma(trials - res + 1.0)
1194 + res * log_small_p
1195 + (trials - res) * log_small_q)
1196 )
1197 {
1198 break;
1199 }
1200 }
1201 }
1202 res = p == small_p ? res : trials - res;
1203 return reinterpret_cast<Scalar&>(res);
1204 }
1205
1206 template<typename Packet, typename Rng>
1207 EIGEN_STRONG_INLINE const Packet packetOp(Rng&& rng, const Packet& _trials, const Packet& p)
1208 {
1209 using namespace Eigen::internal;
1210 using IPacket = decltype(reinterpret_to_int(std::declval<Packet>()));
1211 IPacket valid = pset1<IPacket>(0);
1212 Packet valid_res;
1213
1214 const auto trials = reinterpret_to_int(_trials);
1215 const auto psmall_p = pmin(p, psub(pset1<Packet>(1), p));
1216 const auto ptrials = pcast<IPacket, Packet>(trials);
1217 const auto pmean = pmul(ptrials, psmall_p);
1218 const auto pg1 = plgamma_approx(padd(ptrials, pset1<Packet>(1)));
1219 const auto psqrt_v = psqrt(pmul(pmul(pset1<Packet>(2), pmean), psub(pset1<Packet>(1), psmall_p)));
1220 const auto plog_small_p = plog(psmall_p);
1221 const auto plog_small_q = plog(psub(pset1<Packet>(1), psmall_p));
1222 const auto ppi = pset1<Packet>(constant::pi);
1223 valid = reinterpret_to_int(pcmplt(pmean, pset1<Packet>(1)));
1224 if (predux_any(reinterpret_to_float(valid)))
1225 {
1226 Packet res = pset1<Packet>(0);
1227 Packet val = pset1<Packet>(1), pne_mean = pexp(pnegate(pmean));
1228 for (int _i = 0; ; ++_i)
1229 {
1230 EIGENRAND_CHECK_INFINITY_LOOP();
1231 val = pmul(val, ur.template packetOp<Packet>(rng));
1232 auto c = pand(reinterpret_to_int(pcmplt(pne_mean, val)), valid);
1233 if (!predux_any(reinterpret_to_float(c))) break;
1234 res = padd(res, pcast<IPacket, Packet>(pnegate(c)));
1235 }
1236 valid_res = res;
1237 }
1238
1239 if (!predux_all(reinterpret_to_float(valid)))
1240 {
1241 for (int _i = 0; ; ++_i)
1242 {
1243 EIGENRAND_CHECK_INFINITY_LOOP();
1244 Packet fres, ys, psin, pcos;
1245 psincos(pmul(ppi, ur.template packetOp<Packet>(rng)), psin, pcos);
1246 ys = pdiv(psin, pcos);
1247 fres = ptruncate(padd(pmul(psqrt_v, ys), pmean));
1248
1249 auto p1 = pmul(pmul(pset1<Packet>(1.2), psqrt_v), padd(pset1<Packet>(1), pmul(ys, ys)));
1250 auto p2 = pexp(
1251 padd(padd(psub(
1252 psub(pg1, plgamma_approx(padd(fres, pset1<Packet>(1)))),
1253 plgamma_approx(psub(padd(ptrials, pset1<Packet>(1)), fres))
1254 ), pmul(fres, plog_small_p)), pmul(psub(ptrials, fres), plog_small_q))
1255 );
1256
1257 auto c1 = pand(pcmple(pset1<Packet>(0), fres), pcmple(fres, ptrials));
1258 auto c2 = pcmple(ur.template packetOp<Packet>(rng), pmul(p1, p2));
1259
1260 auto fvalid = pand(c1, c2);
1261 valid_res = pblendv(pnew_andnot(fvalid, reinterpret_to_float(valid)), fres, valid_res);
1262 valid = por(reinterpret_to_int(fvalid), valid);
1263
1264 if (predux_all(reinterpret_to_float(valid)))
1265 {
1266 break;
1267 }
1268 }
1269 }
1270 return pblendv(
1271 pcmpeq(p, psmall_p),
1272 reinterpret_to_float(pcast<Packet, IPacket>(valid_res)),
1273 reinterpret_to_float(psub(trials, pcast<Packet, IPacket>(valid_res)))
1274 );
1275 }
1276 };
1277
1283 template<typename _Scalar>
1284 class GeometricGen : public GenBase<GeometricGen<_Scalar>, _Scalar>
1285 {
1286 static_assert(std::is_same<_Scalar, int32_t>::value, "geomtric needs integral types.");
1287 StdUniformRealGen<float> ur;
1288 double p, rlog_q;
1289
1290 public:
1291 using Scalar = _Scalar;
1292
1298 GeometricGen(double _p = 0.5)
1299 : p{ _p }, rlog_q{ 1 / std::log(1 - p) }
1300 {
1301 }
1302
1303 GeometricGen(const GeometricGen&) = default;
1304 GeometricGen(GeometricGen&&) = default;
1305
1306 GeometricGen& operator=(const GeometricGen&) = default;
1307 GeometricGen& operator=(GeometricGen&&) = default;
1308
1309 template<typename Rng>
1310 EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng)
1311 {
1312 using namespace Eigen::internal;
1313 return (_Scalar)(std::log(1 - ur(std::forward<Rng>(rng))) * rlog_q);
1314 }
1315
1316 template<typename Packet, typename Rng>
1317 EIGEN_STRONG_INLINE const Packet packetOp(Rng&& rng)
1318 {
1319 using namespace Eigen::internal;
1320 using PacketType = decltype(reinterpret_to_float(std::declval<Packet>()));
1321
1322 return pcast<PacketType, Packet>(ptruncate(pmul(plog(
1323 psub(pset1<PacketType>(1), ur.template packetOp<PacketType>(std::forward<Rng>(rng)))
1324 ), pset1<PacketType>(rlog_q))));
1325 }
1326 };
1327
1328
1329 template<typename Derived, typename Urng>
1330 using UniformIntType = CwiseNullaryOp<internal::scalar_rng_adaptor<UniformIntGen<typename Derived::Scalar>, typename Derived::Scalar, Urng, true>, const Derived>;
1331
1345 template<typename Derived, typename Urng>
1346 inline const UniformIntType<Derived, Urng>
1347 uniformInt(Index rows, Index cols, Urng&& urng, typename Derived::Scalar min, typename Derived::Scalar max)
1348 {
1349 return {
1350 rows, cols, { std::forward<Urng>(urng), UniformIntGen<typename Derived::Scalar>{ min, max } }
1351 };
1352 }
1353
1366 template<typename Derived, typename Urng>
1367 inline const UniformIntType<Derived, Urng>
1368 uniformIntLike(Derived& o, Urng&& urng, typename Derived::Scalar min, typename Derived::Scalar max)
1369 {
1370 return {
1371 o.rows(), o.cols(), { std::forward<Urng>(urng), UniformIntGen<typename Derived::Scalar>{ min, max } }
1372 };
1373 }
1374
1375 template<typename Derived, typename Urng>
1376 using DiscreteFType = CwiseNullaryOp<internal::scalar_rng_adaptor<DiscreteGen<typename Derived::Scalar, float>, typename Derived::Scalar, Urng, true>, const Derived>;
1377
1392 template<typename Derived, typename Urng, typename RealIter>
1393 inline const DiscreteFType<Derived, Urng>
1394 discreteF(Index rows, Index cols, Urng&& urng, RealIter first, RealIter last)
1395 {
1396 return {
1397 rows, cols, { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, float>{first, last} }
1398 };
1399 }
1400
1414 template<typename Derived, typename Urng, typename RealIter>
1415 inline const DiscreteFType<Derived, Urng>
1416 discreteFLike(Derived& o, Urng&& urng, RealIter first, RealIter last)
1417 {
1418 return {
1419 o.rows(), o.cols(), { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, float>(first, last) }
1420 };
1421 }
1422
1437 template<typename Derived, typename Urng, typename Real>
1438 inline const DiscreteFType<Derived, Urng>
1439 discreteF(Index rows, Index cols, Urng&& urng, const std::initializer_list<Real>& il)
1440 {
1441 return {
1442 rows, cols, { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, float>{il.begin(), il.end()} }
1443 };
1444 }
1445
1459 template<typename Derived, typename Urng, typename Real>
1460 inline const DiscreteFType<Derived, Urng>
1461 discreteFLike(Derived& o, Urng&& urng, const std::initializer_list<Real>& il)
1462 {
1463 return {
1464 o.rows(), o.cols(), { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, float>(il.begin(), il.end()) }
1465 };
1466 }
1467
1468 template<typename Derived, typename Urng>
1469 using DiscreteDType = CwiseNullaryOp<internal::scalar_rng_adaptor<DiscreteGen<typename Derived::Scalar, double>, typename Derived::Scalar, Urng, true>, const Derived>;
1470
1485 template<typename Derived, typename Urng, typename RealIter>
1486 inline const DiscreteDType<Derived, Urng>
1487 discreteD(Index rows, Index cols, Urng&& urng, RealIter first, RealIter last)
1488 {
1489 return {
1490 rows, cols, { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, double>{first, last} }
1491 };
1492 }
1493
1507 template<typename Derived, typename Urng, typename RealIter>
1508 inline const DiscreteDType<Derived, Urng>
1509 discreteDLike(Derived& o, Urng&& urng, RealIter first, RealIter last)
1510 {
1511 return {
1512 o.rows(), o.cols(), { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, double>{first, last} }
1513 };
1514 }
1515
1530 template<typename Derived, typename Urng, typename Real>
1531 inline const DiscreteDType<Derived, Urng>
1532 discreteD(Index rows, Index cols, Urng&& urng, const std::initializer_list<Real>& il)
1533 {
1534 return {
1535 rows, cols, { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, double>{il.begin(), il.end()} }
1536 };
1537 }
1538
1552 template<typename Derived, typename Urng, typename Real>
1553 inline const DiscreteDType<Derived, Urng>
1554 discreteDLike(Derived& o, Urng&& urng, const std::initializer_list<Real>& il)
1555 {
1556 return {
1557 o.rows(), o.cols(), { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, double>{il.begin(), il.end()} }
1558 };
1559 }
1560
1561 template<typename Derived, typename Urng>
1562 using DiscreteType = CwiseNullaryOp<internal::scalar_rng_adaptor<DiscreteGen<typename Derived::Scalar, int32_t>, typename Derived::Scalar, Urng, true>, const Derived>;
1563
1578 template<typename Derived, typename Urng, typename RealIter>
1579 inline const DiscreteType<Derived, Urng>
1580 discrete(Index rows, Index cols, Urng&& urng, RealIter first, RealIter last)
1581 {
1582 return {
1583 rows, cols, { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, int32_t>{first, last} }
1584 };
1585 }
1586
1600 template<typename Derived, typename Urng, typename RealIter>
1601 inline const DiscreteType<Derived, Urng>
1602 discreteLike(Derived& o, Urng&& urng, RealIter first, RealIter last)
1603 {
1604 return {
1605 o.rows(), o.cols(), { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, int32_t>{first, last} }
1606 };
1607 }
1608
1623 template<typename Derived, typename Urng, typename Real>
1624 inline const DiscreteType<Derived, Urng>
1625 discrete(Index rows, Index cols, Urng&& urng, const std::initializer_list<Real>& il)
1626 {
1627 return {
1628 rows, cols, { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, int32_t>{il.begin(), il.end()} }
1629 };
1630 }
1631
1645 template<typename Derived, typename Urng, typename Real>
1646 inline const DiscreteType<Derived, Urng>
1647 discreteLike(Derived& o, Urng&& urng, const std::initializer_list<Real>& il)
1648 {
1649 return {
1650 o.rows(), o.cols(), { std::forward<Urng>(urng), DiscreteGen<typename Derived::Scalar, int32_t>{il.begin(), il.end()} }
1651 };
1652 }
1653
1654 template<typename Derived, typename Urng>
1655 using PoissonType = CwiseNullaryOp<internal::scalar_rng_adaptor<PoissonGen<typename Derived::Scalar>, typename Derived::Scalar, Urng, true>, const Derived>;
1656
1670 template<typename Derived, typename Urng>
1671 inline const PoissonType<Derived, Urng>
1672 poisson(Index rows, Index cols, Urng&& urng, double mean = 1)
1673 {
1674 return {
1675 rows, cols, { std::forward<Urng>(urng), PoissonGen<typename Derived::Scalar>{mean} }
1676 };
1677 }
1678
1691 template<typename Derived, typename Urng>
1692 inline const PoissonType<Derived, Urng>
1693 poissonLike(Derived& o, Urng&& urng, double mean = 1)
1694 {
1695 return {
1696 o.rows(), o.cols(), { std::forward<Urng>(urng), PoissonGen<typename Derived::Scalar>{mean} }
1697 };
1698 }
1699
1700 template<typename Derived, typename Urng>
1701 using BinomialType = CwiseNullaryOp<internal::scalar_rng_adaptor<BinomialGen<typename Derived::Scalar>, typename Derived::Scalar, Urng, true>, const Derived>;
1702
1717 template<typename Derived, typename Urng>
1718 inline const BinomialType<Derived, Urng>
1719 binomial(Index rows, Index cols, Urng&& urng, typename Derived::Scalar trials = 1, double p = 0.5)
1720 {
1721 return {
1722 rows, cols, { std::forward<Urng>(urng), BinomialGen<typename Derived::Scalar>{trials, p} }
1723 };
1724 }
1725
1739 template<typename Derived, typename Urng>
1740 inline const BinomialType<Derived, Urng>
1741 binomialLike(Derived& o, Urng&& urng, typename Derived::Scalar trials = 1, double p = 0.5)
1742 {
1743 return {
1744 o.rows(), o.cols(), { std::forward<Urng>(urng), BinomialGen<typename Derived::Scalar>{trials, p} }
1745 };
1746 }
1747
1748
1749 template<typename Lhs, typename Rhs, typename Urng>
1750 using BinomialVVType = CwiseHeteroBinaryOp<
1751 internal::scalar_binary_rng_adaptor<BinomialVGen<typename Lhs::Scalar>, typename Lhs::Scalar, typename Lhs::Scalar, typename Rhs::Scalar, Urng, true>,
1752 const Lhs, const Rhs
1753 >;
1754
1768 template<typename Lhs, typename Rhs, typename Urng>
1769 inline const BinomialVVType<Lhs, Rhs, Urng>
1770 binomial(Urng&& urng, const ArrayBase<Lhs>& trials, const ArrayBase<Rhs>& p)
1771 {
1772 return {
1773 static_cast<const Lhs&>(trials), static_cast<const Rhs&>(p),
1774 { std::forward<Urng>(urng), BinomialVGen<typename Lhs::Scalar>{} }
1775 };
1776 }
1777
1778 namespace impl{
1779 template<class Derived, class NewType> struct CastType
1780 {
1781 using type = typename internal::eval<
1782 typename std::remove_reference<
1783 typename internal::cast_return_type<Derived, const CwiseUnaryOp<internal::scalar_cast_op<typename Derived::Scalar, NewType>, const Derived> >::type
1784 >::type
1785 >::type;
1786 };
1787 }
1788
1789 template<typename Lhs, typename Urng>
1790 using BinomialVSType = CwiseHeteroBinaryOp<
1791 internal::scalar_binary_rng_adaptor<BinomialVGen<typename Lhs::Scalar>, typename Lhs::Scalar, typename Lhs::Scalar, float, Urng, true>,
1792 const Lhs, CwiseNullaryOp<internal::scalar_constant_op<float>, const typename impl::CastType<Lhs, float>::type>
1793 >;
1794
1795 template<typename Lhs, typename Urng>
1796 inline const BinomialVSType<Lhs, Urng>
1797 binomial(Urng&& urng, const ArrayBase<Lhs>& trials, float p)
1798 {
1799 return {
1800 static_cast<const Lhs&>(trials),
1801 { trials.rows(), trials.cols(), internal::scalar_constant_op<float>{p} },
1802 { std::forward<Urng>(urng), BinomialVGen<typename Lhs::Scalar>{} }
1803 };
1804 }
1805
1806 namespace impl
1807 {
1808 template<typename Lhs, typename Urng>
1809 using BinomialVType = CwiseBinaryOp<
1810 internal::scalar_binary_rng_adaptor<BinomialVGen<typename Lhs::Scalar>, typename Lhs::Scalar, typename Lhs::Scalar, typename Lhs::Scalar, Urng, true>,
1811 const Lhs, CwiseNullaryOp<internal::scalar_constant_op<float>, const Lhs>
1812 >;
1813
1814 template<typename Lhs, typename Urng>
1815 inline const BinomialVType<Lhs, Urng>
1816 binomial(Urng&& urng, const ArrayBase<Lhs>& trials, float p)
1817 {
1818 static_assert(std::is_same<typename Lhs::Scalar, float>::value, "Lhs must have float scalar type.");
1819 return {
1820 static_cast<const Lhs&>(trials),
1821 { trials.rows(), trials.cols(), internal::scalar_constant_op<float>{p} },
1822 { std::forward<Urng>(urng), BinomialVGen<typename Lhs::Scalar>{} }
1823 };
1824 }
1825 }
1826
1827 template<typename Rhs, typename Urng>
1828 using BinomialSVType = CwiseHeteroBinaryOp<
1829 internal::scalar_binary_rng_adaptor<BinomialVGen<int32_t>, int32_t, int32_t, typename Rhs::Scalar, Urng, true>,
1830 CwiseNullaryOp<internal::scalar_constant_op<int32_t>, const typename impl::CastType<Rhs, int32_t>::type>, const Rhs
1831 >;
1832
1833 template<typename Rhs, typename Urng>
1834 inline const BinomialSVType<Rhs, Urng>
1835 binomial(Urng&& urng, int32_t trials, const ArrayBase<Rhs>& p)
1836 {
1837 return {
1838 { p.rows(), p.cols(), internal::scalar_constant_op<int32_t>{trials} },
1839 static_cast<const Rhs&>(p),
1840 { std::forward<Urng>(urng), BinomialVGen<int32_t>{} }
1841 };
1842 }
1843
1844 template<typename Derived, typename Urng>
1845 using GeometricType = CwiseNullaryOp<internal::scalar_rng_adaptor<GeometricGen<typename Derived::Scalar>, typename Derived::Scalar, Urng, true>, const Derived>;
1846
1860 template<typename Derived, typename Urng>
1861 inline const GeometricType<Derived, Urng>
1862 geometric(Index rows, Index cols, Urng&& urng, double p = 0.5)
1863 {
1864 return {
1865 rows, cols, { std::forward<Urng>(urng), GeometricGen<typename Derived::Scalar>{p} }
1866 };
1867 }
1868
1881 template<typename Derived, typename Urng>
1882 inline const GeometricType<Derived, Urng>
1883 geometricLike(Derived& o, Urng&& urng, double p = 0.5)
1884 {
1885 return {
1886 o.rows(), o.cols(), { std::forward<Urng>(urng), GeometricGen<typename Derived::Scalar>{p} }
1887 };
1888 }
1889 }
1890
1891#ifdef EIGEN_VECTORIZE_NEON
1892 namespace internal
1893 {
1894 template<typename _Scalar, typename Urng, bool _mutable>
1895 struct functor_traits<scalar_rng_adaptor<Rand::DiscreteGen<_Scalar, double>, _Scalar, Urng, _mutable> >
1896 {
1897 enum { Cost = HugeCost, PacketAccess = 0, IsRepeatable = false };
1898 };
1899 }
1900#endif
1901}
1902
1903#endif
Eigen::CwiseHeteroBinaryOp
Generic expression where a coefficient-wise binary operator is applied to two expressions.
Definition: CwiseHeteroBinaryOp.h:87
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:1285
Eigen::Rand::GeometricGen::GeometricGen
GeometricGen(double _p=0.5)
Construct a new Geometric Generator.
Definition: Discrete.h:1298
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:349
Eigen::Rand::StdUniformRealGen
Generator of reals in a range [0, 1)
Definition: Basic.h:495
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:1693
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:1416
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:1347
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:1883
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:1741
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:1602
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:1719
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:1862
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:1672
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:1368
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:1580
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:1394
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:1487
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:1509