| 176 | |
| 177 | template <typename T, typename U> |
| 178 | void fill_poisson(Xoroshiro128plus* rng, U* dst, U* lam, size_t size) { |
| 179 | for (size_t i = 0; i < size; ++i) { |
| 180 | T lambda = static_cast<T>(lam[i]); |
| 181 | T exp_neg_lambda = std::exp(-lambda); |
| 182 | T log_lambda = std::log(lambda), sqrt_lambda = std::sqrt(lambda); |
| 183 | T b = 0.931f + 2.53f * sqrt_lambda; |
| 184 | T a = -0.059f + 0.02483f * b; |
| 185 | T inv_alpha = 1.1239f + 1.1328f / (b - 3.4f); |
| 186 | T vr = 0.9277f - 3.6224f / (b - 2.f); |
| 187 | T u, v, u_shifted, k; |
| 188 | if (lambda == 0) { |
| 189 | dst[i] = U(0); |
| 190 | continue; |
| 191 | } |
| 192 | if (lambda < 10) { |
| 193 | T prod = 1, x = 0; |
| 194 | u = 0; |
| 195 | while (true) { |
| 196 | u = uniform_sample<T>(rng); |
| 197 | prod *= u; |
| 198 | if (prod <= exp_neg_lambda) { |
| 199 | dst[i] = U(x); |
| 200 | break; |
| 201 | } |
| 202 | x += 1; |
| 203 | } |
| 204 | continue; |
| 205 | } |
| 206 | while (true) { |
| 207 | u = uniform_sample<T>(rng) - T(0.5f); |
| 208 | v = uniform_sample<T>(rng); |
| 209 | u_shifted = T(0.5f) - std::abs(u); |
| 210 | k = std::floor((T(2.f) * a / u_shifted + b) * u + lambda + T(0.43f)); |
| 211 | if (u_shifted >= 0.07 && v < vr) { |
| 212 | dst[i] = U(k); |
| 213 | break; |
| 214 | } |
| 215 | if (k < 0 || (u_shifted < T(0.013f) && v > u_shifted)) { |
| 216 | continue; |
| 217 | } |
| 218 | if ((std::log(v) + std::log(inv_alpha) - |
| 219 | std::log(a / (u_shifted * u_shifted) + b)) <= |
| 220 | (-lambda + k * log_lambda - std::lgamma(k + 1))) { |
| 221 | dst[i] = U(k); |
| 222 | break; |
| 223 | } |
| 224 | } |
| 225 | } |
| 226 | } |
| 227 | |
| 228 | template <typename T, typename U> |
| 229 | void fill_beta(Xoroshiro128plus* rng, U* dst, U* alpha, U* beta, size_t size) { |