htswanalysis/MACS/lib/gsl/gsl-1.11/randist/gausszig.c

   1 /* gauss.c - gaussian random numbers, using the Ziggurat method
   2  *
   3  * Copyright (C) 2005  Jochen Voss.
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or
   8  * (at your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful,
  11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13  * GNU General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  18  */
  19
  20 /*
  21  * This routine is based on the following article, with a couple of
  22  * modifications which simplify the implementation.
  23  *
  24  *     George Marsaglia, Wai Wan Tsang
  25  *     The Ziggurat Method for Generating Random Variables
  26  *     Journal of Statistical Software, vol. 5 (2000), no. 8
  27  *     http://www.jstatsoft.org/v05/i08/
  28  *
  29  * The modifications are:
  30  *
  31  * 1) use 128 steps instead of 256 to decrease the amount of static
  32  * data necessary.
  33  *
  34  * 2) use an acceptance sampling from an exponential wedge
  35  * exp(-R*(x-R/2)) for the tail of the base strip to simplify the
  36  * implementation.  The area of exponential wedge is used in
  37  * calculating 'v' and the coefficients in ziggurat table, so the
  38  * coefficients differ slightly from those in the Marsaglia and Tsang
  39  * paper.
  40  *
  41  * See also Leong et al, "A Comment on the Implementation of the
  42  * Ziggurat Method", Journal of Statistical Software, vol 5 (2005), no 7.
  43  *
  44  */
  45
  46
  47 #include <config.h>
  48 #include <math.h>
  49 #include <gsl/gsl_math.h>
  50 #include <gsl/gsl_rng.h>
  51 #include <gsl/gsl_randist.h>
  52
  53 /* position of right-most step */
  54 #define PARAM_R 3.44428647676
  55
  56 /* tabulated values for the heigt of the Ziggurat levels */
  57 static const double ytab[128] = {
  58   1, 0.963598623011, 0.936280813353, 0.913041104253,
  59   0.892278506696, 0.873239356919, 0.855496407634, 0.838778928349,
  60   0.822902083699, 0.807732738234, 0.793171045519, 0.779139726505,
  61   0.765577436082, 0.752434456248, 0.739669787677, 0.727249120285,
  62   0.715143377413, 0.703327646455, 0.691780377035, 0.68048276891,
  63   0.669418297233, 0.65857233912, 0.647931876189, 0.637485254896,
  64   0.62722199145, 0.617132611532, 0.607208517467, 0.597441877296,
  65   0.587825531465, 0.578352913803, 0.569017984198, 0.559815170911,
  66   0.550739320877, 0.541785656682, 0.532949739145, 0.524227434628,
  67   0.515614886373, 0.507108489253, 0.498704867478, 0.490400854812,
  68   0.482193476986, 0.47407993601, 0.466057596125, 0.458123971214,
  69   0.450276713467, 0.442513603171, 0.434832539473, 0.427231532022,
  70   0.419708693379, 0.41226223212, 0.404890446548, 0.397591718955,
  71   0.390364510382, 0.383207355816, 0.376118859788, 0.369097692334,
  72   0.362142585282, 0.355252328834, 0.348425768415, 0.341661801776,
  73   0.334959376311, 0.328317486588, 0.321735172063, 0.31521151497,
  74   0.308745638367, 0.302336704338, 0.29598391232, 0.289686497571,
  75   0.283443729739, 0.27725491156, 0.271119377649, 0.265036493387,
  76   0.259005653912, 0.253026283183, 0.247097833139, 0.241219782932,
  77   0.235391638239, 0.229612930649, 0.223883217122, 0.218202079518,
  78   0.212569124201, 0.206983981709, 0.201446306496, 0.195955776745,
  79   0.190512094256, 0.185114984406, 0.179764196185, 0.174459502324,
  80   0.169200699492, 0.1639876086, 0.158820075195, 0.153697969964,
  81   0.148621189348, 0.143589656295, 0.138603321143, 0.133662162669,
  82   0.128766189309, 0.123915440582, 0.119109988745, 0.114349940703,
  83   0.10963544023, 0.104966670533, 0.100343857232, 0.0957672718266,
  84   0.0912372357329, 0.0867541250127, 0.082318375932, 0.0779304915295,
  85   0.0735910494266, 0.0693007111742, 0.065060233529, 0.0608704821745,
  86   0.056732448584, 0.05264727098, 0.0486162607163, 0.0446409359769,
  87   0.0407230655415, 0.0368647267386, 0.0330683839378, 0.0293369977411,
  88   0.0256741818288, 0.0220844372634, 0.0185735200577, 0.0151490552854,
  89   0.0118216532614, 0.00860719483079, 0.00553245272614, 0.00265435214565
  90 };
  91
  92 /* tabulated values for 2^24 times x[i]/x[i+1],
  93  * used to accept for U*x[i+1]<=x[i] without any floating point operations */
  94 static const unsigned long ktab[128] = {
  95   0, 12590644, 14272653, 14988939,
  96   15384584, 15635009, 15807561, 15933577,
  97   16029594, 16105155, 16166147, 16216399,
  98   16258508, 16294295, 16325078, 16351831,
  99   16375291, 16396026, 16414479, 16431002,
 100   16445880, 16459343, 16471578, 16482744,
 101   16492970, 16502368, 16511031, 16519039,
 102   16526459, 16533352, 16539769, 16545755,
 103   16551348, 16556584, 16561493, 16566101,
 104   16570433, 16574511, 16578353, 16581977,
 105   16585398, 16588629, 16591685, 16594575,
 106   16597311, 16599901, 16602354, 16604679,
 107   16606881, 16608968, 16610945, 16612818,
 108   16614592, 16616272, 16617861, 16619363,
 109   16620782, 16622121, 16623383, 16624570,
 110   16625685, 16626730, 16627708, 16628619,
 111   16629465, 16630248, 16630969, 16631628,
 112   16632228, 16632768, 16633248, 16633671,
 113   16634034, 16634340, 16634586, 16634774,
 114   16634903, 16634972, 16634980, 16634926,
 115   16634810, 16634628, 16634381, 16634066,
 116   16633680, 16633222, 16632688, 16632075,
 117   16631380, 16630598, 16629726, 16628757,
 118   16627686, 16626507, 16625212, 16623794,
 119   16622243, 16620548, 16618698, 16616679,
 120   16614476, 16612071, 16609444, 16606571,
 121   16603425, 16599973, 16596178, 16591995,
 122   16587369, 16582237, 16576520, 16570120,
 123   16562917, 16554758, 16545450, 16534739,
 124   16522287, 16507638, 16490152, 16468907,
 125   16442518, 16408804, 16364095, 16301683,
 126   16207738, 16047994, 15704248, 15472926
 127 };
 128
 129 /* tabulated values of 2^{-24}*x[i] */
 130 static const double wtab[128] = {
 131   1.62318314817e-08, 2.16291505214e-08, 2.54246305087e-08, 2.84579525938e-08,
 132   3.10340022482e-08, 3.33011726243e-08, 3.53439060345e-08, 3.72152672658e-08,
 133   3.8950989572e-08, 4.05763964764e-08, 4.21101548915e-08, 4.35664624904e-08,
 134   4.49563968336e-08, 4.62887864029e-08, 4.75707945735e-08, 4.88083237257e-08,
 135   5.00063025384e-08, 5.11688950428e-08, 5.22996558616e-08, 5.34016475624e-08,
 136   5.44775307871e-08, 5.55296344581e-08, 5.65600111659e-08, 5.75704813695e-08,
 137   5.85626690412e-08, 5.95380306862e-08, 6.04978791776e-08, 6.14434034901e-08,
 138   6.23756851626e-08, 6.32957121259e-08, 6.42043903937e-08, 6.51025540077e-08,
 139   6.59909735447e-08, 6.68703634341e-08, 6.77413882848e-08, 6.8604668381e-08,
 140   6.94607844804e-08, 7.03102820203e-08, 7.11536748229e-08, 7.1991448372e-08,
 141   7.2824062723e-08, 7.36519550992e-08, 7.44755422158e-08, 7.52952223703e-08,
 142   7.61113773308e-08, 7.69243740467e-08, 7.77345662086e-08, 7.85422956743e-08,
 143   7.93478937793e-08, 8.01516825471e-08, 8.09539758128e-08, 8.17550802699e-08,
 144   8.25552964535e-08, 8.33549196661e-08, 8.41542408569e-08, 8.49535474601e-08,
 145   8.57531242006e-08, 8.65532538723e-08, 8.73542180955e-08, 8.8156298059e-08,
 146   8.89597752521e-08, 8.97649321908e-08, 9.05720531451e-08, 9.138142487e-08,
 147   9.21933373471e-08, 9.30080845407e-08, 9.38259651738e-08, 9.46472835298e-08,
 148   9.54723502847e-08, 9.63014833769e-08, 9.71350089201e-08, 9.79732621669e-08,
 149   9.88165885297e-08, 9.96653446693e-08, 1.00519899658e-07, 1.0138063623e-07,
 150   1.02247952126e-07, 1.03122261554e-07, 1.04003996769e-07, 1.04893609795e-07,
 151   1.05791574313e-07, 1.06698387725e-07, 1.07614573423e-07, 1.08540683296e-07,
 152   1.09477300508e-07, 1.1042504257e-07, 1.11384564771e-07, 1.12356564007e-07,
 153   1.13341783071e-07, 1.14341015475e-07, 1.15355110887e-07, 1.16384981291e-07,
 154   1.17431607977e-07, 1.18496049514e-07, 1.19579450872e-07, 1.20683053909e-07,
 155   1.21808209468e-07, 1.2295639141e-07, 1.24129212952e-07, 1.25328445797e-07,
 156   1.26556042658e-07, 1.27814163916e-07, 1.29105209375e-07, 1.30431856341e-07,
 157   1.31797105598e-07, 1.3320433736e-07, 1.34657379914e-07, 1.36160594606e-07,
 158   1.37718982103e-07, 1.39338316679e-07, 1.41025317971e-07, 1.42787873535e-07,
 159   1.44635331499e-07, 1.4657889173e-07, 1.48632138436e-07, 1.50811780719e-07,
 160   1.53138707402e-07, 1.55639532047e-07, 1.58348931426e-07, 1.61313325908e-07,
 161   1.64596952856e-07, 1.68292495203e-07, 1.72541128694e-07, 1.77574279496e-07,
 162   1.83813550477e-07, 1.92166040885e-07, 2.05295471952e-07, 2.22600839893e-07
 163 };
 164
 165
 166 double
 167 gsl_ran_gaussian_ziggurat (const gsl_rng * r, const double sigma)
 168 {
 169   unsigned long int i, j;
 170   int sign;
 171   double x, y;
 172
 173   const unsigned long int range = r->type->max - r->type->min;
 174   const unsigned long int offset = r->type->min;
 175
 176   while (1)
 177     {
 178       if (range >= 0xFFFFFFFF)
 179         {
 180           unsigned long int k = gsl_rng_get(r) - offset;
 181           i = (k & 0xFF);
 182           j = (k >> 8) & 0xFFFFFF;
 183         }
 184       else if (range >= 0x00FFFFFF)
 185         {
 186           unsigned long int k1 = gsl_rng_get(r) - offset;
 187           unsigned long int k2 = gsl_rng_get(r) - offset;
 188           i = (k1 & 0xFF);
 189           j = (k2 & 0x00FFFFFF);
 190         }
 191       else
 192         {
 193           i = gsl_rng_uniform_int (r, 256); /*  choose the step */
 194           j = gsl_rng_uniform_int (r, 16777216);  /* sample from 2^24 */
 195         }
 196
 197       sign = (i & 0x80) ? +1 : -1;
 198       i &= 0x7f;
 199
 200       x = j * wtab[i];
 201
 202       if (j < ktab[i])
 203         break;
 204
 205       if (i < 127)
 206         {
 207           double y0, y1, U1;
 208           y0 = ytab[i];
 209           y1 = ytab[i + 1];
 210           U1 = gsl_rng_uniform (r);
 211           y = y1 + (y0 - y1) * U1;
 212         }
 213       else
 214         {
 215           double U1, U2;
 216           U1 = 1.0 - gsl_rng_uniform (r);
 217           U2 = gsl_rng_uniform (r);
 218           x = PARAM_R - log (U1) / PARAM_R;
 219           y = exp (-PARAM_R * (x - 0.5 * PARAM_R)) * U2;
 220         }
 221
 222       if (y < exp (-0.5 * x * x))
 223         break;
 224     }
 225
 226   return sign * sigma * x;
 227 }