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

long quickdirtyrand(unsigned long *seed);
int parkmillerrand(long *idum);

#define True 1

/*
  Exercise 1.2

  Tests of random number generators. To make the results
comparable, always set the initial seed=7.

   Write a piece of code which calculates its repeat interval
(by brute force). Hint: you can use a "double" (C) or "double
precision" (Fortran) variable to calculate the repeat
interval; this data type almost always nowadays has 15-digit
mantissa, so it is enough to calculate the interval of
any 32-bit integer (which only has up to 10 numbers).


a) 

ANSWER: 5432 (seed 1)
             (seed 7)

b) Repeat the same test for the Park-Miller minimal generator
given in Numerical Recipes (p. 279 in the C edition).

ANSWER: 2.14748e+09 (agrees with NR p. 279)


Note: if you carry out this exercise on a 64-bit computer
(e.g. any Alpha), you need to replace "long" by "int"
in the routines given.

------------------------------------------------------------

*/

main(int argc, char **argv)
{
  /* 
     This code assumes a mantissa in a double greater than the size of an int.
     If this is not true on your machine, you may end up in an infinite loop
  */

  double nbailout;   
  int n;
  int iloop;
  unsigned long first,new;
  unsigned long seedqd;
  long seedpm;


  if (sizeof(long) != 4) {
    printf("integer size for long should be 32 bits\n");
    exit(0);
  }


  /* Part a): test quickanddirty rand */

  seedqd=7;

  /* Get first true seed number in sequence */
  first=seedqd;

  n=0;
  nbailout=0.0;
  while(True) {
    n++; if (n==100000000) { printf("."); fflush(stdout); n=0; }
    nbailout+=1.0;
    new=quickdirtyrand(&seedqd);
    if (new==first) break;
    
    /* Check for bailout limit to avoid infinite looping */
    if (nbailout >= 1e10) break;
  }

  if (nbailout <1e10) {
    printf("\na) Random number generator repeats after %g iterations \n",nbailout);
  }
  else {
    printf("\na) Did not reach repeat interval by 1e10 iterations, bailed out\n");
  }

  /* Part b): test Park-Miller ANSI rand */

  seedpm=7;
  /* Get first number in sequence */
  first=parkmillerrand(&seedpm); 
  n=0;
  nbailout=0.0;
  while(True) {
    n++; if (n==100000000) { printf("."); fflush(stdout); n=0; }
    nbailout+=1.0;
    new=parkmillerrand(&seedpm);
    if (new==first) break;
    
    /* Check for bailout limit to avoid infinite looping */
    if (nbailout >= 1e10) break;
  }

  if (nbailout <1e10) {
    printf("\nb) Random number generator repeats after %g iterations\n",nbailout);
  }
  else {
    printf("\nDid not reach repeat interval by 1e10 iterations, bailed out\n");
  }
    
}

long quickdirtyrand(unsigned long *seed)
{

  long ia=106;
  long ic=1283;
  long im=6075;

  *seed = (*seed*ia + ic)%im;
  return((unsigned long) (*seed));
}




#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define MASK 123459876

/* Special version which returns ans without converting it to a float */
/* Assumes int is 32 bits */
int parkmillerrand(long *idum)
{
        long k;
        int ans;

        *idum ^= MASK;
        k=(*idum)/IQ;
        *idum=IA*(*idum-k*IQ)-IR*k;
        if (*idum < 0) *idum += IM;
        ans=*idum;
        *idum ^= MASK;
        return ans;
}