prng/random.c

   1 #include <stdlib.h>
   2
   3 /*
   4 32bit unsigned is asssumed
   5
   6 original bsd lagged fibonacci generator params:
   7         32byte: 7,3
   8         64byte: 15,1
   9         128byte: 31,3 (default)
  10         256byte: 63,1
  11 only the default 128byte version is implemented here
  12 in case of a smaller state buffer an lcg is used
  13 lcg params and initializer are different from the bsd code
  14
  15 implementing other variants is easy with the current code
  16 only the n,i,j should be set up apropriately in initstate and srandom
  17
  18 but a quick google codesearch shows that noone uses this api so why bother
  19 http://www.google.com/codesearch#search/&q=initstate\+*\%28+case:yes
  20
  21 there are better large period generators eg. well512
  22 this one fails on some dieharder tests
  23 (by comparison the top 32bit of a 64bit lcg passes all tests
  24 so it might make sense to just forget about the larger period
  25 and use a 64bit lcg here)
  26 */
  27
  28 static unsigned init[] = {
  29 0x00000000,0x00000001,0x3c88596c,0x5e8885db,
  30 0x8116017e,0xb4733ac5,0x0cf06d60,0x5e98c13f,
  31 0xc656dd92,0x8e625fc9,0x0438e694,0xa3a5a0e3,
  32 0x401d90e6,0x6c20f30d,0x97377908,0xd64148c7,
  33 0x3c2def7a,0xfb18b891,0xdc6318bc,0x53ae1ceb,
  34 0xaebf0d4e,0x880db455,0x3747f9b0,0x1ac2c14f,
  35 0x120f3e62,0x51a22a59,0x213b8fe4,0xb50e19f3,
  36 0x953816b6,0x611a9e9d,0x43508f58,0xc03b4ad7};
  37
  38 static int n = 31;
  39 static int i = 3;
  40 static int j = 0;
  41 static unsigned *x = init+1;
  42
  43 static unsigned lcg(unsigned x) {
  44         return 1664525*x + 1013904223;
  45 }
  46
  47 static void *savestate() {
  48         x[-1] = ((unsigned)n<<16)|(i<<8)|j;
  49         return x-1;
  50 }
  51
  52 static void loadstate(unsigned *state) {
  53         x = state+1;
  54         n = x[-1]>>16;
  55         i = (x[-1]>>8)&0xff;
  56         j = x[-1]&0xff;
  57 }
  58
  59 void srandom(unsigned seed) {
  60         int k;
  61
  62         i = 3;
  63         j = 0;
  64         x[0] = seed;
  65         for (k = 1; k < n; k++)
  66                 x[k] = lcg(x[k-1]);
  67 }
  68
  69 char *initstate(unsigned seed, char *state, size_t size) {
  70         void *old = savestate();
  71         if (size < 8)
  72                 return 0;
  73         x = (unsigned*)state + 1;
  74         n = size < 128 ? 0 : 31;
  75         srandom(seed);
  76         return old;
  77 }
  78
  79 char *setstate(char *state) {
  80         void *old = savestate();
  81         loadstate((unsigned*)state);
  82         return old;
  83 }
  84
  85 long random(void) {
  86         unsigned k;
  87
  88         if (n == 0) {
  89                 x[0] = lcg(x[0]);
  90                 return x[0]>>1;
  91         }
  92         x[i] += x[j];
  93         k = x[i]>>1;
  94         if (++i == n)
  95                 i = 0;
  96         if (++j == n)
  97                 j = 0;
  98         return k;
  99 }