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

static int ix = 1, iy = 1, iz = 1;  /* 1--30000の任意の整数 */

/*
 * 改良版乱数(Wichmann--Hillの乱数発生法)
 */
void init_rnd(int x, int y, int z)  /* x, y, z = 1..30000 */
{
  ix = x;
  iy = y;
  iz = z;
}

double Rand(void)  /* 0 <= rnd() < 1 */
{
  double r;

  ix = 171 * (ix % 177) -  2 * (ix / 177);
  iy = 172 * (iy % 176) - 35 * (iy / 176);
  iz = 170 * (iz % 178) - 63 * (iz / 178);
  if (ix < 0) ix += 30269;
  if (iy < 0) iy += 30307;
  if (iz < 0) iz += 30323;
  r = ix / 30269.0 + iy / 30307.0 + iz / 30323.0;
  while (r >= 1.0)
    r = r - 1.0;
  
  return r;
}

/*
 * ボックス＝ミューラー法
 */
double NormalRandom(double mu, double sigma)
{
  double t, u, r;
  static int sw = 0;
  static double rs;
  
  if (sw == 0) {
    t = sqrt(-2.0 * log(1.0 - Rand()));
    u = 2.0 * M_PI * Rand();
    r = t * cos(u);
    rs = t * sin(u);
    sw = 1;
  } else {
    r = rs;
    sw = 0;
  }
  
  return sigma * r + mu;
}

int main(void)
{
  int i, n, ix, *histo;
  double mu, sigma, x, s1 = 0.0, s2 = 0.0;
  double *array;
  
  init_rnd((time(NULL) % 30000) + 1, ((time(NULL) + 23456) % 30000) + 1, ((time(NULL) + 45678) % 30000) + 1);  /* 1..30000 の任意の整数で初期化. */
//  printf("個数 n = ");
//  scanf("%d", &n);
  n = 10000;
//  printf("平均 μ = ");
//  scanf("%lf", &mu);
  mu = 1.5;
//  printf("分散 σ = ");
//  scanf("%lf", &sigma);
  sigma = 0.5;
  
  if ((histo = (int *)malloc(sizeof(int) * n)) == NULL)
    exit(1);
  if ((array = (double *)malloc(sizeof(double) * n)) == NULL)
    exit(1);
  
  /*
   * 乱数はarray[]に保存する
   * 度数分布を同時に求める
   */
  for (i = 0; i < n; i++) {
    x = array[i] = NormalRandom(mu, sigma);
    ix = (int)floor(2.0 * x) + 10;
    if (ix >= 0 && ix < 20)
      histo[ix]++;
    s1 += x;
    s2 += x * x;
  }

  /*
   * 平均と標準偏差の計算
   */
  for (i = 0; i < 20; i++)
    printf("%4.1f -- %4.1f: %5.1f%%\n", 0.5 * (i - 10), 0.5 * (i - 9), 100.0 * histo[i] / n);
  s1 /= n;
  s2 = sqrt((s2 - n * s1 * s1) / (n - 1));
  printf("平均 %g  標準偏差 %g\n", s1, s2);
  
  free(array);
  free(histo);
  
  return 0;
}