/* 棄却法によるGauss乱数の生成 */

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

double g(double x)
{
  double pi=3.141592653589793238462643;
  return( exp(-x*x/2.0)/sqrt(2.0*pi) );
} 

main()
{
  unsigned int seed;
  int nrandom;
  double limit;
  int i, j, nin, ntrial;
  double a1;
  double x1, y1;
  float x[1000000]; /* 乱数を保存する配列 */
  int im,m[100];   /* 乱数の分布 */
  double xg, nexact;

  FILE *output;
  
  output = fopen("nm.dat","w");


  printf("Input: Number of random numbers:\n");
  scanf("%u", &nrandom);
  printf("Input: seed\n");
  scanf("%u", &seed);
/*  printf("seed=%u\n",seed); */

  srand(seed);

  limit=3.5;      /* 乱数の範囲 ( -limit < x < limit ) */

  ntrial=0;       /* 試行回数 */
  nin=0;          /* 採用数 */
  i=0;

  for(i=0; i < 100; i++){
    m[i]=0;
  }

  while( nin < nrandom ){
      j=rand();
      a1=j/(RAND_MAX+1.0);  /* [0,1]の一様乱数 */
      x1=-limit+2*limit*a1; /* [-limit,limit]の一様乱数 */
      j=rand();
      a1=j/(RAND_MAX+1.0);  /* [0,1]の一様乱数 */
      y1=a1;
      ntrial=ntrial+1;

      if( y1 < exp(-x1*x1/2.0) ) {
	nin=nin+1;   /* 採用 */
	x[nin]=x1;
/* 	printf("x=%.5f\n",x[nin]); */
	im=(int) floor((10.0*x[nin]));
/* 	printf("x=%f\tim=%d\n",x[nin],im); */
	m[im+50]=m[im+50]+1;
      }
    }

  printf("ntrial=%u\n", ntrial);

  for(im=0; im<100; im++){
    xg=(im-50)/10.0;
    nexact=nrandom*g(xg)*0.1;
/*     printf("%.2f  %d %.2f\n",xg,m[im],nexact); */
    fprintf(output, "%.2f  %d %.2f \n",xg,m[im],nexact);
  }

  return;
}