/* qadrng.h
 * Quick And Dirty Random Number Generator
 * Based on the quick and dirty generator in "Numerical Recipes"
 * Parallelized using a simple leapfrog algorithm
 */

#ifndef _QADRNG_H_
#define _QADRNG_H_ 1

#include <stdio.h>
#include <stdlib.h>

#define QAD_MODULUS      4294967296.0
#define QAD_MULTIPLIER   1664525L
#define QAD_INCREMENT    1013904223L
#define QAD_MAX_PROCS    256

static unsigned long qad_seed[QAD_MAX_PROCS];
static unsigned long qad_mult;
static unsigned long qad_incr;
static int qad_procs = 1;

static void qad_init (int processes, unsigned long seed) {
     int i;

     if (processes < 1 || processes > QAD_MAX_PROCS) {
	  fprintf(stderr, "qad_init: numprocs = %d ", processes);
	  fprintf(stderr, " must be in range [1, %d]\n", QAD_MAX_PROCS);
	  exit(1);
     } else qad_procs = processes;

     for (i = 0; i < processes; i++) {
	  qad_seed[i] = seed;
	  seed = QAD_MULTIPLIER * seed + QAD_INCREMENT;
     }

     qad_mult = QAD_MULTIPLIER;
     qad_incr = QAD_INCREMENT;
     for (i = 1; i < processes; i++)
	  qad_mult *= QAD_MULTIPLIER;
     qad_incr = QAD_INCREMENT * (qad_mult - 1) / (QAD_MULTIPLIER - 1);
}

static unsigned long qad_get_seed (void) {
     return qad_seed[0];
}

static double qad_rand (int process) {
     qad_seed[process] = qad_mult * qad_seed[process] + qad_incr;
     return qad_seed[process] / QAD_MODULUS;
}

#define QAD_RAND(process) \
     ( (qad_seed[(process)] = \
          qad_mult * qad_seed[(process)] + qad_incr ) / QAD_MODULUS )

#endif /* _QADRNG_H_ */