/*    This program by D E Knuth is in the public domain and freely copyable
 *    AS LONG AS YOU MAKE ABSOLUTELY NO CHANGES!
 *    It is explained in Seminumerical Algorithms, 3rd edition, Section 3.6
 *    (or in the errata to the 2nd edition --- see
 *        http://www-cs-faculty.stanford.edu/~knuth/taocp.html
 *    in the changes to Volume 2 on pages 171 and following).              */

/*    N.B. The MODIFICATIONS introduced in the 9th printing (2002) are
      included here; there's no backwards compatibility with the original. */

/*    This version also adopts Brendan McKay's suggestion to
      accommodate naive users who forget to call ran_start(seed).          */

/*    If you find any bugs, please report them immediately to
 *                 taocp@cs.stanford.edu
 *    (and you will be rewarded if the bug is genuine). Thanks!            */

/************ see the book for explanations and caveats! *******************/
/************ in particular, you need two's complement arithmetic **********/

#define KK 100                     /* the long lag */
#define LL  37                     /* the short lag */
#define MM (1L<<30)                 /* the modulus */
#define mod_diff(x,y) (((x)-(y))&(MM-1)) /* subtraction mod MM */

long ran_x[KK];                    /* the generator state */

#ifdef __STDC__
	void ran_array(long aa[], int n)
#else
	void ran_array(aa, n)   /* put n new random numbers in aa */
	long * aa;  /* destination */
	int n;      /* array length (must be at least KK) */
#endif
{
	register int i, j;

	for (j = 0; j < KK; j++) {
		aa[j] = ran_x[j];
	}

	for (; j < n; j++) {
		aa[j] = mod_diff(aa[j - KK], aa[j - LL]);
	}

	for (i = 0; i < LL; i++, j++) {
		ran_x[i] = mod_diff(aa[j - KK], aa[j - LL]);
	}

	for (; i < KK; i++, j++) {
		ran_x[i] = mod_diff(aa[j - KK], ran_x[i - LL]);
	}
}

/* the following routines are from exercise 3.6--15 */
/* after calling ran_start, get new randoms by, e.g., "x=ran_arr_next()" */

#define QUALITY 1009 /* recommended quality level for high-res use */
long ran_arr_buf[QUALITY];
long ran_arr_dummy = -1, ran_arr_started = -1;
long * ran_arr_ptr = &ran_arr_dummy; /* the next random number, or -1 */

#define TT  70   /* guaranteed separation between streams */
#define is_odd(x)  ((x)&1)          /* units bit of x */

#ifdef __STDC__
	void ran_start(long seed)
#else
	void ran_start(seed)    /* do this before using ran_array */
	long seed;            /* selector for different streams */
#endif
{
	register int t, j;
	long x[KK + KK - 1];          /* the preparation buffer */
	register long ss = (seed + 2) & (MM - 2);

	for (j = 0; j < KK; j++) {
		x[j] = ss;                    /* bootstrap the buffer */
		ss <<= 1;

		if (ss >= MM) {
			ss -= MM - 2; /* cyclic shift 29 bits */
		}
	}

	x[1]++;              /* make x[1] (and only x[1]) odd */

#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wcomma"

	for (ss = seed & (MM - 1), t = TT - 1; t; ) {
		for (j = KK - 1; j > 0; j--) {
			x[j + j] = x[j], x[j + j - 1] = 0; /* "square" */
		}

		for (j = KK + KK - 2; j >= KK; j--)
			x[j - (KK - LL)] = mod_diff(x[j - (KK - LL)], x[j]),
							   x[j - KK] = mod_diff(x[j - KK], x[j]);

		if (is_odd(ss)) {              /* "multiply by z" */
			for (j = KK; j > 0; j--) {
				x[j] = x[j - 1];
			}

			x[0] = x[KK];          /* shift the buffer cyclically */
			x[LL] = mod_diff(x[LL], x[KK]);
		}

		if (ss) {
			ss >>= 1;
		} else {
			t--;
		}
	}

#pragma clang diagnostic pop

	for (j = 0; j < LL; j++) {
		ran_x[j + KK - LL] = x[j];
	}

	for (; j < KK; j++) {
		ran_x[j - LL] = x[j];
	}

	for (j = 0; j < 10; j++) {
		ran_array(x, KK + KK - 1); /* warm things up */
	}

	ran_arr_ptr = &ran_arr_started;
}

#define ran_arr_next() (*ran_arr_ptr>=0? *ran_arr_ptr++: ran_arr_cycle())
long ran_arr_cycle() {
	if (ran_arr_ptr == &ran_arr_dummy) {
		ran_start(314159L);    /* the user forgot to initialize */
	}

	ran_array(ran_arr_buf, QUALITY);
	ran_arr_buf[KK] = -1;
	ran_arr_ptr = ran_arr_buf + 1;
	return ran_arr_buf[0];
}

/* Tweaked to include as a library - Fletcher T. Penney */
/*#include <stdio.h>
int main()
{
  register int m; long a[2009];
  ran_start(310952L);
  for (m=0;m<=2009;m++) ran_array(a,1009);
  printf("%ld\n", a[0]);             *//* 995235265 */
/*  ran_start(310952L);
  for (m=0;m<=1009;m++) ran_array(a,2009);
  printf("%ld\n", a[0]);             *//* 995235265 */
/*  printf("%ld\n",ran_arr_next());
  return 0;
} */

long ran_num_next(void) {
	return ran_arr_next();
}