On Sat, Feb 6, 2021 at 2:32 AM <camel-...@protonmail.com> wrote:
> I tried to implement a different implementation of the ziggurat method for
> generating standard normal distributions that is about twice as fast and
> uses 2/3 of the memory than the old one.
> I tested the implementation separately and am very confident it's correct,
> but it does fail 28 test in coverage testing.
> Checking the testing code I found out that all the failed tests are inside
> TestRandomDist which has the goal of "Make[ing] sure the random
> distribution returns the correct value for a given seed". Why would this be
> needed?
> The only explanation I can come up with is that it's standard_normal is,
> in regards to seeding, required to be backwards compatible. If that's the
> case how would, could one even implement a new algorithm?
>
Just for fun, I've attached the (C++) implementation that uses the
exponentially extended base block. Note that the constructor produces the
table.
Chuck
#include <cmath>
#include <stdint.h>
#include <cstdlib>
template <typename RNG>
class RandomNormal {
int32_t m_k[256];
double m_w[256];
double m_f[256];
RNG m_rng;
public:
explicit RandomNormal(RNG & rng) : m_rng(rng)
{
const double a = 4.928760963486943786e-03;
const double two31 = 2147483648.0;
double r = 3.655420419026941481;
m_f[0] = exp(-.5*r*r);
m_w[0] = a/(m_f[0]*two31);
m_k[0] = static_cast<long>(r/m_w[0]);
for(int i = 1; i < 255; ++i) {
const double y = m_f[i-1] + a/r;
const double x = sqrt(-2*log(y));
m_f[i] = y;
m_w[i] = r/two31;
m_k[i] = static_cast<int32_t>(x/m_w[i]);
r = x;
}
m_f[255] = 1.0;
m_w[255] = r/two31;
m_k[255] = 0;
}
double operator()()
{
while (1) {
const uint32_t randu = m_rng();
const uint32_t mask1 = static_cast<uint32_t>(0x000000ffUL);
const uint32_t mask2 = static_cast<uint32_t>(0xffffff00UL);
const int32_t j = static_cast<int32_t>(randu & mask2);
const int32_t i = static_cast<int32_t>(randu & mask1);
const double x = j*m_w[i];
if (abs(j) < m_k[i]) {
return x;
}
if (i > 0) {
const double randd = m_rng.double3();
const double y = m_f[i-1] + randd*(m_f[i] - m_f[i-1]);
if (y < exp(-.5*x*x)) {
return x;
}
}
else {
const double r = 3.655420419026941481;
const double t = -log(m_rng.double3())*(1./r);
const double y = -log(m_rng.double3());
if (2*y > t*t) {
return x > 0 ? (r + t) : -(r + t);
}
}
}
}
};
_______________________________________________
NumPy-Discussion mailing list
NumPy-Discussion@python.org
https://mail.python.org/mailman/listinfo/numpy-discussion
