skyscraper8/GUIs/skyscraper8.UI.ImGui.MonoGame.Bridge/Ziggurat.cs

using System.Diagnostics;
using System.Runtime.CompilerServices;

namespace skyscraper8.UI.MonoGame.Bridge;

//This one is heavily inspired from https://github.com/colgreen/Redzen/blob/main/Redzen/Numerics/Distributions/Double/ZigguratGaussian.cs
//To understand how this algorithm works, please read those Github comments.


public class Ziggurat
{
    const int __blockCount = 128;
    const double __R = 3.442619855899;
    const double __A = 9.91256303526217e-3;
    const ulong __MAXINT = (1UL << 53) - 1;
    const double __INCR = 1.0 / __MAXINT;
    private double[] __x;
    private double[] __y;
    private ulong[] __xComp;
    private double __A_Div_Y0;
    private Random rng;

    public Ziggurat()
        : this(null)
    {
        rng = new Random((int)DateTime.Now.Ticks);
    }
    
    public Ziggurat(int seed)
    : this(null)
    {
        rng = new Random(seed);
    }
    public Ziggurat(Random randomSrc)
    {
        this.rng = randomSrc;
        // Initialise rectangle position data.
        // __x[i] and __y[i] describe the top-right position of Box i.

        // Allocate storage. We add one to the length of _x so that we have an entry at __x[__blockCount], this avoids having
        // to do a special case test when sampling from the top box.
        __x = new double[__blockCount + 1];
        __y = new double[__blockCount];

        // Determine top right position of the base rectangle/box (the rectangle with the Gaussian tale attached).
        // We call this Box 0 or B0 for short.
        // Note. x[0] also describes the right-hand edge of B1. (See diagram).
        __x[0] = __R;
        __y[0] = GaussianPdfDenorm(__R);

        // The next box (B1) has a right hand X edge the same as B0.
        // Note. B1's height is the box area divided by its width, hence B1 has a smaller height than B0 because
        // B0's total area includes the attached distribution tail.
        __x[1] = __R;
        __y[1] = __y[0] + (__A / __x[1]);

        // Calc positions of all remaining rectangles.
        for(int i=2; i < __blockCount; i++)
        {
            __x[i] = GaussianPdfDenormInv(__y[i-1]);
            __y[i] = __y[i-1] + (__A / __x[i]);
        }

        // For completeness we define the right-hand edge of a notional box 6 as being zero (a box with no area).
        __x[__blockCount] = 0.0;

        // Useful precomputed values.
        __A_Div_Y0 = __A / __y[0];
        __xComp = new ulong[__blockCount];

        // Special case for base box. __xComp[0] stores the area of B0 as a proportion of __R
        // (recalling that all segments have area __A, but that the base segment is the combination of B0 and the distribution tail).
        // Thus __xComp[0] is the probability that a sample point is within the box part of the segment.
        __xComp[0] = (ulong)(((__R * __y[0]) / __A) * (double)__MAXINT);

        for(int i=1; i < __blockCount-1; i++)
        {
            __xComp[i] = (ulong)((__x[i+1] / __x[i]) * (double)__MAXINT);
        }
        __xComp[__blockCount-1] = 0;  // Shown for completeness.

        // Sanity check. Test that the top edge of the topmost rectangle is at y=1.0.
        // Note. We expect there to be a tiny drift away from 1.0 due to the inexactness of floating
        // point arithmetic.
        Debug.Assert(System.Math.Abs(1.0 - __y[__blockCount-1]) < 1e-10);
    }
    
    private static double GaussianPdfDenorm(double x)
    {
        return System.Math.Exp(-(x * x * 0.5));
    }
    
    private static double GaussianPdfDenormInv(double y)
    {
        // Operates over the y interval (0,1], which happens to be the y interval of the pdf,
        // with the exception that it does not include y=0, but we would never call with
        // y=0 so it doesn't matter. Note that a Gaussian effectively has a tail going
        // into infinity on the x-axis, hence asking what is x when y=0 is an invalid question
        // in the context of this class.
        return System.Math.Sqrt(-2.0 * System.Math.Log(y));
    }
    
    private void SampleTail(out double x)
    {
        double y;
        do
        {
            // Note. we use NextDoubleNonZero() because Log(0) returns -Infinity and will also tend to be a very slow execution path (when it occurs, which is rarely).
            x = -System.Math.Log(rng.NextDouble()) / __R;
            y = -System.Math.Log(rng.NextDouble());
        }
        while(y + y < x * x);
        x += __R;
    }
    
    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    private static void SetSignBit(ref double x, ref ulong signBit)
    {
        Unsafe.As<double, ulong>(ref x) |= signBit;
    }
    
    public double Sample()
    {
        double x = 0.0;
        while(true)
        {
            // Generate 64 random bits.
            ulong u = (ulong)rng.NextInt64();

            // Note. 61 random bits are required and therefore the lowest three bits are discarded
            // (a typical characteristic of PRNGs is that the least significant bits exhibit lower
            // quality randomness than the higher bits).
            // Note. The bit index values in comments are zero based, i..e the first bit is bit 0.
            // Select a segment (7 bits, bits 3 to 9).
            int s = (int)((u >> 3) & 0x7f);

            // Select the sign bit (bit 10), and shift it to the sign bit position for IEEE754
            // double-precision floating-point format.
            // Previously, the sign bit handling used a conditional branch that optionally multiplied the
            // positive Gaussian sample by -1. However, that approach is considerably slower because modern
            // superscalar CPUs rely heavily on branch prediction, but the sign bit here is randomly generated
            // and therefore entirely unpredictable, i.e. the absolute worse case scenario for a branch
            // predictor!
            ulong signBit = (u & 0x400UL) << 53;

            // Get a uniform random value with interval [0, 2^53-1], or in hexadecimal [0, 0x1f_ffff_ffff_ffff]
            // (i.e. a random 53 bit number) (bits 11 to 63).
            ulong u2 = u >> 11;

            // Special case for the base segment.
            if(s == 0)
            {
                if(u2 < __xComp[0])
                {
                    // Generated x is within R0.
                    x = u2 * __INCR * __A_Div_Y0;
                }
                else
                {
                    // Generated x is in the tail of the distribution.
                    SampleTail(out x);
                }

                SetSignBit(ref x, ref signBit);
                return x;
            }

            // All other segments.
            if(u2 < __xComp[s])
            {
                // Generated x is within the rectangle.
                x = u2 * __INCR * __x[s];
                SetSignBit(ref x, ref signBit);
                return x;
            }

            // Generated x is outside of the rectangle.
            // Generate a random y coordinate and test if our (x,y) is within the distribution curve.
            // This execution path is relatively slow/expensive (makes a call to Math.Exp()) but is relatively rarely executed,
            // although more often than the 'tail' path (above).
            x = u2 * __INCR * __x[s];
            if(__y[s-1] + ((__y[s] - __y[s-1]) * rng.NextDouble()) < GaussianPdfDenorm(x))
            {
                SetSignBit(ref x, ref signBit);
                return x;
            }
        }
    }
}