I was definitely expecting a smoother bell curve so either my algorithm implementation is incorrect or the random numbers are bad (I tried to use a better random routine with 49 bytes of ASM code but that didn't make much difference and I have to call it at least twice per pair of samples since it only generates 8 bits at a time so it ended up slower I think).
My algorithm is supposed to be an implementation of this:
https://en.wikipedia.org/wiki/Box%E2%80 ... _transform
Based on the C++ implementation
Code: Select all
#include <cmath>
#include <limits>
#include <random>
#include <utility>
//"mu" is the mean of the distribution, and "sigma" is the standard deviation.
std::pair<double, double> generateGaussianNoise(double mu, double sigma)
{
constexpr double epsilon = std::numeric_limits<double>::epsilon();
constexpr double two_pi = 2.0 * M_PI;
//initialize the random uniform number generator (runif) in a range 0 to 1
static std::mt19937 rng(std::random_device{}()); // Standard mersenne_twister_engine seeded with rd()
static std::uniform_real_distribution<> runif(0.0, 1.0);
//create two random numbers, make sure u1 is greater than epsilon
double u1, u2;
do
{
u1 = runif(rng);
}
while (u1 <= epsilon);
u2 = runif(rng);
//compute z0 and z1
auto mag = sigma * sqrt(-2.0 * log(u1));
auto z0 = mag * cos(two_pi * u2) + mu;
auto z1 = mag * sin(two_pi * u2) + mu;
return std::make_pair(z0, z1);
}
I did check that the epsilon is good enough so that log(epsilon) does not cause a domain error, not sure if that can happen since I think the range for speccy rand is [0, 65535/65536]?
I should cache 2*pi*u2 since it is used twice as well.
Wondering why it is so jaggedy though I was expecting a nice smooth curve :/
I chose sigma = 128 (middle of screen) and standard deviation so that edge of the screen was +/- 3 std deviations away.
This is one project that is much easier in BASIC than ASM cos of the floating point maths and LN and trig.
I also have a Poisson distribution sampler (it's slow though and breaks - because it is adding teeny weeny floats to a larger float so loses accuracy and may never return an answer - which is not good - if you need loads of iterations).
I have a C# version of an exponential distribution sampler as well (which I actually used in production code for weather [rainfall] simulation).
Obviously a Normal distribution sampler is a super handy piece of code, last time I used one in production code (to give reaction times etc. for AI) it used a big lookup table, so the B-M transform is a big improvement on that, glad I found it (even though it didn't actually get used).
I also converted a golf game (to PS1) which used a Normal distribution to add random noise to CPU player shots (otherwise they got a hole in one every par 3 lol). That also used a lookup table approach.