/*****************************/
/*     乱数の発生            */
/*          coded Y.Suganuma */
/*****************************/
#include <stdio.h>
#include <stdlib.h>

double exp_d(double);
double norm_d(double, double);
double unifm_d();

main()
{
	double a, b, m, r, mean, sd;
	int i1, n = 10000, seed;

	printf("乱数の初期値を入力してください ");
	scanf("%d", &seed);

	srand(seed);   // 乱数の初期設定
/*
          [a b] 区間の一様乱数
*/
	mean = 0.0;
	sd   = 0.0;
	a    = 2.0;
	b    = 4.0;

	for (i1 = 0; i1 < n; i1++) {
		r     = a + (b - a) * unifm_d();
		mean += r;
		sd   += r * r;
	}

	mean = mean / n;
	sd   = sqrt(fabs(n /(n - 1.0) * (sd / n - mean * mean)));

	printf("一様分布\n");
	printf("   平均 %f 標準偏差 %f\n", mean, sd);
/*
          平均 m の指数分布
*/
	mean = 0.0;
	sd   = 0.0;
	m    = 0.5;

	for (i1 = 0; i1 < n; i1++) {
		r     = exp_d(m);
		mean += r;
		sd   += r * r;
	}

	mean = mean / n;
	sd   = sqrt(fabs(n /(n - 1.0) * (sd / n - mean * mean)));

	printf("指数分布\n");
	printf("   平均 %f 標準偏差 %f\n", mean, sd);
/*
          標準正規分布
*/
	mean = 0.0;
	sd   = 0.0;

	for (i1 = 0; i1 < n; i1++) {
		r     = norm_d(0.0, 1.0);
		mean += r;
		sd   += r * r;
	}

	mean = mean / n;
	sd   = sqrt(fabs(n /(n - 1.0) * (sd / n - mean * mean)));

	printf("標準正規分布\n");
	printf("   平均 %f 標準偏差 %f\n", mean, sd);

	return 0;
}

/**********************************/
/*     [0, 1]区間の一様乱数の発生 */
/*          rerutn : 乱数         */
/**********************************/
double unifm_d()
{
   double x;

   while ((x = (double)rand() / RAND_MAX) == 0.0)
      ;

   return x;
}

/*********************************/
/*     平均値 m の指数乱数の発生 */
/*          m : 平均             */
/*********************************/
#include <math.h>

double exp_d(double m)
{
	return -m * log(unifm_d());
}

/***************************************/
/*     正規分布変量の発生              */
/*          m : 平均                   */
/*          s : 標準偏差               */
/*               return : 正規分布変量 */
/***************************************/
double norm_d(double m, double s)
{
	double x;
	int i1;

	x = 0.0;

	for (i1 = 0; i1 < 12; i1++)
		x += unifm_d();

	x = s * (x - 6.0) + m;

	return x;
}