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

#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 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 on input were not right,
 *            if 0, then x and y 
 *   x,y      two independent normally distributed variables. */

double map_to_normal ( double mu, double sigma, double x );
/*
 * DESCRIPTION :
 *   Given a normally distributed number x with average zero
 *   and standard deviation one, returns a normally distributed
 *   number y with average mu and standard deviation sigma. */

int normald ( double mu, double sigma, double *x, double *y );
/*
 * DESCRIPTION :
 *   Returns in x and y two normally distributed numbers
 *   with mean mu and standard deviation sigma. */

double mtbf ( int n, int m, double *mu, double *sigma );
/*
 * DESCRIPTION :
 *   Returns the expected life span of a composite product
 *   of m parts, where part i has expected life span mu[i],
 *   with standard deviation sigma[i], using n simulations. */

int main ( void )
{
   printf("MTBF problem with SPRNG ...\n");
   int seed = make_sprng_seed();  
   init_sprng(0,seed,SPRNG_DEFAULT); 

   int i,n,m;
   printf("Give number of components : ");
   scanf("%d",&m);
   double *mu = (double*)calloc(m,sizeof(double));
   double *sigma = (double*)calloc(m,sizeof(double));

   for(i=0; i<m; i++)
   {
      printf("average life span for part %d ? ",i);
      scanf("%lf",&mu[i]);
      printf("standard deviation for part %d ? ",i);
      scanf("%lf",&sigma[i]);
   }
   for(i=0; i<m; i++)
      printf("mu[%d] = %.3lf  sigma[%d] = %.3lf\n",i,mu[i],i,sigma[i]);
   printf("Give number of simulations : ");
   scanf("%d",&n);
 
   double lifespan = mtbf(n,m,mu,sigma);

   printf("expected life span : %.15lf\n",lifespan);

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

double map_to_normal ( double mu, double sigma, double x )
{
   return mu + sigma*x;
}

int normald ( double mu, double sigma, double *x, double *y )
{
   int fail = 1;
   while(fail == 1) 
   {
      *x = sprng();
      *y = sprng();
      fail = normal(x,y);
   }
   *x = map_to_normal(mu,sigma,*x);
   *y = map_to_normal(mu,sigma,*y);
   return 0;
}

double min ( double a, double b )
{
   if(a<b)
      return a;
   else
      return b;
}

double mtbf ( int n, int m, double *mu, double *sigma )
{
   int i,cnt=0;
   double s[n+1];
   do
   {
      normald(mu[0],sigma[0],&s[cnt],&s[cnt+1]);
      for(i=1; i<m; i++)
      {
         double x,y;
         normald(mu[i],sigma[i],&x,&y);
         s[cnt] = min(s[cnt],x);
         s[cnt+1] = min(s[cnt+1],y);
      }
      cnt = cnt + 2;
   } while (cnt < n);

   double sum = 0.0;
   for(i=0; i<cnt; i++) sum = sum + s[i]; 
   return sum/cnt;
}