code/src/Util/rng.c

   1 // ISAAC PRNG
   2 // Initial code by Bob Jenkins, March 1996 - modified in 2008
   3 // Reference: http://burtleburtle.net/bob/rand/isaacafa.html
   4 // Further slightly modified by Benjamin Auder, 2013
   5
   6 #include "Util/rng.h"
   7 #include <stdlib.h>
   8 #include <time.h>
   9
  10 // Internal state
  11 static uint32_t mm[256];
  12 static uint32_t aa, bb, cc;
  13
  14 // initialize the (pseudo-)random generator
  15 void init_rng(int flag)
  16 {
  17     // 'bootstrap' isaac with basic bad PRNG
  18     srand(time(NULL));
  19
  20     aa = rand() & 0xffffffff;
  21     bb = rand() & 0xffffffff;
  22     cc = rand() & 0xffffffff;
  23
  24     uint32_t a, b, c, d, e, f, g, h, i;
  25     // The golden ratio:
  26     a = b = c = d = e = f = g = h = 0x9e3779b9;
  27     // Scramble it
  28     for (i = 0; i < 4; ++i)
  29         mix(a, b, c, d, e, f, g, h);
  30
  31     // Fill in mm[] with messy stuff (the seed)
  32     uint32_t seedArray[256];
  33     for (i = 0; i < 256; ++i) mm[i] = seedArray[i] = rand() & 0xffffffff;
  34     for (i = 0; i < 256; i += 8)
  35     {
  36         if (flag)
  37         {
  38             // Use all the information in the seed
  39             a += seedArray[i  ];
  40             b += seedArray[i + 1];
  41             c += seedArray[i + 2];
  42             d += seedArray[i + 3];
  43             e += seedArray[i + 4];
  44             f += seedArray[i + 5];
  45             g += seedArray[i + 6];
  46             h += seedArray[i + 7];
  47         }
  48         mix(a, b, c, d, e, f, g, h);
  49         mm[i  ] = a;
  50         mm[i + 1] = b;
  51         mm[i + 2] = c;
  52         mm[i + 3] = d;
  53         mm[i + 4] = e;
  54         mm[i + 5] = f;
  55         mm[i + 6] = g;
  56         mm[i + 7] = h;
  57     }
  58
  59     if (flag)
  60     {
  61         // Do a second pass to make all of the seed affect all of mm
  62         for (i = 0; i < 256; i += 8)
  63         {
  64             a += mm[i  ];
  65             b += mm[i + 1];
  66             c += mm[i + 2];
  67             d += mm[i + 3];
  68             e += mm[i + 4];
  69             f += mm[i + 5];
  70             g += mm[i + 6];
  71             h += mm[i + 7];
  72             mix(a, b, c, d, e, f, g, h);
  73             mm[i  ] = a;
  74             mm[i + 1] = b;
  75             mm[i + 2] = c;
  76             mm[i + 3] = d;
  77             mm[i + 4] = e;
  78             mm[i + 5] = f;
  79             mm[i + 6] = g;
  80             mm[i + 7] = h;
  81         }
  82     }
  83 }
  84
  85 // return a (pseudo-)random integer
  86 uint32_t get_rand_int()
  87 {
  88     // TODO: register variables ? (x,y,i)
  89     uint32_t x, y;
  90     static uint32_t i = 0;
  91
  92     if (i == 0)
  93     {
  94         cc = cc + 1;  // cc just gets incremented once per 256 results
  95         bb = bb + cc; // then combined with bb
  96     }
  97
  98     x = mm[i];
  99     switch (i % 4)
 100     {
 101     case 0:
 102         aa = aa^(aa << 13);
 103         break;
 104     case 1:
 105         aa = aa^(aa >> 6);
 106         break;
 107     case 2:
 108         aa = aa^(aa << 2);
 109         break;
 110     case 3:
 111         aa = aa^(aa >> 16);
 112         break;
 113     }
 114
 115     // NOTE: bits 2..9 are chosen from x but 10..17 are chosen
 116     //       from y.  The only important thing here is that 2..9 and 10..17
 117     //       don't overlap.  2..9 and 10..17 were then chosen for speed in
 118     //       the optimized version (rand.c) */
 119     // See http://burtleburtle.net/bob/rand/isaac.html
 120     // for further explanations and analysis.
 121
 122     aa        = mm[(i + 128) % 256] + aa;
 123     mm[i] = y = mm[(x >> 2) % 256]  + aa + bb;
 124     bb        = mm[(y >> 10) % 256] + x;
 125
 126     i = (i + 1) % 256;
 127
 128     return bb;
 129 }
 130
 131 // return a (pseudo-)random real number in [0,1]
 132 Real get_rand_real()
 133 {
 134     return (Real) (INV_RANDMAX * (double) get_rand_int());
 135 }