/* Test normal.   v. 020327 */

#include <math.h>

  double seed = 86329.0;

int main()
{

  int    i;

  double normal1(),
         val, sum=0,
         NumValues = 5000;

  /* Generate normals */
     for ( i = 0; i<NumValues; i++ )
       {
       val =  normal1( 0.75, 4.5 );
       sum += val;
       printf( "%6.5f\n", val );
       }
     printf( "Avg: %6.5f\n", sum/NumValues );
}

double normal1( double mean, double sd )
  {
  int i;
  double sum = -6.0;

  double random();

  for ( i=0; i<12; i++ )
    sum += random( &seed );
  return( mean + sd*sum );

  }

/* Park & Miller function (see CACM Oct. 88) real version */
/* of their generator.                                    */
double random( seed )
  double *seed;
{
	double  a = 16807.0,            /* 7^5 */
	        m = 2147483647.0,       /* 2^31 - 1 */
	        q = 127773.0,           /* m / a (int divide) */
	      	r = 2836.0,		/* m % a */
		lo, hi, test;

	hi = floor( *seed / q );
	lo = *seed - q * hi;
	test = a * lo - r * hi;
	*seed = test > 0.0 ? test : test + m;
	return( (double ) *seed / m );
}