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(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; } } } }