#include <limits.h>

#include "random.h"

/* assumes 32-bit ints in rnd() */
/* assumes IEEE doubles in frnd() */

#define N 63 /* number of unsigned longs of state */

static unsigned long int state[N]
 = { 0x9a319039, 0x32d9c024, 0x9b663182, 0x5da1f342,
     0xde3b81e0, 0xdf0a6fb5, 0xf103bc02, 0x48f340fb,
     0x7449e56b, 0xbeb1dbb0, 0xab5c5918, 0x946554fd,
     0x8c2e680f, 0xeb3d799f, 0xb11ee0b7, 0x2d436b86,
     0xda672e2a, 0x1588ca88, 0xe369735d, 0x904f35f7,
     0xd7158fd6, 0x6fa6f051, 0x616e6b96, 0xac94efdc,
     0x36413f93, 0xc622c298, 0xf5a42ab8, 0x8a88d77b,
     0xf5ad9d0e, 0x8999220b, 0x27fb47b9 };
static int sp = 0;
static int h1 = 0;
static int h2 = 1;
static int stirred = 0;

static unsigned long int r(void)
{
 register int rv;

 rv = state[h2] += state[h1];
 if (++h1 >= N) h1 = 0;
 if (++h2 >= N) h2 = 0;
 return(rv);
}

static void stir(void)
{
 register int n;

 if (stirred) return;
 h2 = (h1 + h2) % N;
 for (n=1024;n>0;n--) r();
 stirred = 1;
}

int rndint(void)
{
 stir();
 return(((unsigned int)r())>>1);
}

MUFQUAD rndquad(void)
{
 MUFQUAD q;
 int i;

 stir();
 for (i=0;i<sizeof(q);i++) q = (q ^ r()) << CHAR_BIT;
 return(q);
}

unsigned int rnd(register unsigned int n)
{
 register unsigned int i;

 i = r();
 if (n & (n-1)) return(i%n);
 if (n & 0x0000ffff) i >>= 16;
 if (n & 0x00ff00ff) i >>= 8;
 if (n & 0x0f0f0f0f) i >>= 4;
 if (n & 0x33333333) i >>= 2;
 if (n & 0x55555555) i >>= 1;
 return(i);
}

double frnd(void)
{
 double d;

 stir();
 ((unsigned long int *)&d)[0] = r();
 ((unsigned long int *)&d)[1] = r();
 ((unsigned char *)&d)[0] = 0x3f;
 ((unsigned char *)&d)[1] |= 0xf0;
 return(d-1);
}

void seed_random(const void *data, register unsigned int nb)
{
 register const unsigned char *dp;
 register unsigned long int t;
 register int x;

 dp = data;
 x = sp;
 for (;nb>0;nb--)
  { t = state[x] ^ *dp++;
    state[x] = (t & 1) ? ((t>>1) ^ 0xedb88320) : (t>>1);
    if (++x >= N) x = 0;
  }
 sp = x;
}