/*2kadai02-2.c*/

/*ran1の定義*/
#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)

double ran1(long *idum)
{
int j;
long k;
static long iy=0;
static long iv[NTAB];
double temp;

if (*idum <= 0 || !iy){
if (-(*idum) < 1) *idum=1;
else *idum = -(*idum);
for (j=NTAB+7;j>=0;j--) {
k=(*idum)/IQ;
*idum=IA*(*idum-k*IQ)-IR*k;
if (*idum <0) *idum += IM;
if (j <NTAB) iv[j] = *idum;
}
iy=iv[0];
}
k=(*idum)/IQ;
*idum=IA*(*idum-k*IQ)-IR*k;
if (*idum < 0) *idum += IM;
j=iy/NDIV;
iy=iv[j];
iv[j] = *idum;
if ((temp=AM*iy) > RNMX) return RNMX;
else return temp;
}


/*gaussの定義*/
#include<math.h>
double gasdev(long *idum)
{
double ran1(long *idum);
static int iset=0;
static double gset;
double fac,rsq,v1,v2;

if (iset ==0){
do{
v1=2.0*ran1(idum)-1.0;
v2=2.0*ran1(idum)-1.0;
rsq=v1*v1+v2*v2;
}while (rsq >= 1.0 || rsq == 0.0);
fac=sqrt(-2.0*log(rsq)/rsq);
gset=v1*fac;
iset=1;
return v2*fac;
}else{
iset=0;
return gset;
}
}

#include<stdio.h>
#include<stdlib.h>
double ran1(long *);
double gasdev(long *);
main(){
long *idum;
int total;
int k;
double g,m,n,a,b,c,d1,d2,e1,e2,f1,f2,g1;
*idum=-1;
printf("乱数を何回発生させますか？:");
scanf("%d",&total);
m=0.0;
a=0.0;
for(k=0;k<total;k++)
{g=gasdev(idum);
m=m+g;
a=a+g*g;
}
m=m/(double)total;
n=m*m;
a=a/(double)total;
b=a-n;
c=sqrt(b);

*idum=-1;
d1=0.0;
e1=0.0;
f1=0.0;

for(k=0;k<total;k++)
{g=gasdev(idum);
g1=fabs(g);
if(g1<=c){d1=d1+1.0;}
if(g1<=2*c){e1=e1+1.0;}
if(g1<=3*c){f1=f1+1.0;}
}
d2=d1/(double)total*100;
e2=e1/(double)total*100;
f2=f1/(double)total*100;

/*各範囲における乱数の数を出力*/
printf("\n各範囲における乱数の数を出力\n");
printf("-σ=<g=<σの間には%lf個の乱数が入っています。\n",d1);
printf("-2σ=<g=<2σの間には%lf個の乱数が入っています。\n",e1);
printf("-3σ=<g=<3σの間には%lf個の乱数が入っています。\n",f1);

/*各範囲における乱数が入る確率を出力*/
printf("\n各範囲における乱数が入る確率を出力\n");
printf("-σ=<g=<σの間に乱数が入る確率は%lf％です。\n",d2);
printf("-2σ=<g=<2σの間に乱数が入る確率は%lf％です。\n",e2);
printf("-3σ=<g=<3σの間に乱数が入る確率は%lf％です。\n",f2);
}