/*-------------------------------------------------------------

Copyright (C) 2000 Peter Clote
All Rights Reserved.

Permission to use, copy, modify, and distribute this
software and its documentation for NON-COMMERCIAL purposes
and without fee is hereby granted provided that this
copyright notice appears in all copies.


THE AUTHOR AND PUBLISHER MAKE NO REPRESENTATIONS OR
WARRANTIES ABOUT THE SUITABILITY OF THE SOFTWARE, EITHER
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE, OR NON-INFRINGEMENT. THE AUTHORS
AND PUBLISHER SHALL NOT BE LIABLE FOR ANY DAMAGES SUFFERED
BY LICENSEE AS A RESULT OF USING, MODIFYING OR DISTRIBUTING
THIS SOFTWARE OR ITS DERIVATIVES.

-------------------------------------------------------------*/



/*------------------------------------------------

	Program uniformDist.c

	This program illustrates the linear congruential pseudorandom
	number generator. To avoid overflow in multiplying large
	numbers, a multiplication routine (as in multiplying polynomials)
	is used.

	WARNING: a potential source of error lies in confusing long int
	and int. This program uses long integers, together with the
	conversion sequence
		%ld
	rather than
		%d
	If in main() one prints out the random numbers using %d conversion
	then negative numbers arise. If one uses the routines in this
	program within other programs, be sure to do a type conversion
		(int) random(max)
	or rewrite the routines using integers in place of long integers.


	USAGE: uniformDist <n>
	  <n> is the number of samples to compute mean and variance
 
-------------------------------------------------*/

#include <stdio.h>
#include <stdlib.h> // for RAND_MAX

#define DEBUG 0

/* see mult() and randomLong() for the 
   meaning of the following constants */
#define M 100000000
#define K 10000
#define A 31415926
#define B 1


/***************************************
	Idea: Global variable x is initialized to a seed.
	Program could be modified to allow user to define seed.
		x_{n+1} = (A * x_n + B) mod M
	M is 10^8, and K is 10^4, the square root of M.
***************************************/

long x = 7293641; /* initialize global x to seed */

long mult(long p, long q)
{
	long p0, p1, q0, q1;
	p0 = p % K; p1 = p/K;
	q0 = q % K; q1 = q/K;
	return( (((p0*q1 + p1*q0) % K)*K + p0*q0) % M);

}

long randomLong(long max)
{
	/* randomLong returns value in range 0..max-1  */
	long mult(long,long);  /* prototype of function called */

	x = (mult(A,x)+B) % M;  /* WARNING: x is global */
	return( x % max );

	/* Statistical properties of x % max are not as good as
		(long)((double)x/M * max)
	   but it is easier to write. */
}

/*
  main()
  
  approximate and print mean and variance for new defined and builtin
  random number generators. For further comparison print theoretical values.
*/
main(int argc, char *argv[])
{
        int i,n;
	double sum,sumSq,sample,sampleMean,sampleVariance;

	if (argc != 2) {
	   fprintf(stderr,"%s num\n",argv[0]);
	   exit(1); 
	}

	n = atoi(argv[1]);

	/* approximate and print mean and variance for randomLong (LCG
	   method random number generator) */

	sum = sumSq = 0;
	for (i=0;i<n;++i){

#if DEBUG
		/*----- for debugging ---*/
		printf("%ld\n",randomLong(M));
		sample = (double) randomLong(M)/M;
		printf("%lf\n",sample);
#endif

		/*----- gather statistics on unif generator ---*/
		sample = (double) randomLong(M)/M;
		sum += sample;
		sumSq += sample*sample;
	}
	sampleMean = sum/n;
	sampleVariance = sumSq/n - sampleMean*sampleMean;
	
	printf("\n\nUniform distribution in (0,1) - LCG method\n");
	printf("Mean\tVariance");
	printf("\n\n%lf\t%lf\n\n\n",sampleMean,sampleVariance);

	
	/* do the same for builtin random number generator */

	sum = sumSq = 0;  // reinitialize to 0
	for (i=0;i<n;++i){

#if DEBUG
		/*----- for debugging ---*/
		printf("%ld\n",rand());
		sample = (double) rand()/(RAND_MAX+1.0);
		printf("%lf\n",sample);
#endif

		/*----- gather statistics on unif generator ---*/
		sample = (double) rand()/(RAND_MAX+1.0);
		sum += sample;
		sumSq += sample*sample;
	}
	sampleMean = sum/n;
	sampleVariance = sumSq/n - sampleMean*sampleMean;
	
	printf("\n\nUniform distribution in (0,1) - built-in method.\n");
	printf("Mean\tVariance");
	printf("\n\n%lf\t%lf\n\n\n",sampleMean,sampleVariance);

	
	/* print theoretical values */
	
	printf("\n\nUniform distribution in (0,1) - theoretical values.\n");
	printf("Mean\tVariance");
	printf("\n\n%lf\t%lf\n\n\n", 1.0/2.0, 1.0/12.0);

}