#pragma warning (disable:4996)
#include<stdio.h>
#include<math.h>

#define M_PI (double)3.14159265358979323846e0
#define NN (int)1000

int TNHSS;	//Total Number of the High School Student
int TNUE;	//Total Number of the University Examinee
int TNUE_C;	//Calculated Total Number of the University Examinee

double DVHSS[NN];	//Deviation Value of the High School Student
int RO_HSDV[NN];	//Rank Order of the High School Student (HSDV)
double ARUE_HSDV[NN];//Application Rate of the University Examinee (HSDV)
int ROUE_HSDV[NN];	//Rank Order of the University Examinee (HSDV)
double DVUE[NN];	//Deviation Value of the University Examinee
int ROUE_UEDV[NN];	//Rank Order of the University Examinee (UEDV)

void INIT(void);
double NDF5010(double bDV);	//Normal Distribution Function, MV is 50, SD is 10
double ARUE(double bDV);	//Applicant Rate of the University Examinee
double INTEGRAND(int mode,double bDV);//integrand function
double SIMPSON(int mode,double bDDV,double bLDV,double bUDV);
//integration of NDF5010 with Simpson's rule

int main(void)
{
	int i;
	double DELTA_DV,LOWER_DV,UPPER_DV;
	double buf0;
	FILE *fp;

	INIT();

	LOWER_DV=0.0e0;
	UPPER_DV=1.0e2;
	DELTA_DV=pow(2.0e0,-4.0e0);

	TNHSS=1055807;	//total number of the high school student
	TNUE=631267;	//total number of the university examinee
	TNUE_C=(int)SIMPSON(1,DELTA_DV,LOWER_DV,UPPER_DV);
	//total number of the university examinee

	for(i=0;i<NN;i++){
		buf0=UPPER_DV-(double)i*1.0e-1;
		DVHSS[i]=buf0;
		RO_HSDV[i]=(int)SIMPSON(0,DELTA_DV,buf0,UPPER_DV);
		ARUE_HSDV[i]=1.0e2*ARUE(buf0);
		ROUE_HSDV[i]=(int)SIMPSON(1,DELTA_DV,buf0,UPPER_DV);
		DVUE[i]=buf0;
		ROUE_UEDV[i]=(int)SIMPSON(2,DELTA_DV,buf0,UPPER_DV);
	}

	fp=fopen("log-data.txt","w");
	fprintf(fp,"Basic Parameter DELTA_DV : %12.4e\n",DELTA_DV);
	fprintf(fp,"                LOWER_DV : %12.4e\n",LOWER_DV);
	fprintf(fp,"                UPPER_DV : %12.4e\n\n",UPPER_DV);
	fprintf(fp,"Total Number of the High School Student               : %8d\n",TNHSS);
	fprintf(fp,"                the University Examinee               : %8d\n",TNUE);
	fprintf(fp,"                the University Examinee (calculation) : %8d\n\n",TNUE_C);
	fprintf(fp,"i,DVHSS,RO_HSDV,ARUE_HSDV[%],ROUE_HSDV,DVUE,ROUE_UEDV\n");
	for(i=0;i<NN;i++){
		fprintf(fp,"%d,%8.2f,%d,%8.2f,%d,%8.2f,%d\n",
			i,DVHSS[i],RO_HSDV[i],ARUE_HSDV[i],ROUE_HSDV[i],DVUE[i],ROUE_UEDV[i]);
	}
	fclose(fp);

	return 0;
}

void INIT(void)
{
	int i;

	for(i=0;i<NN;i++){
		DVHSS[i]=0.0e0;
		RO_HSDV[i]=0;
		ARUE_HSDV[i]=0.0e0;
		ROUE_HSDV[i]=0;
		DVUE[i]=0.0e0;
		ROUE_UEDV[i]=0;
	}

	return;
}

double NDF5010(double bDV)
{
	double value,bMV,bSD;
	double buf0,buf1;

	bMV=5.0e1;
	bSD=1.0e1;
	buf0=1.0e0/(sqrt(2.0e0*M_PI)*bSD);
	buf1=-((bDV-bMV)*(bDV-bMV))/(2.0e0*bSD*bSD);
	value=buf0*exp(buf1);

	return value;
}

double ARUE(double bDV)
{
	double value;

	value=1.0e-2*(1.015e-1*bDV*bDV-6.09e0*bDV+1.085e2);
	if(bDV<3.0e1){value=1.7e-1;}
	if(bDV>=5.8e1){value=1.0e0;}

	return value;
}

double INTEGRAND(int mode,double bDV)
{
	double value;
	double buf0,buf1;

	switch (mode)
	{
	case 0:
		buf0=NDF5010(bDV);
		value=TNHSS*buf0;
		break;
	case 1:
		buf0=ARUE(bDV);
		buf1=NDF5010(bDV);
		value=TNHSS*buf0*buf1;
		break;
	case 2:
		buf0=NDF5010(bDV);
		value=TNUE_C*buf0;
		break;
	default:
		value=-1.0e0;
		break;
	}

	return value;
}

double SIMPSON(int mode,double bDDV,double bLDV,double bUDV)
{
	int i,bn;
	double value;
	double buf0,buf1,buf2,buf3;

	if((bDDV<=0.0e0)||(bUDV<=bLDV)){return 0.0e0;}
	bn=(int)((bUDV-bLDV)/bDDV);

	if(bn==0){
		buf0=5.0e-1*(bLDV+bUDV);
		value=INTEGRAND(mode,bLDV)+4.0e0*INTEGRAND(mode,buf0)+INTEGRAND(mode,bUDV);
		value*=(bUDV-bLDV)/6.0e0;
	}else{
		value=0.0e0;
		buf0=5.0e-1*bDDV;
		for(i=1;i<=bn;i++){
			buf1=bLDV+(double)(i-1)*bDDV;
			buf2=buf1+buf0;
			buf3=buf1+bDDV;
			value+=INTEGRAND(mode,buf1)+4.0e0*INTEGRAND(mode,buf2)+INTEGRAND(mode,buf3);
		}
		value*=bDDV/6.0e0;
		buf1=buf3;
		buf2=5.0e-1*(buf1+bUDV);
		buf3=bUDV;
		buf0=INTEGRAND(mode,buf1)+4.0e0*INTEGRAND(mode,buf2)+INTEGRAND(mode,buf3);
		buf0*=(bUDV-buf1)/6.0e0;
		value+=buf0;
	}

	return value;
}