104 lines
2.5 KiB
C++
104 lines
2.5 KiB
C++
|
/*
|
|||
|
|
* This Source Code Form is subject to the terms of the Mozilla Public
|
|||
|
|
* License, v. 2.0. If a copy of the MPL was not distributed with this
|
|||
|
|
* file, You can obtain one at http://mozilla.org/MPL/2.0/.
|
|||
|
|
*
|
|||
|
|
* Copyright 2020, Danny Robson <danny@nerdcruft.net>
|
|||
|
|
*/
|
|||
|
|
|
|||
|
|
#pragma once
|
|||
|
|
|
|||
|
|
#include "uniform.hpp"
|
|||
|
|
|
|||
|
|
#include <type_traits>
|
|||
|
|
|
|||
|
|
|
|||
|
|
namespace cruft::rand::distribution {
|
|||
|
|
template <typename ResultT>
|
|||
|
|
requires (std::is_floating_point_v<ResultT>)
|
|||
|
|
class normal {
|
|||
|
|
public:
|
|||
|
|
using result_type = ResultT;
|
|||
|
|
struct param_type {
|
|||
|
|
result_type mean;
|
|||
|
|
result_type stddev;
|
|||
|
|
};
|
|||
|
|
|
|||
|
|
normal ():
|
|||
|
|
normal (0)
|
|||
|
|
{ ; }
|
|||
|
|
|
|||
|
|
explicit normal (result_type mean, result_type stddev = 1)
|
|||
|
|
: m_param {
|
|||
|
|
.mean = mean,
|
|||
|
|
.stddev = stddev
|
|||
|
|
}
|
|||
|
|
{ ; }
|
|||
|
|
|
|||
|
|
|
|||
|
|
explicit normal (param_type const &_param)
|
|||
|
|
: m_param (_param)
|
|||
|
|
{ ; }
|
|||
|
|
|
|||
|
|
|
|||
|
|
void reset (void)
|
|||
|
|
{
|
|||
|
|
m_live = false;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
template <typename GeneratorT>
|
|||
|
|
result_type
|
|||
|
|
operator() (GeneratorT &&g)
|
|||
|
|
{
|
|||
|
|
return (*this)(g, m_param);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
// We use the Box–Muller transform to convert pairs of uniform reals
|
|||
|
|
// to normally distributed reals.
|
|||
|
|
template <typename GeneratorT>
|
|||
|
|
result_type
|
|||
|
|
operator() (GeneratorT &&g, param_type const ¶ms)
|
|||
|
|
{
|
|||
|
|
if (m_live) {
|
|||
|
|
m_live = false;
|
|||
|
|
return m_prev * params.stddev + params.mean;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
auto [u, v, s] = find_uvs (g);
|
|||
|
|
result_type z0 = u * std::sqrt (-2 * std::log (s) / s);
|
|||
|
|
result_type z1 = v * std::sqrt (-2 * std::log (s) / s);
|
|||
|
|
|
|||
|
|
m_prev = z1;
|
|||
|
|
m_live = true;
|
|||
|
|
|
|||
|
|
return z0 * params.stddev + params.mean;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
private:
|
|||
|
|
template <typename GeneratorT>
|
|||
|
|
std::tuple<result_type, result_type, result_type>
|
|||
|
|
find_uvs (GeneratorT &&g)
|
|||
|
|
{
|
|||
|
|
while (1) {
|
|||
|
|
uniform_real_distribution<result_type> unit (-1, 1);
|
|||
|
|
result_type u = unit (g);
|
|||
|
|
result_type v = unit (g);
|
|||
|
|
result_type s = u * u + v * v;
|
|||
|
|
|
|||
|
|
#pragma GCC diagnostic push
|
|||
|
|
#pragma GCC diagnostic ignored "-Wfloat-equal"
|
|||
|
|
if (s != 0 && s < 1)
|
|||
|
|
return { u, v, s };
|
|||
|
|
#pragma GCC diagnostic pop
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
param_type m_param;
|
|||
|
|
bool m_live = false;
|
|||
|
|
result_type m_prev;
|
|||
|
|
};
|
|||
|
|
}
|