1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192 |
- function Ziggurat(seeder = 123456789){
- var jsr = seeder;
- var wn = Array(128);
- var fn = Array(128);
- var kn = Array(128);
- function RNOR(){
- var hz = SHR3();
- var iz = hz & 127;
- return (Math.abs(hz) < kn[iz]) ? hz * wn[iz] : nfix(hz, iz);
- }
- this.nextGaussian = function(){
- return RNOR();
- };
- function nfix(hz, iz){
- var r = 3.442619855899;
- var r1 = 1.0 / r;
- var x;
- var y;
- while(true){
- x = hz * wn[iz];
- if( iz === 0 ){
- x = (-Math.log(UNI()) * r1);
- y = -Math.log(UNI());
- while( y + y < x * x){
- x = (-Math.log(UNI()) * r1);
- y = -Math.log(UNI());
- }
- return ( hz > 0 ) ? r+x : -r-x;
- }
- if( fn[iz] + UNI() * (fn[iz-1] - fn[iz]) < Math.exp(-0.5 * x * x) ){
- return x;
- }
- hz = SHR3();
- iz = hz & 127;
-
- if( Math.abs(hz) < kn[iz]){
- return (hz * wn[iz]);
- }
- }
- }
- function SHR3(){
- var jz = jsr;
- var jzr = jsr;
- jzr ^= (jzr << 13);
- jzr ^= (jzr >>> 17);
- jzr ^= (jzr << 5);
- jsr = jzr;
- return (jz+jzr) | 0;
- }
- function UNI(){
- return 0.5 * (1 + SHR3() / -Math.pow(2,31));
- }
- function zigset(){
- // seed generator based on current time
- //jsr ^= new Date().getTime();
- var m1 = 2147483648.0;
- var dn = 3.442619855899;
- var tn = dn;
- var vn = 9.91256303526217e-3;
-
- var q = vn / Math.exp(-0.5 * dn * dn);
- kn[0] = Math.floor((dn/q)*m1);
- kn[1] = 0;
- wn[0] = q / m1;
- wn[127] = dn / m1;
- fn[0] = 1.0;
- fn[127] = Math.exp(-0.5 * dn * dn);
- for(var i = 126; i >= 1; i--){
- dn = Math.sqrt(-2.0 * Math.log( vn / dn + Math.exp( -0.5 * dn * dn)));
- kn[i+1] = Math.floor((dn/tn)*m1);
- tn = dn;
- fn[i] = Math.exp(-0.5 * dn * dn);
- wn[i] = dn / m1;
- }
- }
- zigset();
- }
- module.exports = Ziggurat;
|