util/random/normal.h


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38

#pragma once

#include <cmath>

// sometimes we need stateless normal distribution...

/*
 * normal distribution with Box-Muller transform
 * http://www.design.caltech.edu/erik/Misc/Gaussian.html
 */
template <typename T, typename TRng>
static inline T StdNormalDistribution(TRng&& rng) noexcept {
    T x;
    T y;
    T r;

    do {
        x = (T)rng.GenRandReal1() * T(2) - T(1);
        y = (T)rng.GenRandReal1() * T(2) - T(1);
        r = x * x + y * y;
    } while (r > T(1) || r <= T(0));

    return x * std::sqrt(-T(2) * std::log(r) / r);
}

template <typename T, typename TRng>
static inline T NormalDistribution(TRng&& rng, T m, T d) noexcept {
    return StdNormalDistribution<T>(rng) * d + m;
}

// specialized for float, double, long double
template <class T>
T StdNormalRandom() noexcept;

template <class T>
static inline T NormalRandom(T m, T d) noexcept {
    return StdNormalRandom<T>() * d + m;
}