commit last state

[ppam-mpi.git] / code / src / Util / rng.c

// ISAAC PRNG
// Initial code by Bob Jenkins, March 1996 - modified in 2008
// Reference: http://burtleburtle.net/bob/rand/isaacafa.html
// Further slightly modified by Benjamin Auder, 2013

#include "Util/rng.h"
#include <stdlib.h>
#include <time.h>

// Internal state
static uint32_t mm[256];
static uint32_t aa, bb, cc;

// initialize the (pseudo-)random generator
void init_rng(int flag)
{
	// 'bootstrap' isaac with basic bad PRNG
	srand(time(NULL));

	aa = rand() & 0xffffffff;
	bb = rand() & 0xffffffff;
	cc = rand() & 0xffffffff;

	uint32_t a, b, c, d, e, f, g, h, i;
	// The golden ratio:
	a = b = c = d = e = f = g = h = 0x9e3779b9;
	// Scramble it
	for (i = 0; i < 4; ++i)
		mix(a, b, c, d, e, f, g, h);

	// Fill in mm[] with messy stuff (the seed)
	uint32_t seedArray[256];
	for (i = 0; i < 256; ++i) mm[i] = seedArray[i] = rand() & 0xffffffff;
	for (i = 0; i < 256; i += 8)
	{
		if (flag)
		{
			// Use all the information in the seed
			a += seedArray[i  ];
			b += seedArray[i + 1];
			c += seedArray[i + 2];
			d += seedArray[i + 3];
			e += seedArray[i + 4];
			f += seedArray[i + 5];
			g += seedArray[i + 6];
			h += seedArray[i + 7];
		}
		mix(a, b, c, d, e, f, g, h);
		mm[i  ] = a;
		mm[i + 1] = b;
		mm[i + 2] = c;
		mm[i + 3] = d;
		mm[i + 4] = e;
		mm[i + 5] = f;
		mm[i + 6] = g;
		mm[i + 7] = h;
	}

	if (flag)
	{
		// Do a second pass to make all of the seed affect all of mm
		for (i = 0; i < 256; i += 8)
		{
			a += mm[i  ];
			b += mm[i + 1];
			c += mm[i + 2];
			d += mm[i + 3];
			e += mm[i + 4];
			f += mm[i + 5];
			g += mm[i + 6];
			h += mm[i + 7];
			mix(a, b, c, d, e, f, g, h);
			mm[i  ] = a;
			mm[i + 1] = b;
			mm[i + 2] = c;
			mm[i + 3] = d;
			mm[i + 4] = e;
			mm[i + 5] = f;
			mm[i + 6] = g;
			mm[i + 7] = h;
		}
	}
}

// return a (pseudo-)random integer
uint32_t get_rand_int()
{
	// TODO: register variables ? (x,y,i)
	uint32_t x, y;
	static uint32_t i = 0;

	if (i == 0)
	{
		cc = cc + 1;  // cc just gets incremented once per 256 results
		bb = bb + cc; // then combined with bb
	}

	x = mm[i];
	switch (i % 4)
	{
	case 0:
		aa = aa^(aa << 13);
		break;
	case 1:
		aa = aa^(aa >> 6);
		break;
	case 2:
		aa = aa^(aa << 2);
		break;
	case 3:
		aa = aa^(aa >> 16);
		break;
	}

	// NOTE: bits 2..9 are chosen from x but 10..17 are chosen
	//       from y.  The only important thing here is that 2..9 and 10..17
	//       don't overlap.  2..9 and 10..17 were then chosen for speed in
	//       the optimized version (rand.c) */
	// See http://burtleburtle.net/bob/rand/isaac.html
	// for further explanations and analysis.

	aa        = mm[(i + 128) % 256] + aa;
	mm[i] = y = mm[(x >> 2) % 256]  + aa + bb;
	bb        = mm[(y >> 10) % 256] + x;

	i = (i + 1) % 256;

	return bb;
}

// return a (pseudo-)random real number in [0,1]
Real get_rand_real()
{
	return (Real) (INV_RANDMAX * (double) get_rand_int());
}