/* L-6 MCS 572 Mon 23 Jan 2012 : sprng_normal.c
 * We use the uniformly distributed numbers by SPRNG
 * in the Box-Muller method to generate normally distributed numbers
 * with mean zero and standard deviation one. */

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

#define SIMPLE_SPRNG
#include "sprng.h"

int normal ( double *x, double *y );
/*
 * DESCRIPTION :
 *   Generates two independent normally distributed
 *   variables x and y, following Algorithm P (Knuth Vol 2), 
 *   the polar method is due to Box and Muller.
 *
 * ON ENTRY :
 *   x,y      two independent random variables,
 *            uniformly distributed in [0,1].
 *
 * ON RETURN : fail = normal(&x,&y)
 *   fail     if 1, then x and y are outside the unit disk,
 *            if 0, then x and y 
 *   x,y      two independent normally distributed variables. */

int check_normal ( int n );
/*
 * DESCRIPTION :
 *   Calls normal n times and the computes average mu,
 *   the standard deviation sigma, and 
 *   the number of samples in the interval [mu-sigma,mu+sigma].
 *   The integer on return is the number of samples generated. */

int main ( void )
{
   printf("normal variables with SPRNG ...\n");
   /* make new seed each time program is run */
   int seed = make_sprng_seed();  
   /* initialize the stream */
   init_sprng(0,seed,SPRNG_DEFAULT); 

   double x = sprng();
   double y = sprng();
   int fail = normal(&x,&y);
   if(fail != 0)
      printf("please try again...\n");
   else
   {
      printf("a normal random variable : %.15f\n",x);
      printf("a normal random variable : %.15f\n",y);
   }

   int n;
   printf("give number of samples : "); scanf("%d",&n);
   int nb = check_normal(n);
   printf("generated %d normal random numbers\n",nb);

   return 0;
}

int normal ( double *x, double *y )
{
   double s;

   *x = 2.0*(*x) - 1.0;
   *y = 2.0*(*y) - 1.0;
   s = (*x)*(*x) + (*y)*(*y);
   
   if(s >= 1.0)
      return 1;
   else
   {
      double ln_s = log(s);
      double rt_s = sqrt(-2.0*ln_s/s);
      *x = (*x)*rt_s;
      *y = (*y)*rt_s;
      return 0;
   }
}

int check_normal ( int n )
{
   double *data,x,y;
   int i,cnt = 0;
   
   data = (double*)calloc(2*n,sizeof(double));
   for(i=0; i<n; i++)
   {
      x = sprng(); y = sprng();
      if(normal(&x,&y) == 0)
      {
         data[cnt] = x; data[cnt+1] = y;
         cnt = cnt + 2;
      }
   }
   double mu = 0.0;
   for(i=0; i<cnt; i++) mu = mu + data[i];
   mu = mu/cnt;
   double sigma = 0.0;
   for(i=0; i<cnt; i++)
      sigma = sigma + (data[i]-mu)*(data[i]-mu);
   sigma = sqrt(sigma/cnt);
   printf("mu = %.15f, sigma = %.15f\n",mu,sigma);

   double a = mu-sigma;
   double b = mu+sigma;
   int nrmcnt = 0;
   for(i=0; i<cnt; i++)
      if(data[i] >= mu-sigma)
         if(data[i] <= mu+sigma) nrmcnt++;

   double r = ((double) nrmcnt)/cnt;
   printf("ratio of #samples in [%.2f,%.2f] : %.4f\n",a,b,r);
   return cnt;
}