/* 確率分布の simulation
 * shige
 * sym1.c
 * 06/21 2001
 */
/* 1)	[0,1) を一様分布する y=drand48() から、ある確率分布に従う確率変
 * 	数を作る。そしてそれを使って中心極限定理を simulation する。
 * 2)	各確率分布に対して、x=F^{-1}(y) とすればよい。そのために、
 *	F^{-1}(y) を用意する必要がある。 
 */
/* gcc -Wall -O3 -o sym1 -DTEST0 sym1.c -lm ; ./sym1 > dat-0
 * gcc -Wall -O3 -o sym1 -DTEST1 sym1.c -lm
 * foreach N ( 1 3 5 10 20 50 ) 
 *   ./sym1 > dat-N$N
 * end
 */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/* 次のような密度関数 f(x):
 * f(x)=0 (x<a), =2(x-a)/(b-a)^2 (a<x<b), =0 (x>b)
 * (F(x)=0 (x<a), =(x-a)^2/(b-a)^2 (a<x<b), =1 (x>b))
 * ==> F^{-1}(y)=a+(b-a)*sqrt(y) (0<y<1)
 */
double Finv(double y,double a,double b)
{
    return a+(b-a)*sqrt(y);
}

#define DIV 1000	/* 度数分布を作るための分割数 */
#define A 5.0	/* x の最小値 */
#define B 6.0	/* x の最大値 */
#define SAMPLES 100000	/* 度数分布のサンプル数 */
#define class(x) (int)(floor(((x)-A)*DIV/(B-A)))
	/* 度数分布の階級を計算するマクロ */

/*#define N 50 */

#ifdef TEST1
/* 中心極限定理 */
int main(int argc,char *argv)
{
    int hist[DIV],j,index;
    int k,N;
    double x,dx;

    if(argc<2){ fprintf(stderr,"Usage: sym1 [N]\n"); return 1; }
    if((N=atoi(argv[1]))<=0){
	fprintf(stderr,"Specify the positive N.\n"); return 1;
    }
    for(j=1;j<=DIV;j++) hist[j-1]=0;
    for(j=1;j<=SAMPLES;j++){
	/* N 個の x の平均を取る */
	x=0.0;
	for(k=1;k<=N;k++) x+=Finv(drand48(),A,B);
	x/=N;
	index=class(x);
	if(index<0) index=0;
	else if(index>=DIV) index=DIV-1;
	hist[index]++;
    }
    dx=(B-A)/DIV;
    x=A;
    for(j=1;j<=DIV;j++){
	printf("%f %f\n",x,(double)hist[j-1]*DIV/SAMPLES);
	x+=dx;
    }
    return 0;
}
#endif

#ifdef TEST0
/* test 0 */
int main(void)
{
    int hist[DIV],j,index;
    double x,dx;

    for(j=1;j<=DIV;j++) hist[j-1]=0;
    for(j=1;j<=SAMPLES;j++){
	x=Finv(drand48(),A,B);
#ifdef DEBUG
fprintf(stderr,"%d: %f\n",j,x);
#endif
	index=class(x);
	if(index<0) index=0;
	else if(index>=DIV) index=DIV-1;
	hist[index]++;
    }
    dx=(B-A)/DIV;
    x=A;
    for(j=1;j<=DIV;j++){
	printf("%f %f\n",x,(double)hist[j-1]*DIV/SAMPLES);
	x+=dx;
    }
    return 0;
}
#endif