int comp_doub_ascend(const void * a, const void * b) {
if (*(double*)a > *(double*)b) return 1;
else if (*(double*)a < *(double*)b) return -1;
else return 0;
}
double erfinv(double x){
int k,m,n = 20;
double s = 0,y = 0,q;
double c[n];
c[0] = 1;
for(k=1;k<n;k++){
s = 0;
for(m=0;m<k;m++)s += (c[m]*c[k-1-m])/((m+1)*(2*m+1));
c[k] = s;
}//k
q = sqrt(_PI)*x/2.;
for(k=0;k<n;k++) y += ((c[k])/(2*k+1))*pow(q,2*k+1);
return y;
}//erfinv()
double erf(double x){
int i,n = 200;
double t = 0,s = 0,dt = x/double(n);
for(i=0;i<n;i++){
s += exp(-t*t)*dt;
t += dt;
}
return s * 2./sqrt(_PI);
}//erf()
void exp_vars(int n,double *y,double u,unsigned long long seed){
// array of random numbers from an exonentual distribution
// URAND random number algorithm
// p represents random probabilities >0 <1
double p;
int i;
static unsigned long long s = 0,IY = 1;
if(s != seed){IY = seed; s = seed;}
unsigned long long m = pow(2,32)-1;
unsigned long long IC = 2*floor(m*(.5 - sqrt(3)/6.)) + 1;
unsigned long long IA = 8*m*ceil(atan(1.)/8.) + 5;
IY = IY * IA + IC;
for(i=0;i<n;i++){
p = (double)(IY%(m+1))/(double)(m+1);
IY = IY * IA + IC;
// quantile - probability to variable for a given mean
y[i] = -1.* log(1. - p) * u;
}
}//exp_vars()
void norm_vars(long long n,double *y,double u,double sigma,unsigned long long seed){
// array of random numbers from a normal distribution
// URAND random number algorithm
// p represents random probabilities >0 <1
double q,p,einv,k1;
int i,k,nc = 20;
static unsigned long long s = 0,IY = 1;
if(s != seed){IY = seed; s = seed;}
// random number gnerator constants
unsigned long long m = pow(2,32);
unsigned long long IC = 2*floor(m*(.5 - sqrt(3)/6.)) + 1;
unsigned long long IA = 8*m*ceil(atan(1.)/8.) + 5;
IY = IY * IA + IC; //first random number
//inverse error function coefficients
double c[nc] =
{1.,1.,1.16667,1.41111,1.73373,2.14858,
2.67717,3.34815,4.9149,5.27542,6.639,
8.3655,10.5518,13.3208,16.8285,21.2729,
26.9053,34.0446,43.0957,54.5727};
k1 = sqrt(2)*sigma;
for(i=0;i<n;i++){
//probability >0 < 1
p = (double)(IY%(m+1))/(double)(m+1);
IY = IY * IA + IC; //next random number
// inverse error function
p = 2.*p-1.;
q = sqrt(_PI)*p/2.;
einv = 0.;
for(k=0;k<nc;k++) einv += ((c[k])/(2*k+1))*pow(q,2*k+1);
//normal distribution quantile
y[i] = u + einv*k1;
}
}//norm_vars()
int save(int n,double x[]){
FILE *ptr = fopen("data.txt","w");
if(!ptr)return(1);
fprintf(ptr,"%d\n",n);
for(int i = 0;i <n;i++)
fprintf(ptr,"%.8f\n",x[i]);
fclose(ptr);
cout<<"SAVE DONE"<<endl;
return(0);
}//save()
void ave(int n,double y[],double *av,double *std){
// mean and standard deviation
int i;
double asum = 0, ssum = 0;
for(i=0;i<n;i++)asum += y[i];
*av = asum/double(n);
for(i=0;i<n;i++)ssum += pow(*av-y[i],2);
*std = sqrt(ssum/double(n));
}//ave()
double min(int n,double y[]){
int i;
double miny = y[0];
for(i=0;i<n;i++)if(y[i]<miny)miny = y[i];
return miny;
}//min()
double max(int n,double y[]){
int i;
double maxy = y[0];
for(i=0;i<n;i++)if(y[i]>maxy)maxy = y[i];
return maxy;
}//max()
int main() {
unsigned seed = (time(NULL));
srand(seed);
int i, n = 1000;
double u = 10.,av,std = 3.,lo,hi;
double y[n];// = new double[n];
//exp_vars(n,y,u,seed);
norm_vars(n,y,u,std,seed);
ave(n,y,&av,&std);
qsort(y, n, sizeof(double), comp_doub_ascend);
//for(i=0;i<n;i++)printf("%10.6f\n",y[i]);
save(n,y);
lo = min(n,y);
hi = max(n,y);
printf("min %10.4f \tmax %10.4f\n",lo,hi);
printf("ave %10.4f \tstd %10.4f\n",av,std);
double z1,z2;
for(i=1;i<10;i++){
z1 = double(i)/10.;
z2 = erf(erfinv(z1));
printf("%.6f %.6f\n",z1,z2);
}
return(0);
}//main()
