binaries/src/clustalo/src/squid/sre_random.c

   1 /* sre_random.c
   2  *
   3  * Portable random number generator, and sampling routines.
   4  *
   5  * SRE, Tue Oct  1 15:24:11 2002 [St. Louis]
   6  * CVS $Id: sre_random.c,v 1.1 2002/10/09 14:26:09 eddy Exp)
   7  */
   8
   9 #include <stdio.h>
  10 #include <stdlib.h>
  11 #include <math.h>
  12 #include "sre_random.h"
  13
  14 static int sre_randseed = 42;   /* default seed for sre_random()   */
  15
  16 /* Function: sre_random()
  17  *
  18  * Purpose:  Return a uniform deviate x, 0.0 <= x < 1.0.
  19  *
  20  *           sre_randseed is a static variable, set
  21  *           by sre_srandom(). When it is non-zero,
  22  *           we re-seed.
  23  *
  24  *           Implements L'Ecuyer's algorithm for combining output
  25  *           of two linear congruential generators, plus a Bays-Durham
  26  *           shuffle. This is essentially ran2() from Numerical Recipes,
  27  *           sans their nonhelpful Rand/McNally-esque code obfuscation.
  28  *
  29  *           Overflow errors are avoided by Schrage's algorithm:
  30  *               az % m = a(z%q) - r(z/q) (+m if <0)
  31  *           where q=m/a, r=m%a
  32  *
  33  *           Requires that long int's have at least 32 bits.
  34  *           This function uses statics and is NOT THREADSAFE.
  35  *
  36  * Reference: Press et al. Numerical Recipes in C, 1992.
  37  *
  38  * Reliable and portable, but slow. Benchmarks on wrasse,
  39  * using Linux gcc and Linux glibc rand() (see randspeed, in Testsuite):
  40  *     sre_random():    0.5 usec/call
  41  *     rand():          0.2 usec/call
  42  */
  43 double
  44 sre_random(void)
  45 {
  46   static long  rnd1;            /* random number from LCG1 */
  47   static long  rnd2;            /* random number from LCG2 */
  48   static long  rnd;             /* random number we return */
  49   static long  tbl[64];         /* table for Bays/Durham shuffle */
  50   long x,y;
  51   int i;
  52
  53   /* Magic numbers a1,m1, a2,m2 from L'Ecuyer, for 2 LCGs.
  54    * q,r derive from them (q=m/a, r=m%a) and are needed for Schrage's algorithm.
  55    */
  56   long a1 = 40014;
  57   long m1 = 2147483563;
  58   long q1 = 53668;
  59   long r1 = 12211;
  60
  61   long a2 = 40692;
  62   long m2 = 2147483399;
  63   long q2 = 52774;
  64   long r2 = 3791;
  65
  66   if (sre_randseed > 0)
  67     {
  68       rnd1 = sre_randseed;
  69       rnd2 = sre_randseed;
  70                                 /* Fill the table for Bays/Durham */
  71       for (i = 0; i < 64; i++) {
  72         x    = a1*(rnd1%q1);   /* LCG1 in action... */
  73         y    = r1*(rnd1/q1);
  74         rnd1 = x-y;
  75         if (rnd1 < 0) rnd1 += m1;
  76
  77         x    = a2*(rnd2%q2);   /* LCG2 in action... */
  78         y    = r2*(rnd2/q2);
  79         rnd2 = x-y;
  80         if (rnd2 < 0) rnd2 += m2;
  81
  82         tbl[i] = rnd1-rnd2;
  83         if (tbl[i] < 0) tbl[i] += m1;
  84       }
  85       sre_randseed = 0;         /* drop the flag. */
  86     }/* end of initialization*/
  87
  88
  89   x    = a1*(rnd1%q1);   /* LCG1 in action... */
  90   y    = r1*(rnd1/q1);
  91   rnd1 = x-y;
  92   if (rnd1 < 0) rnd1 += m1;
  93
  94   x    = a2*(rnd2%q2);   /* LCG2 in action... */
  95   y    = r2*(rnd2/q2);
  96   rnd2 = x-y;
  97   if (rnd2 < 0) rnd2 += m2;
  98
  99                         /* Choose our random number from the table... */
 100   i   = (int) (((double) rnd / (double) m1) * 64.);
 101   rnd = tbl[i];
 102                         /* and replace with a new number by L'Ecuyer. */
 103   tbl[i] = rnd1-rnd2;
 104   if (tbl[i] < 0) tbl[i] += m1;
 105
 106   return ((double) rnd / (double) m1);
 107 }
 108
 109 /* Function: sre_srandom()
 110  *
 111  * Purpose:  Initialize with a random seed. Seed must be
 112  *           >= 0 to work; we silently enforce this.
 113  */
 114 void
 115 sre_srandom(int seed)
 116 {
 117   if (seed < 0)  seed = -1 * seed;
 118   if (seed == 0) seed = 42;
 119   sre_randseed = seed;
 120 }
 121
 122 /* Function: sre_random_positive()
 123  * Date:     SRE, Wed Apr 17 13:34:32 2002 [St. Louis]
 124  *
 125  * Purpose:  Assure 0 < x < 1 (positive uniform deviate)
 126  */
 127 double
 128 sre_random_positive(void)
 129 {
 130   double x;
 131   do { x = sre_random(); } while (x == 0.0);
 132   return x;
 133 }
 134
 135 /* Function: ExponentialRandom()
 136  * Date:     SRE, Mon Sep  6 21:24:29 1999 [St. Louis]
 137  *
 138  * Purpose:  Pick an exponentially distributed random variable
 139  *           0 > x >= infinity
 140  *
 141  * Args:     (void)
 142  *
 143  * Returns:  x
 144  */
 145 double
 146 ExponentialRandom(void)
 147 {
 148   double x;
 149
 150   do x = sre_random(); while (x == 0.0);
 151   return -log(x);
 152 }
 153
 154 /* Function: Gaussrandom()
 155  *
 156  * Pick a Gaussian-distributed random variable
 157  * with some mean and standard deviation, and
 158  * return it.
 159  *
 160  * Based on RANLIB.c public domain implementation.
 161  * Thanks to the authors, Barry W. Brown and James Lovato,
 162  * University of Texas, M.D. Anderson Cancer Center, Houston TX.
 163  * Their implementation is from Ahrens and Dieter, "Extensions
 164  * of Forsythe's method for random sampling from the normal
 165  * distribution", Math. Comput. 27:927-937 (1973).
 166  *
 167  * Impenetrability of the code is to be blamed on its FORTRAN/f2c lineage.
 168  *
 169  */
 170 double
 171 Gaussrandom(double mean, double stddev)
 172 {
 173   static double a[32] = {
 174     0.0,3.917609E-2,7.841241E-2,0.11777,0.1573107,0.1970991,0.2372021,0.2776904,    0.3186394,0.36013,0.4022501,0.4450965,0.4887764,0.5334097,0.5791322,
 175     0.626099,0.6744898,0.7245144,0.7764218,0.8305109,0.8871466,0.9467818,
 176     1.00999,1.077516,1.150349,1.229859,1.318011,1.417797,1.534121,1.67594,
 177     1.862732,2.153875
 178   };
 179   static double d[31] = {
 180     0.0,0.0,0.0,0.0,0.0,0.2636843,0.2425085,0.2255674,0.2116342,0.1999243,
 181     0.1899108,0.1812252,0.1736014,0.1668419,0.1607967,0.1553497,0.1504094,
 182     0.1459026,0.14177,0.1379632,0.1344418,0.1311722,0.128126,0.1252791,
 183     0.1226109,0.1201036,0.1177417,0.1155119,0.1134023,0.1114027,0.1095039
 184   };
 185   static double t[31] = {
 186     7.673828E-4,2.30687E-3,3.860618E-3,5.438454E-3,7.0507E-3,8.708396E-3,
 187     1.042357E-2,1.220953E-2,1.408125E-2,1.605579E-2,1.81529E-2,2.039573E-2,
 188     2.281177E-2,2.543407E-2,2.830296E-2,3.146822E-2,3.499233E-2,3.895483E-2,
 189     4.345878E-2,4.864035E-2,5.468334E-2,6.184222E-2,7.047983E-2,8.113195E-2,
 190     9.462444E-2,0.1123001,0.136498,0.1716886,0.2276241,0.330498,0.5847031
 191   };
 192   static double h[31] = {
 193     3.920617E-2,3.932705E-2,3.951E-2,3.975703E-2,4.007093E-2,4.045533E-2,
 194     4.091481E-2,4.145507E-2,4.208311E-2,4.280748E-2,4.363863E-2,4.458932E-2,
 195     4.567523E-2,4.691571E-2,4.833487E-2,4.996298E-2,5.183859E-2,5.401138E-2,
 196     5.654656E-2,5.95313E-2,6.308489E-2,6.737503E-2,7.264544E-2,7.926471E-2,
 197     8.781922E-2,9.930398E-2,0.11556,0.1404344,0.1836142,0.2790016,0.7010474
 198   };
 199   static long i;
 200   static double snorm,u,s,ustar,aa,w,y,tt;
 201
 202   u = sre_random();
 203   s = 0.0;
 204   if(u > 0.5) s = 1.0;
 205   u += (u-s);
 206   u = 32.0*u;
 207   i = (long) (u);
 208   if(i == 32) i = 31;
 209   if(i == 0) goto S100;
 210   /*
 211    * START CENTER
 212    */
 213   ustar = u-(double)i;
 214   aa = *(a+i-1);
 215 S40:
 216   if(ustar <= *(t+i-1)) goto S60;
 217   w = (ustar-*(t+i-1))**(h+i-1);
 218 S50:
 219   /*
 220    * EXIT   (BOTH CASES)
 221    */
 222   y = aa+w;
 223   snorm = y;
 224   if(s == 1.0) snorm = -y;
 225   return (stddev*snorm + mean);
 226 S60:
 227   /*
 228    * CENTER CONTINUED
 229    */
 230   u = sre_random();
 231   w = u*(*(a+i)-aa);
 232   tt = (0.5*w+aa)*w;
 233   goto S80;
 234 S70:
 235   tt = u;
 236   ustar = sre_random();
 237 S80:
 238   if(ustar > tt) goto S50;
 239   u = sre_random();
 240   if(ustar >= u) goto S70;
 241   ustar = sre_random();
 242   goto S40;
 243 S100:
 244   /*
 245    * START TAIL
 246    */
 247   i = 6;
 248   aa = *(a+31);
 249   goto S120;
 250 S110:
 251   aa += *(d+i-1);
 252   i += 1;
 253 S120:
 254   u += u;
 255   if(u < 1.0) goto S110;
 256   u -= 1.0;
 257 S140:
 258   w = u**(d+i-1);
 259   tt = (0.5*w+aa)*w;
 260   goto S160;
 261 S150:
 262   tt = u;
 263 S160:
 264   ustar = sre_random();
 265   if(ustar > tt) goto S50;
 266   u = sre_random();
 267   if(ustar >= u) goto S150;
 268   u = sre_random();
 269   goto S140;
 270 }
 271
 272
 273 /* Functions: DChoose(), FChoose()
 274  *
 275  * Purpose:   Make a random choice from a normalized distribution.
 276  *            DChoose() is for double-precision vectors;
 277  *            FChoose() is for single-precision float vectors.
 278  *            Returns the number of the choice.
 279  */
 280 int
 281 DChoose(double *p, int N)
 282 {
 283   double roll;                  /* random fraction */
 284   double sum;                   /* integrated prob */
 285   int    i;                     /* counter over the probs */
 286
 287   roll    = sre_random();
 288   sum     = 0.0;
 289   for (i = 0; i < N; i++)
 290     {
 291       sum += p[i];
 292       if (roll < sum) return i;
 293     }
 294   return (int) (sre_random() * N);         /* bulletproof */
 295 }
 296 int
 297 FChoose(float *p, int N)
 298 {
 299   float roll;                   /* random fraction */
 300   float sum;                    /* integrated prob */
 301   int   i;                      /* counter over the probs */
 302
 303   roll    = sre_random();
 304   sum     = 0.0;
 305   for (i = 0; i < N; i++)
 306     {
 307       sum += p[i];
 308       if (roll < sum) return i;
 309     }
 310   return (int) (sre_random() * N);           /* bulletproof */
 311 }
 312
 313