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

unsigned long seed=7;

void ansistandardsrand(unsigned int seed_set);
double ansistandardrand(void);
double parkmillerrand(long *idum);

/*
  Exercise 1.3, see README

*/

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
  */

  int i;
  double x,y;
  FILE *fp;

  ansistandardsrand(7);

  fp=fopen("ansixy.dat","w");
  for (i=0;i<10000;i++) {
    x=ansistandardrand();
    y=ansistandardrand();
    fprintf(fp,"%g %g\n",x,y);
  }
  fclose(fp);

  seed=7;
  fp=fopen("parkmillerxy.dat","w");
  for (i=0;i<10000;i++) {
    x=parkmillerrand(&seed);
    y=parkmillerrand(&seed);
    fprintf(fp,"%g %g\n",x,y);
  }
  fclose(fp);


}

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

  seed = seed*a + c;
  return (double) ((seed/65536) % div)/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

double parkmillerrand(long *idum)
{
        long k;
        double ans;

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