34#include <Random123/philox.h>
74template <RNGSalt salt>
77 using rng_type = r123::Philox4x64;
78 using ctr_type = rng_type::ctr_type;
79 using key_type = rng_type::key_type;
81 const ctr_type c{{counter, 0
u, 0
u, 0
u}};
83 auto const id1 =
static_cast<uint32_t
>(key1);
84 auto const id2 =
static_cast<uint32_t
>(key2);
88 return rng_type{}(c, k);
110template <
RNGSalt salt, std::size_t
N = 3,
111 std::enable_if_t<(
N > 1) and (
N <= 4),
int> = 0>
112auto noise_uniform(uint64_t counter, uint32_t seed,
int key1,
int key2 = 0) {
114 auto const integers = philox_4_uint64s<salt>(counter, seed, key1, key2);
116 std::transform(integers.begin(), integers.begin() +
N, noise.begin(),
117 [](std::size_t value) { return Utils::uniform(value) - 0.5; });
121template <RNGSalt salt, std::
size_t N, std::enable_if_t<N == 1,
int> = 0>
122auto noise_uniform(uint64_t counter, uint32_t seed,
int key1,
int key2 = 0) {
124 auto const integers = philox_4_uint64s<salt>(counter, seed, key1, key2);
148template <
RNGSalt salt, std::size_t
N = 3,
149 class = std::enable_if_t<(
N >= 1) and (
N <= 4)>>
152 auto const integers = philox_4_uint64s<salt>(counter, seed, key1, key2);
153 static const double epsilon = std::numeric_limits<double>::min();
155 constexpr std::size_t M = (
N <= 2) ? 2 : 4;
157 std::transform(integers.begin(), integers.begin() + M,
u.begin(),
158 [](std::size_t value) {
159 auto u = Utils::uniform(value);
160 return (u < epsilon) ? epsilon : u;
168 constexpr double two_pi = 2.0 *
Utils::pi();
170 auto const modulo = sqrt(-2.0 * log(
u[0]));
171 auto const angle = two_pi *
u[1];
172 noise[0] = modulo * cos(angle);
174 noise[1] = modulo * sin(angle);
178 auto const modulo = sqrt(-2.0 * log(
u[2]));
179 auto const angle = two_pi *
u[3];
180 noise[2] = modulo * cos(angle);
182 noise[3] = modulo * sin(angle);
195template <
typename T> std::mt19937
mt19937(T &&seed) {
196 std::mt19937 generator(seed);
197 generator.discard(1'000'000);
Vector implementation and trait types for boost qvm interoperability.
auto philox_4_uint64s(uint64_t counter, uint32_t seed, int key1, int key2=0)
get 4 random uint 64 from the Philox RNG
auto noise_uniform(uint64_t counter, uint32_t seed, int key1, int key2=0)
Generator for random uniform noise.
auto noise_gaussian(uint64_t counter, uint32_t seed, int key1, int key2=0)
Generator for Gaussian noise.
std::mt19937 mt19937(T &&seed)
Mersenne Twister with warmup.
DEVICE_QUALIFIER constexpr T pi()
Ratio of diameter and circumference of a circle.
constexpr double uniform(uint64_t in)
Uniformly map unsigned integer to double.
constexpr uint64_t u32_to_u64(uint32_t high, uint32_t low)