randmtzig.c

/*
   A C-program for MT19937, with initialization improved 2002/2/10.
   Coded by Takuji Nishimura and Makoto Matsumoto.
   This is a faster version by taking Shawn Cokus's optimization,
   Matthe Bellew's simplification, Isaku Wada's real version.
   David Bateman added normal and exponential distributions following
   Marsaglia and Tang's Ziggurat algorithm.

   Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
   Copyright (C) 2004, David Bateman
   All rights reserved.

   Redistribution and use in source and binary forms, with or without
   modification, are permitted provided that the following conditions
   are met:

     1. Redistributions of source code must retain the above copyright
        notice, this list of conditions and the following disclaimer.

     2. Redistributions in binary form must reproduce the above copyright
        notice, this list of conditions and the following disclaimer in the
        documentation and/or other materials provided with the distribution.

     3. The names of its contributors may not be used to endorse or promote
        products derived from this software without specific prior written
        permission.

   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER
   OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.


   Any feedback is very welcome.
   http://www.math.keio.ac.jp/matumoto/emt.html
   email: matumoto@math.keio.ac.jp
*/

/*
  Modified by Viral B. Shah for julia to support dsfmt and only __LP64__
  systems. 52-bits of randomness are used from the mantissa of random double
  precision numbers generated by dsfmt.
 */

#include <math.h>
#include <stdio.h>
#include <stddef.h>
#include <time.h>
#ifndef _MSC_VER
#include <sys/time.h>
#endif

#include <stdlib.h>
#include <string.h>

#define DSFMT_DO_NOT_USE_OLD_NAMES
#define DSFMT_MEXP 19937
#include <dSFMT.h>

typedef ptrdiff_t randmtzig_idx_type;
typedef signed char randmtzig_int8_t;
typedef unsigned char randmtzig_uint8_t;
typedef short randmtzig_int16_t;
typedef unsigned short randmtzig_uint16_t;
typedef int randmtzig_int32_t;
typedef unsigned int randmtzig_uint32_t;
typedef long long randmtzig_int64_t;
typedef unsigned long long randmtzig_uint64_t;

/* Declarations */

extern double randmtzig_randn (dsfmt_t *dsfmt);
extern double randmtzig_gv_randn (void);
extern double randmtzig_gv_exprnd (void);

/* ===== Uniform generators ===== */

inline static randmtzig_uint64_t gv_randi (void)
{
    double r = dsfmt_gv_genrand_close1_open2();
    return *((uint64_t *) &r) & 0x000fffffffffffff;
}

/* generates a random number on (0,1) with 53-bit resolution */
inline static double gv_randu (void)
{
    return dsfmt_gv_genrand_open_open();
}

inline static randmtzig_uint64_t randi (dsfmt_t *dsfmt)
{
    double r = dsfmt_genrand_close1_open2(dsfmt);
    return *((uint64_t *) &r) & 0x000fffffffffffff;
}

/* generates a random number on (0,1) with 53-bit resolution */
inline static double randu (dsfmt_t *dsfmt)
{
    return dsfmt_genrand_open_open(dsfmt);
}

/* ===== Ziggurat normal and exponential generators ===== */
# define ZIGINT randmtzig_uint64_t
# define EMANTISSA 4503599627370496  /* 52 bit mantissa */
# define ERANDI gv_randi() /* 52 bits for mantissa */
# define NMANTISSA 2251799813685248
# define NRANDI gv_randi() /* 51 bits for mantissa + 1 bit sign */
# define RANDU gv_randu()

#define ZIGGURAT_TABLE_SIZE 256

#define ZIGGURAT_NOR_R 3.6541528853610088
#define ZIGGURAT_NOR_INV_R 0.27366123732975828
#define NOR_SECTION_AREA 0.00492867323399

#define ZIGGURAT_EXP_R 7.69711747013104972
#define ZIGGURAT_EXP_INV_R 0.129918765548341586
#define EXP_SECTION_AREA 0.0039496598225815571993


/*
This code is based on the paper Marsaglia and Tsang, "The ziggurat method
for generating random variables", Journ. Statistical Software. Code was
presented in this paper for a Ziggurat of 127 levels and using a 32 bit
integer random number generator. This version of the code, uses the
Mersenne Twister as the integer generator and uses 256 levels in the
Ziggurat. This has several advantages.

  1) As Marsaglia and Tsang themselves states, the more levels the few
     times the expensive tail algorithm must be called
  2) The cycle time of the generator is determined by the integer
     generator, thus the use of a Mersenne Twister for the core random
     generator makes this cycle extremely long.
  3) The license on the original code was unclear, thus rewriting the code
     from the article means we are free of copyright issues.
  4) Compile flag for full 53-bit random mantissa.

It should be stated that the authors made my life easier, by the fact that
the algorithm developed in the text of the article is for a 256 level
ziggurat, even if the code itself isn't...

One modification to the algorithm developed in the article, is that it is
assumed that 0 <= x < Inf, and "unsigned long"s are used, thus resulting in
terms like 2^32 in the code. As the normal distribution is defined between
-Inf < x < Inf, we effectively only have 31 bit integers plus a sign. Thus
in Marsaglia and Tsang, terms like 2^32 become 2^31. We use NMANTISSA for
this term.  The exponential distribution is one sided so we use the
full 32 bits.  We use EMANTISSA for this term.

It appears that I'm slightly slower than the code in the article, this
is partially due to a better generator of random integers than they
use. But might also be that the case of rapid return was optimized by
inlining the relevant code with a #define. As the basic Mersenne
Twister is only 25% faster than this code I suspect that the main
reason is just the use of the Mersenne Twister and not the inlining,
so I'm not going to try and optimize further.
*/


// void randmtzig_create_ziggurat_tables (void)
// {
//   int i;
//   double x, x1;

//   /* Ziggurat tables for the normal distribution */
//   x1 = ZIGGURAT_NOR_R;
//   wi[255] = x1 / NMANTISSA;
//   fi[255] = exp (-0.5 * x1 * x1);

//   /* Index zero is special for tail strip, where Marsaglia and Tsang
//    * defines this as
//    * k_0 = 2^31 * r * f(r) / v, w_0 = 0.5^31 * v / f(r), f_0 = 1,
//    * where v is the area of each strip of the ziggurat.
//    */
//   ki[0] = (ZIGINT) (x1 * fi[255] / NOR_SECTION_AREA * NMANTISSA);
//   wi[0] = NOR_SECTION_AREA / fi[255] / NMANTISSA;
//   fi[0] = 1.;

//   for (i = 254; i > 0; i--)
//     {
//       /* New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus
//        * need inverse operator of y = exp(-0.5*x*x) -> x = sqrt(-2*ln(y))
//        */
//       x = sqrt(-2. * log(NOR_SECTION_AREA / x1 + fi[i+1]));
//       ki[i+1] = (ZIGINT)(x / x1 * NMANTISSA);
//       wi[i] = x / NMANTISSA;
//       fi[i] = exp (-0.5 * x * x);
//       x1 = x;
//     }

//   ki[1] = 0;

//   /* Zigurrat tables for the exponential distribution */
//   x1 = ZIGGURAT_EXP_R;
//   we[255] = x1 / EMANTISSA;
//   fe[255] = exp (-x1);

//   /* Index zero is special for tail strip, where Marsaglia and Tsang
//    * defines this as
//    * k_0 = 2^32 * r * f(r) / v, w_0 = 0.5^32 * v / f(r), f_0 = 1,
//    * where v is the area of each strip of the ziggurat.
//    */
//   ke[0] = (ZIGINT) (x1 * fe[255] / EXP_SECTION_AREA * EMANTISSA);
//   we[0] = EXP_SECTION_AREA / fe[255] / EMANTISSA;
//   fe[0] = 1.;

//   for (i = 254; i > 0; i--)
//     {
//       /* New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus
//        * need inverse operator of y = exp(-x) -> x = -ln(y)
//        */
//       x = - log(EXP_SECTION_AREA / x1 + fe[i+1]);
//       ke[i+1] = (ZIGINT)(x / x1 * EMANTISSA);
//       we[i] = x / EMANTISSA;
//       fe[i] = exp (-x);
//       x1 = x;
//     }
//   ke[1] = 0;
// }

// Tables for randn
static ZIGINT ki[ZIGGURAT_TABLE_SIZE] =
    {2104047571230236,               0,1693657211688499,1919380038164751,
     2015384402142420,2068365869415708,2101878624030987,2124958784087614,
     2141808670783638,2154644611559370,2164744887580145,2172897953690771,
     2179616279367521,2185247251864556,2190034623104318,2194154434518163,
     2197736978772008,2200880740889623,2203661538008543,2206138681107245,
     2208359231804928,2210361007256700,2212174742387166,2213825672703393,
     2215334711001466,2216719334486539,2217994262138197,2219171977964129,
     2220263139537873,2221276900116549,2222221164932202,2223102796828387,
     2223927782546019,2224701368169460,2225428170203747,2226112267247709,
     2226757276104752,2227366415327922,2227942558554233,2228488279492093,
     2229005890046815,2229497472774805,2229964908626691,2230409900758245,
     2230833995044249,2231238597815812,2231624991249884,2231994346765634,
     2232347736722468,2232686144665663,2233010474325699,2233321557544631,
     2233620161275830,2233906993781039,2234182710130112,2234447917093281,
     2234703177502812,2234949014149981,2235185913274123,2235414327692697,
     2235634679614740,2235847363174420,2236052746716668,2236251174862705,
     2236442970379808,2236628435876608,2236807855342616,2236981495548416,
     2237149607321006,2237312426707072,2237470176035519,2237623064889274,
     2237771290995262,2237915041040474,2238054491421185,2238189808931596,
     2238321151397547,2238448668260322,2238572501115061,2238692784207837,
     2238809644895031,2238923204068302,2239033576548092,2239140871448347,
     2239245192514865,2239346638439450,2239445303151863,2239541276091355,
     2239634642459413,2239725483455210,2239813876495104,2239899895417414,
     2239983610673598,2240065089506859,2240144396119109,2240221591827156,
     2240296735208897,2240369882240222,2240441086423317,2240510398906937,
     2240577868599239,2240643542273660,2240707464668327,2240769678579424,
     2240830224948918,2240889142947021,2240946470049710,2241002242111632,
     2241056493434688,2241109256832545,2241160563691345,2241210444026824,
     2241258926538069,2241306038658085,2241351806601384,2241396255408737,
     2241439408989263,2241481290159988,2241521920683014,2241561321300414,
     2241599511766981,2241636510880914,2241672336512567,2241707005631317,
     2241740534330669,2241772937851645,2241804230604542,2241834426189118,
     2241863537413270,2241891576310240,2241918554154426,2241944481475803,
     2241969368073032,2241993223025259,2242016054702647,2242037870775672,
     2242058678223188,2242078483339294,2242097291739004,2242115108362739,
     2242131937479636,2242147782689690,2242162646924702,2242176532448058,
     2242189440853303,2242201373061504,2242212329317384,2242222309184204,
     2242231311537365,2242239334556685,2242246375717338,2242252431779384,
     2242257498775863,2242261571999386,2242264645987166,2242266714504423,
     2242267770526080,2242267806216682,2242266812908434,2242264781077261,
     2242261700316790,2242257559310117,2242252345799249,2242246046552055,
     2242238647326588,2242230132832599,2242220486690050,2242209691384432,
     2242197728218658,2242184577261284,2242170217290794,2242154625735654,
     2242137778609814,2242119650443302,2242100214207531,2242079441234882,
     2242057301132111,2242033761687055,2242008788768083,2241982346215658,
     2241954395725333,2241924896721420,2241893806220494,2241861078683807,
     2241826665857576,2241790516600019,2241752576693859,2241712788642894,
     2241671091451056,2241627420382213,2241581706698751,2241533877376746,
     2241483854795259,2241431556397014,2241376894317324,2241319774977796,
     2241260098640839,2241197758920517,2241132642244683,2241064627262631,
     2240993584191722,2240919374095516,2240841848084869,2240760846432212,
     2240676197587764,2240587717084761,2240495206318733,2240398451183547,
     2240297220544145,2240191264522592,2240080312570135,2239964071293311,
     2239842221996510,2239714417896679,2239580280957705,2239439398282173,
     2239291317986176,2239135544468183,2238971532964959,2238798683265249,
     2238616332424332,2238423746288075,2238220109591870,2238004514345197,
     2237775946143192,2237533267957802,2237275200846732,2237000300869931,
     2236706931309079,2236393229029127,2236057063479481,2235695986373225,
     2235307169458838,2234887326941556,2234432617919425,2233938522519742,
     2233399683022654,2232809697779175,2232160850599794,2231443750584617,
     2230646845562145,2229755753817960,2228752329126507,2227613325162477,
     2226308442121145,2224797391720369,2223025347823800,2220915633329775,
     2218357446086993,2215184158448627,2211132412537323,2205758503851011,
     2198248265654920,2186916352102052,2167562552481677,2125549880839429};

static double wi[ZIGGURAT_TABLE_SIZE] =
    {17367254121656703e-31,9558660348275145e-32,12708704832820278e-32,
     14909740960986864e-32,16658733630346416e-32,18136120809053487e-32,
     1942972015219358e-31,20589500627632916e-32,21646860576118966e-32,
     2262294039150043e-31,23532718913376864e-32,24387234556800803e-32,
     25194879828681465e-32,2596219977196592e-31,26694407473112964e-32,
     2739572968463095e-31,280696460019946e-30,28719058903642897e-32,
     29346417484275224e-32,29953809336344285e-32,30543030006769113e-32,
     3111563633851158e-31,3167298801818414e-31,3221628035016365e-31,
     32746570407564125e-32,33264798116476e-29,337718034169968e-30,
     34268340352771636e-32,34755088731390227e-32,3523266384567022e-31,
     3570162463362898e-31,3616248057128073e-31,36615697529342477e-32,
     3706170277693123e-31,37500889278448874e-32,3793361940125627e-31,
     38360228129389374e-32,3878102586096749e-31,3919630085297984e-31,
     39606321365983254e-32,40011337552278087e-32,4041158312387907e-31,
     4080727683070036e-31,4119862377455137e-31,41585816580575855e-32,
     41969036444492247e-32,4234845407127582e-31,42724230518658345e-32,
     43096517956924877e-32,4346546035489394e-31,4383119410062289e-31,
     4419384856424202e-31,4455354660935343e-31,4491040505860591e-31,
     4526453511835132e-31,45616042766683e-29,4596502910863464e-31,
     4631159070186941e-31,4665581985579899e-31,469978049067346e-30,
     4733763047137822e-31,4767537768070579e-31,4801112439606964e-31,
     4834494540915173e-31,4867691262722585e-31,4900709524503576e-31,
     4933555990446197e-31,4966237084303158e-31,499875900322208e-30,
     5031127730640677e-31,5063349048324261e-31,5095428547615612e-31,
     5127371639960692e-31,5159183566767805e-31,5190869408652579e-31,
     5222434094116442e-31,52538824077020155e-32,5285218997665102e-31,
     5316448383199491e-31,5347574961247755e-31,5378603012928409e-31,
     5409536709607314e-31,5440380118638932e-31,5471137208800966e-31,
     550181185544408e-30,5532407845376661e-31,5562928881503102e-31,
     5593378587232605e-31,5623760510674315e-31,5654078128633358e-31,
     5684334850421336e-31,5714534021493849e-31,5744678926926726e-31,
     5774772794741848e-31,5804818799092685e-31,5834820063319006e-31,
     5864779662879593e-31,589470062817121e-30,5924585947241581e-31,
     5954438568403615e-31,598426140275769e-30,601405732662843e-30,
     6043829183921996e-31,6073579788409578e-31,6103311925942512e-31,
     6133028356604082e-31,6162731816802865e-31,6192425021312213e-31,
     6222110665260248e-31,6251791426074554e-31,6281469965385542e-31,
     6311148930892342e-31,6340830958194888e-31,6370518672595733e-31,
     640021469087503e-30,6429921623041988e-31,645964207406601e-30,
     648937864559066e-30,6519133937633505e-31,6548910550274845e-31,
     6578711085338253e-31,6608538148065851e-31,6638394348791179e-31,
     6668282304612498e-31,6698204641069389e-31,6728163993825439e-31,
     6758163010359885e-31,6788204351671041e-31,681829069399439e-30,
     6848424730538249e-31,6878609173239948e-31,6908846754545526e-31,
     6939140229215998e-31,696949237616333e-30,6999906000319335e-31,
     7030383934540792e-31,7060929041554193e-31,7091544215943653e-31,
     7122232386185626e-31,7152996516734219e-31,7183839610161045e-31,
     7214764709353755e-31,7245774899777502e-31,72768733118038725e-32,
     7308063123111988e-31,7339347561166714e-31,7370729905779203e-31,
     7402213491755235e-31,7433801711637146e-31,7465498018545449e-31,
     7497305929126601e-31,7529229026613742e-31,7561270964007667e-31,
     7593435467385694e-31,7625726339346621e-31,7658147462600412e-31,
     7690702803711903e-31,7723396417008341e-31,7756232448661274e-31,
     778921514095401e-30,7822348836746627e-31,7855637984151357e-31,
     7889087141432085e-31,7922700982142658e-31,7956484300519808e-31,
     7990442017147628e-31,8024579184911813e-31,8058900995263265e-31,
     8093412784812165e-31,812812004227522e-30,8163028415800651e-31,
     8198143720697359e-31,8233471947596931e-31,8269019271079405e-31,
     8304792058796362e-31,834079688112767e-30,8377040521411316e-31,
     8413529986789175e-31,8450272519715296e-31,8487275610177406e-31,
     85245470086869e-29,8562094740097588e-31,8599927118319072e-31,
     86380527619967175e-32,8676480611237092e-31,8715219945465259e-31,
     8754280402508787e-31,8793671999012706e-31,8833405152300122e-31,
     88734907038049e-29,8913939944215902e-31,8954764640486935e-31,
     8995977064883017e-31,9037590026252085e-31,9079616903732087e-31,
     9122071683126914e-31,9164968996211253e-31,9208324163254476e-31,
     9252153239087913e-31,9296473063078686e-31,9341301313417584e-31,
     938665656617903e-30,9432558359669126e-31,9479027264644209e-31,
     95260849610588e-29,957375432209002e-30,962205950628746e-30,
     9671026058815726e-31,972068102289435e-30,9771053062699983e-31,
     9822172599183368e-31,9874071960473548e-31,9926785548800904e-31,
     9980350026176626e-31,10034804521429213e-31,10090190861630543e-31,
     10146553831460223e-31,10203941464676316e-31,1026240537260681e-30,
     10322001115479755e-31,10382788623508751e-31,10444832675993878e-31,
     10508203448348659e-31,1057297713900341e-30,10639236690670377e-31,
     10707072623626628e-31,107765840026618e-29,10847879564397177e-31,
     10921079038143372e-31,109963147017795e-29,11073733224929686e-31,
     11153497865847152e-31,11235791107104895e-31,11320817840158973e-31,
     11408809242576976e-31,1150002753783406e-30,11594771891443527e-31,
     11693385786905373e-31,1179626635295029e-30,11903876299277459e-31,
     1201675939253847e-30,12135560818661637e-31,12261054417445396e-31,
     12394179789158183e-31,12536093926597603e-31,1268824481425016e-30,
     12852479319091384e-31,13031206634685398e-31,13227655770190893e-31,
     13446300925006917e-31,13693606835124475e-31,13979436672771461e-31,
     14319989869657897e-31,14744848603594667e-31,1531787274160907e-30,
     16227698675312968e-31};

static double fi[ZIGGURAT_TABLE_SIZE] =
    {1.,.9771017012827331,.9598790918124174,.9451989534530794,
     .9320600759689914,.9199915050483614,.9087264400605639,.898095921906305,
     .8879846607634008,.8783096558161477,.869008688043794,.8600336212030095,
     .8513462584651245,.842915653118442,.8347162929929313,.8267268339520951,
     .8189291916094156,.8113078743182208,.8038494831763903,.7965423304282554,
     .7893761435711993,.7823418326598627,.775431304986139,.7686373158033355,
     .7619533468415471,.7553735065117552,.7488924472237273,.7425052963446368,
     .7362075981312672,.729995264565803,.7238645334728822,.717811932634902,
     .711834248882359,.7059285013367979,.7000919181404905,.694321916130033,
     .6886160830085275,.6829721616487918,.6773880362225135,.6718617199007669,
     .6663913439123812,.6609751477802419,.6556114705832252,.650298743114295,
     .6450354808242524,.6398202774564395,.6346517992909606,.6295287799281287,
     .6244500155502747,.6194143606090396,.6144207238920772,.6094680649288958,
     .6045553907005499,.599681752622168,.5948462437709915,.5900479963357922,
     .5852861792663006,.5805599961036837,.5758686829752109,.5712115067380753,
     .5665877632589521,.5619967758172782,.5574378936214867,.5529104904285204,
     .5484139632579217,.5439477311926505,.5395112342595453,.5351039323830201,
     .5307253044061945,.5263748471741873,.5220520746747954,.5177565172322012,
     .5134877207497434,.5092452459981365,.5050286679458292,.5008375751284826,
     .49667156905479676,.4925302636461491,.4884132847077125,.48432026942891204,
     .48025086591125016,.4762047327216842,.4721815384698837,.46818096140782267,
     .4642026890502793,.460246417814924,.45631185268077407,.4523987068638829,
     .4485067015092144,.4446355653977281,.4407850346677702,.43695485254992955,
     .43314476911457434,.42935454103134185,.42558393133990086,.4218327092313535,
     .41810064983968476,.4143875340427069,.41069314827198344,.40701728433124823,
     .4033597392228692,.399720314981932,.3960988185175474,.39249506146101104,
     .3889088600204649,.38534003484173424,.3817884108750316,.37825381724723833,
     .37473608713949164,.37123505766982157,.3677505697805964,.36428246813054976,
     .36083060099117586,.3573948201472906,.35397498080156936,.35057094148288126,
     .34718256395825153,.3438097131482915,.3404522570459456,.33711006663841303,
     .33378301583210873,.3304709813805373,.32717384281495887,.32389148237773235,
     .3206237849582305,.3173706380312227,.3141319315976305,.310907558127564,
     .307697412505554,.30450139197789644,.30131939610203423,.29815132669790145,
     .2949970878011627,.2918565856182811,.28872972848335393,.28561642681665805,
     .28251659308484933,.27943014176276515,.2763569892967811,.27329705406967564,
     .2702502563669598,.26721651834463167,.2641957639983174,.2611879191337636,
     .25819291133864797,.2552106699556771,.25224112605694377,.2492842124195167,
     .24633986350223877,.24340801542371202,.24048860594144916,.23758157443217368,
     .2346868618732527,.23180441082524855,.22893416541557743,.22607607132326474,
     .22323007576478943,.22039612748101145,.21757417672517823,.2147641752520084,
     .21196607630785277,.20917983462193548,.20640540639867916,.20364274931112133,
     .20089182249543117,.19815258654653795,.1954250035148854,.19270903690432864,
     .19000465167119293,.18731181422451676,.18463049242750437,.18196065560021638,
     .17930227452353026,.17665532144440646,.17401977008249914,.17139559563815535,
     .16878277480185,.1661812857651097,.1635911082329826,.16101222343811727,
     .1584446141565199,.15588826472506426,.15334316106083742,.15080929068240986,
     .1482866427331284,.14577520800653765,.14327497897404687,.14078594981496803,
     .138308116449064,.13584147657175705,.13338602969216254,.13094177717412792,
     .12850872228047336,.12608687022065,.12367622820205106,.12127680548523516,
     .11888861344334545,.11651166562603685,.11414597782825504,.1117915681642454,
     .10944845714720981,.10711666777507266,.10479622562286683,.10248715894230599,
     .10018949876917177,.09790327903921535,.09562853671335306,.09336531191302634,
     .09111364806670041,.08887359206859394,.08664519445086755,.08442850957065445,
     .0822235958134955,.08003051581494733,.07784933670237201,.07568013035919481,
     .07352297371424082,.07137794905914183,.06924514439725017,.06712465382802392,
     .06501657797147035,.06292102443797778,.06083810834975175,.05876795292113793,
     .056710690106399425,.05466646132507786,.05263541827697361,.05061772386112175,
     .04861355321603513,.04662309490208967,.044646552251446515,.042684144916619336,
     .04073611065607874,.0388027074046569,.03688421568869112,.034980941461833046,
     .033093219458688684,.031221417192023686,.02936593975823011,.027527235669693315,
     .02570580400863265,.023902203305873237,.022117062707379908,.020351096230109344,
     .018605121275783343,.016880083152595836,.015177088307982065,.013497450601780796,
     .0118427578579431,.0102149714397311,.008616582769422912,.007050875471392109,
     .005522403299264755,.0040379725933718715,.002609072746106362,.0012602859304985975};

// Tables for exprnd
static ZIGINT ke[ZIGGURAT_TABLE_SIZE] =
    {3985772928715748,               0,2742928985168065,3438700186803721,
     3744780257810519,3914896975372863,4022625697542798,4096776410635450,
     4150853606149210,4192001604687417,4224344877584101,4250427292531740,
     4271901371161554,4289886428824118,4305167164135199,4318309783140431,
     4329732973408940,4339752937704679,4348612900760388,4356502988721768,
     4363573953227346,4369946852445020,4375720012348349,4380974119031481,
     4385776001930298,4390181484145305,4394237557465219,4397984061535398,
     4401454994146430,4404679543790856,4407682910787985,4410486965794400,
     4413110782053579,4415571068741702,4417882526198713,4420058138987325,
     4422109419110772,4424046609003130,4425878851844253,4427614335173868,
     4429260412563040,4430823707156475,4432310200160197,4433725306767517,
     4435073941555377,4436360575016074,4437589282595121,4438763787369085,
     4439887497305303,4440963537889317,4441994780778252,4442983869033585,
     4443933239400428,4444845142028910,4445721657973833,4446564714759241,
     4447376100252993,4448157475061632,4448910383626429,4449636264176642,
     4450336457674983,4451012215872352,4451664708573597,4452295030203006,
     4452904205747010,4453493196141906,4454062903166143,4454614173889474,
     4455147804725090,4455664545125435,4456165100957688,4456650137590828,
     4457120282722585,4457576128971459,4458018236256245,4458447133983073,
     4458863323057847,4459267277740095,4459659447352586,4460040257859578,
     4460410113325310,4460769397263133,4461118473884710,4461457689257740,
     4461787372379910,4462107836175980,4462419378424319,4462722282618581,
     4463016818769709,4463303244152965,4463581804004301,4463852732169940,
     4464116251712773,4464372575478779,4464621906626490,4464864439122178,
     4465100358203284,4465329840812355,4465553056003596,4465770165323939,
     4465981323170417,4466186677125455,4466386368271563,4466580531486827,
     4466769295722448,4466952784263502,4467131114974006,4467304400527265,
     4467472748622447,4467636262188208,4467795039574164,4467949174730939,
     4468098757379442,4468243873170018,4468384603832024,4468521027314373,
     4468653217917530,4468781246417428,4468905180181701,4469025083278642,
     4469141016579234,4469253037852582,4469361201855066,4469465560413474,
     4469566162502383,4469663054316032,4469756279334881,4469845878387080,
     4469931889704995,4470014348976986,4470093289394551,4470168741694984,
     4470240734199652,4470309292847996,4470374441227332,4470436200598525,
     4470494589917605,4470549625853344,4470601322800852,4470649692891185,
     4470694745996980,4470736489734116,4470774929459349,4470810068263924,
     4470841906963074,4470870444081369,4470895675833821,4470917596102651,
     4470936196409614,4470951465883737,4470963391224346,4470971956659198,
     4470977143897542,4470978932077904,4470977297710362,4470972214613072,
     4470963653842747,4470951583618802,4470935969240827,4470916772999009,
     4470893954077117,4470867468447603,4470837268758338,4470803304210460,
     4470765520426769,4470723859310029,4470678258890503,4470628653161980,
     4470574971905457,4470517140499614,4470455079717082,4470388705505446,
     4470317928751818,4470242655029689,4470162784326669,4470078210751556,
     4469988822219058,4469894500110287,4469795118907000,4469690545797298,
     4469580640250319,4469465253557163,4469344228335006,4469217397991048,
     4469084586142556,4468945605988875,4468800259630802,4468648337332217,
     4468489616718259,4468323861903709,4468150822544456,4467970232804102,
     4467781810226787,4467585254506222,4467380246139658,4467166444954116,
     4466943488490515,4466710990229518,4466468537640691,4466215690034133,
     4465951976190801,4465676891744455,4465389896284247,4465090410142477,
     4464777810826750,4464451429049612,4464110544301482,4463754379904174,
     4463382097472202,4462992790697122,4462585478355953,4462159096427753,
     4461712489182116,4461244399078944,4460753455289386,4460238160612098,
     4459696876515553,4459127805983956,4458528973779075,4457898203649722,
     4457233091920646,4456530976767892,4455788902331217,4455003576616607,
     4454171321891082,4453288015951104,4452349022232651,4451349106194827,
     4450282334707462,4449141954247903,4447920242480611,4446608326137821,
     4445195955871677,4443671225661690,4442020220072463,4440226566619900,
     4438270861888260,4436129927556552,4433775834104270,4431174602388627,
     4428284451100006,4425053392146958,4421415870372502,4417287970124084,
     4412560416174562,4407088078325945,4400673742272494,4393042098597073,
     4383796248451589,4372341169422858,4357740343059956,4338425130125967,
     4311541827049177,4271262897902398,4203411844498905,4061213381260384};

static double we[ZIGGURAT_TABLE_SIZE] =
    {19311480126418366e-31,1417802848791084e-32,23278824993382457e-33,
     30487830247064326e-33,3666569771447489e-32,4217930218928974e-32,
     4722256155686277e-32,51911915446217885e-33,5632347108395505e-32,
     6051008260642765e-32,645101650967275e-31,6835264680370054e-32,
     7205993957468906e-32,7564981553739299e-32,7913664396195108e-32,
     8253223556351894e-32,8584643616885051e-32,8908755486564743e-32,
     9226267962966373e-32,9537791450529272e-32,9843856087455926e-32,
     10144925809006294e-32,10441409405585343e-32,10733669323436384e-32,
     1102202874567019e-31,11306777346479334e-32,11588176009705533e-32,
     11866460730417886e-32,1214184586569436e-31,12414526862326387e-32,
     12684682560606153e-32,12952477151912284e-32,1321806185153881e-31,
     13481576335745447e-32,13743149982367625e-32,14002902946807862e-32,
     14260947099321287e-32,14517386844829297e-32,14772319842763584e-32,
     15025837641447456e-32,15278026239101652e-32,15528966581595696e-32,
     1577873500545958e-31,1602740363335091e-31,16275040728083524e-32,
     16521711010420076e-32,16767475945078279e-32,17012393998770646e-32,
     17256520873568226e-32,17499909718432365e-32,17742611321380505e-32,
     17984674284430714e-32,18226145183195818e-32,18467068712763576e-32,
     18707487821298258e-32,18947443832625902e-32,19186976558915997e-32,
     19426124404443042e-32,19664924461299023e-32,19903412597830144e-32,
     20141623540485899e-32,20379590949693882e-32,2061734749030844e-31,
     2085492489712377e-31,21092354035891528e-32,21329664960238294e-32,
     21566886964838972e-32,2180404863516701e-31,22041177894111562e-32,
     2227830204572395e-31,2251544781633135e-31,22752641393233694e-32,
     22989908461180186e-32,23227274236804366e-32,23464763501180916e-32,
     2370240063065339e-31,23940209626069303e-32,2417821414054771e-31,
     24416437505894123e-32,24654902757768304e-32,2489363265970225e-31,
     2513264972605797e-31,2537197624400795e-31,2561163429461499e-31,
     2585164577308239e-31,26092032408240577e-32,2633281578133145e-31,
     2657401734414762e-31,2681565843657999e-31,2705776030362351e-31,
     27300344111887955e-32,27543430965657624e-32,2778704192254128e-31,
     2803119800875143e-31,28275920234049704e-32,2852122960639331e-31,
     28767147146315804e-32,29013693901073754e-32,29260890958589514e-32,
     29508759461219033e-32,2975732061937252e-31,3000659572501474e-31,
     3025660616507079e-31,3050737343476251e-31,3075891915089994e-31,
     31011265065151543e-32,3126443307731675e-31,31518445248623523e-32,
     31773323815073683e-32,32029091200858335e-32,32285770031865573e-32,
     3254338314930261e-31,3280195362345436e-31,3306150476760074e-31,
     3332206015211484e-31,33583643618764577e-32,33846279295240445e-32,
     34109991609932597e-32,34374805306980633e-32,34640745461620167e-32,
     3490783749585068e-31,3517610719444983e-31,3544558072136013e-31,
     3571628463647465e-31,35988245912849274e-32,3626149195437003e-31,
     36536050613905045e-32,36811950211971757e-32,3708921955595139e-31,
     37367887959883854e-32,3764798526487784e-31,37929541860172334e-32,
     3821258870488753e-31,38497157350504876e-32,3878327996411799e-31,
     39070989352498183e-32,3936031898702075e-31,3965130302950038e-31,
     3994397635898684e-31,40238374599574693e-32,40534534149283966e-32,
     4083249221007178e-31,41132286819038357e-32,4143395688089474e-31,
     417375422017632e-30,42043083524385856e-32,4235062256482152e-31,
     4266020205071558e-31,42971865761233266e-32,43285658568752094e-32,
     4360162648241568e-31,43919816693657415e-32,4424027762380992e-31,
     4456305897392361e-31,4488821177692617e-31,4521578845226347e-31,
     4554584286317242e-31,4587843037674623e-31,4621360792696427e-31,
     4655143408087069e-31,4689196910809916e-31,4723527505395548e-31,
     4758141581628553e-31,4793045722637247e-31,4828246713412587e-31,
     4863751549784512e-31,489956744788614e-30,4935701854138577e-31,
     4972162455791703e-31,5008957192059114e-31,5046094265888434e-31,
     5083582156411624e-31,5121429632123542e-31,5159645764841062e-31,
     5198239944499494e-31,5237221894847848e-31,5276601690109886e-31,
     531638977268369e-30,535659697195905e-30,5397234524338979e-31,
     5438314094559637e-31,547984779841163e-30,5521848226975234e-31,
     5564328472492872e-31,5607302156013967e-31,5650783456960506e-31,
     5694787144776348e-31,5739328612839635e-31,5784423914835991e-31,
     5830089803810586e-31,5876343774140057e-31,5923204106690931e-31,
     5970689917460091e-31,6018821210025236e-31,6067618932170007e-31,
     6117105037089722e-31,616730254963062e-30,6218235638068533e-31,
     6269929691993326e-31,6322411406934211e-31,6375708876439426e-31,
     6429851692413595e-31,6484871054618903e-31,6540799890364481e-31,
     6597672985544566e-31,6655527128343343e-31,6714401267106488e-31,
     677433668409101e-30,6835377187051274e-31,6897569320906848e-31,
     6960962602074885e-31,7025609778445959e-31,7091567118449584e-31,
     7158894733208553e-31,7227656936438121e-31,7297922647529085e-31,
     7369765844191243e-31,7443266072160415e-31,7518509020832513e-31,
     7595587175337749e-31,7674600557578427e-31,7755657571215791e-31,
     7838875968622858e-31,792438396157355e-30,8012321502113083e-31,
     8102841765913146e-31,8196112877806125e-31,8292319928581809e-31,
     8391667344146798e-31,849438168364877e-30,8600714963334941e-31,
     8710948629387904e-31,882539833807214e-30,8944419748519865e-31,
     9068415597131669e-31,9197844409811865e-31,9333231329422952e-31,
     9475181706524984e-31,9624398345658476e-31,978170365478442e-30,
     994806847238388e-30,1012465014428832e-30,10312843657756166e-31,
     1051435160404455e-30,10731281954224043e-31,10966288068517408e-31,
     1122277490935032e-30,11505212963006663e-31,11819635283304206e-31,
     12174462832361815e-31,12581958069755114e-31,13060984107128082e-31,
     13642786158057857e-31,14384889932178723e-31,15412190700064194e-31,
     17091034077168055e-31};

static double fe[ZIGGURAT_TABLE_SIZE] =
    {              1.0,.9381436808621746,.9004699299257464,.8717043323812036,
     .8477855006239896,.8269932966430503,.8084216515230084,.7915276369724956,
     .7759568520401156,.7614633888498963,.7478686219851951,.7350380924314235,
     .722867659593572,.711274760805076,.7001926550827882,.689566496117078,
     .6793505722647654,.6695063167319247,.6600008410789997,.650805833414571,
     .6418967164272661,.6332519942143661,.6248527387036659,.6166821809152077,
     .608725382079622,.6009689663652322,.5934009016917334,.586010318477268,
     .578787358602845,.5717230486648258,.5648091929124002,.5580382822625874,
     .5514034165406413,.5448982376724396,.5385168720028619,.5322538802630432,
     .5261042139836197,.5200631773682336,.5141263938147486,.5082897764106429,
     .5025495018413477,.49690198724154955,.49134386959403253,.4858719873418849,
     .4804833639304542,.4751751930373774,.46994482528396,.4647897562504262,
     .4597076156421377,.45469615747461545,.449753251162755,.4448768734145485,
     .4400651008423539,.4353161032156366,.43062813728845883,.42599954114303434,
     .4214287289976166,.4169141864330029,.4124544659971612,.4080481831520324,
     .4036940125305303,.3993906844752311,.39513698183329016,.3909317369847971,
     .38677382908413765,.38266218149600983,.3785957594095808,.37457356761590216,
     .370594648435146,.36665807978151416,.3627629733548178,.3589084729487498,
     .35509375286678746,.35131801643748334,.347580494621637,.3438804447045024,
     .34021714906678,.33658991402867755,.332998068761809,.3294409642641363,
     .3259179723935562,.3224284849560891,.31897191284495724,.31554768522712895,
     .31215524877417955,.3087940669345602,.30546361924459026,.3021634006756935,
     .2988929210155818,.2956517042812612,.2924392881618926,.28925522348967775,
     .2860990737370768,.28297041453878075,.27986883323697287,.27679392844851736,
     .27374530965280297,.27072259679906,.2677254199320448,.2647534188350622,
     .2618062426893629,.25888354974901623,.2559850070304154,.25311029001562946,
     .2502590823688623,.24743107566532763,.2446259691318921,.24184346939887721,
     .23908329026244918,.23634515245705964,.23362878343743335,.2309339171696274,
     .2282602939307167,.22560766011668407,.22297576805812017,.2203643758433595,
     .21777324714870053,.21520215107537868,.21265086199297828,.21011915938898826,
     .20760682772422204,.2051136562938377,.20263943909370902,.20018397469191127,
     .19774706610509887,.19532852067956322,.19292814997677132,.1905457696631954,
     .1881811994042543,.1858342627621971,.18350478709776746,.1811926034754963,
     .1788975465724783,.17661945459049488,.1743581691713535,.17211353531532006,
     .16988540130252766,.1676736186172502,.165478041874936,.16329852875190182,
     .16113493991759203,.1589871389693142,.15685499236936523,.15473836938446808,
     .15263714202744286,.1505511850010399,.1484803756438668,.14642459387834494,
     .14438372216063478,.1423576454324722,.14034625107486245,.1383494288635802,
     .13636707092642886,.13439907170221363,.13244532790138752,.13050573846833077,
     .12858020454522817,.12666862943751067,.12477091858083096,.12288697950954514,
     .12101672182667483,.11916005717532768,.11731689921155557,.11548716357863353,
     .11367076788274431,.1118676316700563,.11007767640518538,.1083008254510338,
     .10653700405000166,.10478613930657017,.10304816017125772,.10132299742595363,
     .09961058367063713,.0979108533114922,.0962237425504328,.09454918937605586,
     .09288713355604354,.09123751663104016,.08960028191003286,.08797537446727022,
     .08636274114075691,.08476233053236812,.08317409300963238,.08159798070923742,
     .0800339475423199,.07848194920160642,.0769419431704805,.07541388873405841,
     .07389774699236475,.07239348087570874,.07090105516237183,.06942043649872875,
     .0679515934219366,.06649449638533977,.06504911778675375,.06361543199980733,
     .062193415408540995,.06078304644547963,.059384305633420266,.05799717563120066,
     .05662164128374288,.05525768967669704,.05390531019604609,.05256449459307169,
     .05123523705512628,.04991753428270637,.0486113855733795,.04731679291318155,
     .04603376107617517,.04476229773294328,.04350241356888818,.042254122413316234,
     .04101744138041482,.039792391023374125,.03857899550307486,.03737728277295936,
     .03618728478193142,.03500903769739741,.03384258215087433,.032687963508959535,
     .03154523217289361,.030414443910466604,.029295660224637393,.028188948763978636,
     .0270943837809558,.026012046645134217,.024942026419731783,.02388442051155817,
     .02283933540638524,.02180688750428358,.020787204072578117,.019780424338009743,
     .01878670074469603,.01780620041091136,.016839106826039948,.015885621839973163,
     .014945968011691148,.014020391403181938,.013109164931254991,.012212592426255381,
     .011331013597834597,.010464810181029979,.00961441364250221,.008780314985808975,
     .00796307743801704,.007163353183634984,.006381905937319179,.005619642207205483,
     .004877655983542392,.004157295120833795,.003460264777836904,.002788798793574076,
     .0021459677437189063,.0015362997803015724,.0009672692823271745,.00045413435384149677};


/*
 * Here is the guts of the algorithm. As Marsaglia and Tsang state the
 * algorithm in their paper
 *
 * 1) Calculate a random signed integer j and let i be the index
 *     provided by the rightmost 8-bits of j
 * 2) Set x = j * w_i. If j < k_i return x
 * 3) If i = 0, then return x from the tail
 * 4) If [f(x_{i-1}) - f(x_i)] * U < f(x) - f(x_i), return x
 * 5) goto step 1
 *
 * Where f is the functional form of the distribution, which for a normal
 * distribution is exp(-0.5*x*x)
 */

/* NOTE: This is identical to randmtzig_gv_randn() below except for the random number generation */
double randmtzig_randn (dsfmt_t *dsfmt)
{
  while (1)
    {
        /* arbitrary mantissa (selected by randi, with 1 bit for sign) */
        const randmtzig_uint64_t r = randi(dsfmt);
        const randmtzig_int64_t rabs=r>>1;
        const int idx = (int)(rabs&0xFF);
        const double x = ( r&1 ? -rabs : rabs) * wi[idx];

        if (rabs < (randmtzig_int64_t)ki[idx]) {
            return x;        /* 99.3% of the time we return here 1st try */
        } else if (idx == 0) {
            /* As stated in Marsaglia and Tsang
             *
             * For the normal tail, the method of Marsaglia[5] provides:
             * generate x = -ln(U_1)/r, y = -ln(U_2), until y+y > x*x,
             * then return r+x. Except that r+x is always in the positive
             * tail!!!! Any thing random might be used to determine the
             * sign, but as we already have r we might as well use it
             *
             * [PAK] but not the bottom 8 bits, since they are all 0 here!
             */
            double xx, yy;
            do {
                xx = - ZIGGURAT_NOR_INV_R * log (randu(dsfmt));
                yy = - log (randu(dsfmt));
            }
            while ( yy+yy <= xx*xx);
            return (rabs&0x100 ? -ZIGGURAT_NOR_R-xx : ZIGGURAT_NOR_R+xx);
        } else if ((fi[idx-1] - fi[idx]) * randu(dsfmt) + fi[idx] < exp(-0.5*x*x)) {
            return x;
        }

    }
}

/* NOTE: This is identical to randmtzig_randn() above except for the random number generation */
double randmtzig_gv_randn (void)
{
  while (1)
    {
        /* arbitrary mantissa (selected by NRANDI, with 1 bit for sign) */
        const randmtzig_uint64_t r = NRANDI;
        const randmtzig_int64_t rabs=r>>1;
        const int idx = (int)(rabs&0xFF);
        const double x = ( r&1 ? -rabs : rabs) * wi[idx];

        if (rabs < (randmtzig_int64_t)ki[idx]) {
            return x;        /* 99.3% of the time we return here 1st try */
        } else if (idx == 0) {
            /* As stated in Marsaglia and Tsang
             *
             * For the normal tail, the method of Marsaglia[5] provides:
             * generate x = -ln(U_1)/r, y = -ln(U_2), until y+y > x*x,
             * then return r+x. Except that r+x is always in the positive
             * tail!!!! Any thing random might be used to determine the
             * sign, but as we already have r we might as well use it
             *
             * [PAK] but not the bottom 8 bits, since they are all 0 here!
             */
            double xx, yy;
            do {
                xx = - ZIGGURAT_NOR_INV_R * log (RANDU);
                yy = - log (RANDU);
            }
            while ( yy+yy <= xx*xx);
            return (rabs&0x100 ? -ZIGGURAT_NOR_R-xx : ZIGGURAT_NOR_R+xx);
        } else if ((fi[idx-1] - fi[idx]) * RANDU + fi[idx] < exp(-0.5*x*x)) {
            return x;
        }

    }
}

double randmtzig_gv_exprnd (void)
{
    while (1)
    {
        ZIGINT ri = ERANDI;
        const int idx = (int)(ri & 0xFF);
        const double x = ri * we[idx];
        if (ri < ke[idx])
            return x;        // 98.9% of the time we return here 1st try
        else if (idx == 0)
        {
            /* As stated in Marsaglia and Tsang
             *
             * For the exponential tail, the method of Marsaglia[5] provides:
             * x = r - ln(U);
             */
            return ZIGGURAT_EXP_R - log(RANDU);
    }
    else if ((fe[idx-1] - fe[idx]) * RANDU + fe[idx] < exp(-x))
        return x;
    }
}

#ifdef STANDALONE

int main(int ac, char *av[]) {
    if (ac == 1) {
        printf("Usage: randmtzig <n>\n");
        return (-1);
    }

    int n = atoi(av[1]);
    time_t t1;

    dsfmt_gv_init_gen_rand(0);

    double *p; posix_memalign((void **)&p, 16, n*sizeof(double));
    uint32_t *u; posix_memalign((void **)&u, 16, 2*n*sizeof(uint32_t));

    t1 = clock();
    dsfmt_gv_fill_array_close_open(p, n);
    printf("Uniform fill (n): %f\n", (clock() - t1) / (double) CLOCKS_PER_SEC);

    t1 = clock();
    for (int i = 0; i < n; i++)  p[i] = dsfmt_gv_genrand_close_open();
    printf("Uniform (n): %f\n", (clock() - t1) / (double) CLOCKS_PER_SEC);

    t1 = clock();
    for (int i = 0; i < 2*n; i++)  u[i] = dsfmt_gv_genrand_uint32();
    printf("Uniform 32-bit ints (2*n): %f\n", (clock() - t1) / (double) CLOCKS_PER_SEC);

    memset((void *)p, 0, n*sizeof(double));
    t1 = clock();
    for (int i = 0; i < n; i++)  p[i] = randmtzig_gv_randn();
    printf("Normal (n): %f\n", (clock() - t1) / (double) CLOCKS_PER_SEC);
    for (int i = 0; i < 10; i++)  printf("%lf\n", p[i]);

    return 0;
}

#endif