Skip to content

Instantly share code, notes, and snippets.

@thynson
Last active March 23, 2023 05:24
Show Gist options
  • Save thynson/b25b8fe8f706958eaf21620280957def to your computer and use it in GitHub Desktop.
Save thynson/b25b8fe8f706958eaf21620280957def to your computer and use it in GitHub Desktop.
The Javascript(Typescript) port of Ziggurat method of generating a Gauss distributed variable.
interface RandomEngine {
randomUint32(): number;
}
const PARAM_R = 3.44428647676;
/* tabulated values for the heigt of the Ziggurat levels */
const ytab = [
1, 0.963598623011, 0.936280813353, 0.913041104253, 0.892278506696, 0.873239356919, 0.855496407634, 0.838778928349,
0.822902083699, 0.807732738234, 0.793171045519, 0.779139726505, 0.765577436082, 0.752434456248, 0.739669787677,
0.727249120285, 0.715143377413, 0.703327646455, 0.691780377035, 0.68048276891, 0.669418297233, 0.65857233912,
0.647931876189, 0.637485254896, 0.62722199145, 0.617132611532, 0.607208517467, 0.597441877296, 0.587825531465,
0.578352913803, 0.569017984198, 0.559815170911, 0.550739320877, 0.541785656682, 0.532949739145, 0.524227434628,
0.515614886373, 0.507108489253, 0.498704867478, 0.490400854812, 0.482193476986, 0.47407993601, 0.466057596125,
0.458123971214, 0.450276713467, 0.442513603171, 0.434832539473, 0.427231532022, 0.419708693379, 0.41226223212,
0.404890446548, 0.397591718955, 0.390364510382, 0.383207355816, 0.376118859788, 0.369097692334, 0.362142585282,
0.355252328834, 0.348425768415, 0.341661801776, 0.334959376311, 0.328317486588, 0.321735172063, 0.31521151497,
0.308745638367, 0.302336704338, 0.29598391232, 0.289686497571, 0.283443729739, 0.27725491156, 0.271119377649,
0.265036493387, 0.259005653912, 0.253026283183, 0.247097833139, 0.241219782932, 0.235391638239, 0.229612930649,
0.223883217122, 0.218202079518, 0.212569124201, 0.206983981709, 0.201446306496, 0.195955776745, 0.190512094256,
0.185114984406, 0.179764196185, 0.174459502324, 0.169200699492, 0.1639876086, 0.158820075195, 0.153697969964,
0.148621189348, 0.143589656295, 0.138603321143, 0.133662162669, 0.128766189309, 0.123915440582, 0.119109988745,
0.114349940703, 0.10963544023, 0.104966670533, 0.100343857232, 0.0957672718266, 0.0912372357329, 0.0867541250127,
0.082318375932, 0.0779304915295, 0.0735910494266, 0.0693007111742, 0.065060233529, 0.0608704821745, 0.056732448584,
0.05264727098, 0.0486162607163, 0.0446409359769, 0.0407230655415, 0.0368647267386, 0.0330683839378, 0.0293369977411,
0.0256741818288, 0.0220844372634, 0.0185735200577, 0.0151490552854, 0.0118216532614, 0.00860719483079,
0.00553245272614, 0.00265435214565,
];
/* tabulated values for 2^24 times x[i]/x[i+1],
* used to accept for U*x[i+1]<=x[i] without any floating point operations */
const ktab = [
0, 12590644, 14272653, 14988939, 15384584, 15635009, 15807561, 15933577, 16029594, 16105155, 16166147, 16216399,
16258508, 16294295, 16325078, 16351831, 16375291, 16396026, 16414479, 16431002, 16445880, 16459343, 16471578,
16482744, 16492970, 16502368, 16511031, 16519039, 16526459, 16533352, 16539769, 16545755, 16551348, 16556584,
16561493, 16566101, 16570433, 16574511, 16578353, 16581977, 16585398, 16588629, 16591685, 16594575, 16597311,
16599901, 16602354, 16604679, 16606881, 16608968, 16610945, 16612818, 16614592, 16616272, 16617861, 16619363,
16620782, 16622121, 16623383, 16624570, 16625685, 16626730, 16627708, 16628619, 16629465, 16630248, 16630969,
16631628, 16632228, 16632768, 16633248, 16633671, 16634034, 16634340, 16634586, 16634774, 16634903, 16634972,
16634980, 16634926, 16634810, 16634628, 16634381, 16634066, 16633680, 16633222, 16632688, 16632075, 16631380,
16630598, 16629726, 16628757, 16627686, 16626507, 16625212, 16623794, 16622243, 16620548, 16618698, 16616679,
16614476, 16612071, 16609444, 16606571, 16603425, 16599973, 16596178, 16591995, 16587369, 16582237, 16576520,
16570120, 16562917, 16554758, 16545450, 16534739, 16522287, 16507638, 16490152, 16468907, 16442518, 16408804,
16364095, 16301683, 16207738, 16047994, 15704248, 15472926,
];
/* tabulated values of 2^{-24}*x[i] */
const wtab = [
1.62318314817e-8, 2.16291505214e-8, 2.54246305087e-8, 2.84579525938e-8, 3.10340022482e-8, 3.33011726243e-8,
3.53439060345e-8, 3.72152672658e-8, 3.8950989572e-8, 4.05763964764e-8, 4.21101548915e-8, 4.35664624904e-8,
4.49563968336e-8, 4.62887864029e-8, 4.75707945735e-8, 4.88083237257e-8, 5.00063025384e-8, 5.11688950428e-8,
5.22996558616e-8, 5.34016475624e-8, 5.44775307871e-8, 5.55296344581e-8, 5.65600111659e-8, 5.75704813695e-8,
5.85626690412e-8, 5.95380306862e-8, 6.04978791776e-8, 6.14434034901e-8, 6.23756851626e-8, 6.32957121259e-8,
6.42043903937e-8, 6.51025540077e-8, 6.59909735447e-8, 6.68703634341e-8, 6.77413882848e-8, 6.8604668381e-8,
6.94607844804e-8, 7.03102820203e-8, 7.11536748229e-8, 7.1991448372e-8, 7.2824062723e-8, 7.36519550992e-8,
7.44755422158e-8, 7.52952223703e-8, 7.61113773308e-8, 7.69243740467e-8, 7.77345662086e-8, 7.85422956743e-8,
7.93478937793e-8, 8.01516825471e-8, 8.09539758128e-8, 8.17550802699e-8, 8.25552964535e-8, 8.33549196661e-8,
8.41542408569e-8, 8.49535474601e-8, 8.57531242006e-8, 8.65532538723e-8, 8.73542180955e-8, 8.8156298059e-8,
8.89597752521e-8, 8.97649321908e-8, 9.05720531451e-8, 9.138142487e-8, 9.21933373471e-8, 9.30080845407e-8,
9.38259651738e-8, 9.46472835298e-8, 9.54723502847e-8, 9.63014833769e-8, 9.71350089201e-8, 9.79732621669e-8,
9.88165885297e-8, 9.96653446693e-8, 1.00519899658e-7, 1.0138063623e-7, 1.02247952126e-7, 1.03122261554e-7,
1.04003996769e-7, 1.04893609795e-7, 1.05791574313e-7, 1.06698387725e-7, 1.07614573423e-7, 1.08540683296e-7,
1.09477300508e-7, 1.1042504257e-7, 1.11384564771e-7, 1.12356564007e-7, 1.13341783071e-7, 1.14341015475e-7,
1.15355110887e-7, 1.16384981291e-7, 1.17431607977e-7, 1.18496049514e-7, 1.19579450872e-7, 1.20683053909e-7,
1.21808209468e-7, 1.2295639141e-7, 1.24129212952e-7, 1.25328445797e-7, 1.26556042658e-7, 1.27814163916e-7,
1.29105209375e-7, 1.30431856341e-7, 1.31797105598e-7, 1.3320433736e-7, 1.34657379914e-7, 1.36160594606e-7,
1.37718982103e-7, 1.39338316679e-7, 1.41025317971e-7, 1.42787873535e-7, 1.44635331499e-7, 1.4657889173e-7,
1.48632138436e-7, 1.50811780719e-7, 1.53138707402e-7, 1.55639532047e-7, 1.58348931426e-7, 1.61313325908e-7,
1.64596952856e-7, 1.68292495203e-7, 1.72541128694e-7, 1.77574279496e-7, 1.83813550477e-7, 1.92166040885e-7,
2.05295471952e-7, 2.22600839893e-7,
];
// double
// gsl_ran_gaussian_ziggurat (gsl_rng *r, double sigma)
// {
// unsigned long U, sign, i, j;
// double x, y;
//
// while (1) {
// U = gsl_rng_uint32 (r);
// i = U & 0x0000007F; /* 7 bit to choose the step */
// sign = U & 0x00000080; /* 1 bit for the sign */
// j = U>>8; /* 24 bit for the x-value */
//
// x = j*wtab[i];
// if (j < ktab[i]) break;
//
// if (i<127) {
// double y0, y1;
// y0 = ytab[i];
// y1 = ytab[i+1];
// y = y1+(y0-y1)*gsl_rng_uniform(r);
// } else {
// x = PARAM_R - log(1.0-gsl_rng_uniform(r))/PARAM_R;
// y = exp(-PARAM_R*(x-0.5*PARAM_R))*gsl_rng_uniform(r);
// }
// if (y < exp(-0.5*x*x)) break;
// }
// return sign ? sigma*x : -sigma*x;
// }
export function uniformDistributed(randomEngine: RandomEngine) {
return randomEngine.randomUint32() * 2.3283064365386963e-10; // 2.3283064365386963e-10 is 2**-32
}
export function normalDistributed(randomEngine: RandomEngine, sigma: number): number {
let x;
let y;
let sign;
for (;;) {
const u = randomEngine.randomUint32();
const i = u & 0x7f;
sign = u & 0x80;
const j = u >>> 8;
x = j * wtab[i]!;
if (j < ktab[i]!) break;
const e = Math.exp(-0.5 * x * x);
if (i < 127) {
const y0 = ytab[i]!;
const y1 = ytab[i + 1]!;
y = y1 + (y0 - y1) * uniformDistributed(randomEngine);
} else {
x = PARAM_R - Math.log(1.0 - uniformDistributed(randomEngine)) / PARAM_R;
y = Math.exp(-PARAM_R * (x - 0.5 * PARAM_R)) * uniformDistributed(randomEngine);
}
if (y < e) {
break;
}
}
return sign ? sigma * x : -sigma * x;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment