| 90 | |
| 91 | template <typename T, typename U> |
| 92 | void fill_gamma(Xoroshiro128plus* rng, U* dst, size_t size, U* shape, U* scale) { |
| 93 | for (size_t i = 0; i < size; ++i) { |
| 94 | T a = static_cast<T>(shape[i]); |
| 95 | T b = static_cast<T>(scale[i]); |
| 96 | T scale = b; |
| 97 | bool a_less_one = a < 1.f ? true : false; |
| 98 | if (a <= 0) { |
| 99 | dst[i] = U(0.0f); |
| 100 | continue; |
| 101 | }; |
| 102 | T d = a + (a_less_one ? 2.0f / 3.0f : -1.0f / 3.0f); |
| 103 | T c = 1.0f / std::sqrt(9.0f * d); |
| 104 | while (true) { |
| 105 | T x, y; |
| 106 | x = normal_sample<T>(rng); |
| 107 | y = 1.0f + c * x; |
| 108 | if (y <= 0) |
| 109 | continue; |
| 110 | T v = y * y * y; |
| 111 | T u = uniform_sample<T>(rng); |
| 112 | T xx = x * x; |
| 113 | if ((u < 1.0f - 0.0331f * xx * xx) || |
| 114 | std::log(u) < 0.5f * xx + d * (1.0f - v + std::log(v))) { |
| 115 | dst[i] = U(scale * d * v); |
| 116 | if (a_less_one) |
| 117 | dst[i] *= U(std::pow(uniform_sample<T>(rng), T(1.f / a))); |
| 118 | break; |
| 119 | } |
| 120 | } |
| 121 | } |
| 122 | } |
| 123 | |
| 124 | template <typename U> |
| 125 | void fill_multinomial_without_replacement( |