12 #ifndef EIGENRAND_DISTS_GAMMAPOISSON_H
13 #define EIGENRAND_DISTS_GAMMAPOISSON_H
28 template<
typename _Scalar>
31 static_assert(std::is_same<_Scalar, int32_t>::value,
"negativeBinomial needs integral types.");
35 using Scalar = _Scalar;
44 :
gamma{ (float)_trials, (
float)((1 - _p) / _p) }
49 template<
typename Rng>
50 EIGEN_STRONG_INLINE
const _Scalar operator() (Rng&& rng)
52 using namespace Eigen::internal;
54 return PoissonGen<_Scalar>{v}(rng);
57 template<
typename Packet,
typename Rng>
58 EIGEN_STRONG_INLINE
const Packet packetOp(Rng&& rng)
60 using namespace Eigen::internal;
61 using ur_base = UniformRealGen<float>;
62 using PacketType = decltype(reinterpret_to_float(std::declval<Packet>()));
64 auto mean =
gamma.template packetOp<PacketType>(rng);
65 auto res = pset1<Packet>(0);
66 PacketType val = pset1<PacketType>(1), pne_mean = pexp(pnegate(mean));
67 if (pmovemask(pcmplt(pset1<PacketType>(12), mean)) == 0)
71 val = pmul(val, ur.template packetOp<PacketType>(rng));
72 auto c = reinterpret_to_int(pcmplt(pne_mean, val));
73 if (pmovemask(c) == 0)
break;
74 res = padd(res, pnegate(c));
80 auto& cm = Rand::detail::CompressMask<
sizeof(Packet)>::get_inst();
81 const PacketType ppi = pset1<PacketType>(constant::pi),
82 psqrt_tmean = psqrt(pmul(pset1<PacketType>(2), mean)),
83 plog_mean = plog(mean),
84 pg1 = psub(pmul(mean, plog_mean), plgamma(padd(mean, pset1<PacketType>(1))));
87 PacketType fres, yx, psin, pcos;
88 psincos(pmul(ppi, ur.template packetOp<PacketType>(rng)), psin, pcos);
89 yx = pdiv(psin, pcos);
90 fres = ptruncate(padd(pmul(psqrt_tmean, yx), mean));
92 auto p1 = pmul(padd(pmul(yx, yx), pset1<PacketType>(1)), pset1<PacketType>(0.9));
93 auto p2 = pexp(psub(psub(pmul(fres, plog_mean), plgamma(padd(fres, pset1<PacketType>(1)))), pg1));
95 auto c1 = pcmple(pset1<PacketType>(0), fres);
96 auto c2 = pcmple(ur.template packetOp<PacketType>(rng), pmul(p1, p2));
100 gamma.cache_rest_cnt = cm.compress_append(cands, pand(c1, c2),
101 gamma.template get<PacketType>(),
gamma.cache_rest_cnt, full);
102 if (full)
return pcast<PacketType, Packet>(cands);
108 template<
typename Derived,
typename Urng>
109 using NegativeBinomialType = CwiseNullaryOp<internal::scalar_rng_adaptor<NegativeBinomialGen<typename Derived::Scalar>,
typename Derived::Scalar, Urng,
true>,
const Derived>;
125 template<
typename Derived,
typename Urng>
126 inline const NegativeBinomialType<Derived, Urng>
127 negativeBinomial(Index rows, Index cols, Urng&& urng,
typename Derived::Scalar trials = 1,
double p = 0.5)
147 template<
typename Derived,
typename Urng>
148 inline const NegativeBinomialType<Derived, Urng>