Ziggurat algorithm for Gaussian random number generation.

Ziggurat.ahk

;* RandomNormal(Mean, Deviation)
;* Description:
	;* https://en.wikipedia.org/wiki/Ziggurat_algorithm
RandomNormal(vMean := 0, vDeviation := 1.0) {
	Static __K := (v := __RandomNormal()).k, __W := v.w, __F := v.f  ;* Populate the lookup tables.

	Loop {
		u := RandomUniform(-0x80000000, 0x7FFFFFFF), i := (u & 0xFF) + 1

		If (Abs(u) < __K[i]) {  ;* Rectangle. This will be the case for 99.33% of values (512 rectangles would be 99.64%).
			Return, (u*__W[i]*vDeviation + vMean)
		}

		x := u*__W[i]

		If (i == 1) {  ;* Base segment. Sample using a ratio of uniforms.
			While (2*y <= x**2) {
				x := -Ln(RandomUniform())*.27366123732975828, y := -Ln(RandomUniform())  ;? .27366123732975828 = 1/r
			}

			Return, (((u > 0)*2 - 1)*(3.6541528853610088 + x)*vDeviation + vMean)
		}

		If ((__F[i - 1] - __F[i])*RandomUniform() + __F[i] < Exp(-.5*x**2)) {  ;* Wedge.
			Return, (x*vDeviation + vMean)
		}

		;* The wedge was missed; start again.
	}
}
__RandomNormal() {
	r := 3.6541528853610088, v := 0.00492867323399, q := Exp(-.5*r**2)  ;? r = start of the tail, v = area of each rectangle
		, k := [(r*(q/v*2147483648.0)), 0], w := [(v/q)/2147483648.0], w[256] := r/2147483648.0, f := [1.0], f[256] := q  ;* Index zero is for the base segment, where Marsaglia and Tsang define this as k[0] = 2^31*r*f(r)/v, w[0] = .5^31*v/f(r), f[0] = 1.0.

	i := 256
	While (--i > 1) {
		x := Sqrt(-2.0*Ln(v/r + f[i + 1]))

		k[i + 1] := ((x/r)*2147483648.0), w[i] := x/2147483648.0, f[i] := Exp(-.5*x**2)

		r := x
	}

	Return, ({"k": k, "w": w, "f": f})
}

RandomUniform(vMin := 0, vMax := 1.0) {
	Random, r, vMin, vMax

	Return, (r)
}