/* 
 * simulate normals via rejection 
 *
 * Math 5040
 */

/* Compile line:
 * gcc normal.c -o normal -lm    
 */

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

int main(void) {
	
	float U,V,W,X,Y;
	int i;
	int n=0;
	
	srand(time(NULL));  /* set seed for pseudo random number generator */
	
	for(i=0;i<100;i++) 
	{
	while(1) {
	  U = rand()/(float) RAND_MAX; /* U uniform */
	  V = rand()/(float) RAND_MAX; /* V uniform */
	  Y = -log(1-V); /* Y exponential */
	  /* check if U < f(Y)/cg(Y) */
	  if (U < exp(-(Y-1)*(Y-1)*0.5)) {
	    W = rand()/(float) RAND_MAX; /* W uniform */
	    /* Determine sign of X */
	    if (W < 0.5) {X = -Y;}
	    else {X = Y;}
	    break;
	  }
	}
	if ((X > -1)&&(X < 1)) n = n + 1; /* increment if |X| < 1 */
	}
	printf("fraction of 100 independent normals in (-1,1) is %.2f.\n",
	       n/(float) 100);
	
	return(0);
}