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

unsigned long seed=7;

void ansistandardsrand(unsigned int seed_set);
int ansistandardrand(void);
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.

a) Test of why it is _not_ a good idea to use the system standard
random number generator. Below (also given in the lectures)
is the ANSI example standard random number generator.
  
   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).

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;
  int mya[10],myc[10];

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


  /* Part a): test ANSI rand */

  ansistandardsrand(7);

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

  n=0;
  nbailout=0.0;
  while(True) {
    n++; if (n==100000000) { printf("."); fflush(stdout); n=0; }
    nbailout+=1.0;
    ansistandardrand();
    new=seed;
    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 */

  seed=7;
  /* Get first number in sequence */
  first=parkmillerrand(&seed); 
  n=0;
  nbailout=0.0;
  while(True) {
    n++; if (n==100000000) { printf("."); fflush(stdout); n=0; }
    nbailout+=1.0;
    new=parkmillerrand(&seed);
    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");
  }
    
}

int ansistandardrand(void)
{
  long a=1103515245;
  long c=12345;
  long div=32768;

  seed = seed*a + c;
  return (unsigned int) (seed/65536) % div;
}

void ansistandardsrand(unsigned int seed_set)
{
  seed=seed_set;
}



#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;
}