/* starasc.c */

#include <stdio.h>
#include <stdlib.h>
#include <limits.h>
#include <tgmath.h>
#include <values.h>
#include "starasc.h"


int main ( void )
{
  int i, n = 0;
  double r;
  for (i = 0; i < 100; i ++) {
    r = gaussrand();
    if (abs(r) < 1.0) ++ n;
    /*
    printf("%e\n", r);
    */
  }
  printf("%d/%d\n", n, i);
  exit(0);
}

/*******************************************************************/

/* http://c-faq.com/lib/gaussian.html */

double gaussrand()
{
  
  static double V1, V2, S;
  static int phase = 0;
  double X;
  
  if (phase == 0) 
    {
      do
	{
	  double U1 = (double) rand() / RAND_MAX;
	  double U2 = (double) rand() / RAND_MAX;

	  V1 = 2 * U1 - 1;
	  V2 = 2 * U2 - 1;
	  S = V1 * V1 + V2 * V2;
	} while ( S >= 1 || S == 0);
      
      X = V1 * sqrt(-2 * log(S) / S);
    }
  else
    X = V2 * sqrt(-2 * log(S) / S);
  
  phase = 1 - phase;
  
  return X;
}