/* Box-Muller法による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;
  int i, j, nin;
  double r1, r2, fact, x1, x2;
  float x[2000001]; /* 乱数を保存する配列 */
  int im,m[100];   /* 乱数の分布 */
  double xg, nexact;
  double pi=3.141592653589793238462643;

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

  nin=0;          /* 採用数 */
  i=0;

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

  while( nin < nrandom ){
      j=rand();
      r1=j/(RAND_MAX+1.0);  /* [0,1]の一様乱数 */
      j=rand();
      r2=j/(RAND_MAX+1.0);  /* [0,1]の一様乱数 */

      fact=sqrt(-2.0*log(r1));
 
      x1=fact*cos(2.0*pi*r2);
      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;

      x2=fact*sin(2.0*pi*r2);
      nin=nin+1;
      x[nin]=x2;

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

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