/***************************************************/
/* gccTES1_LP3.c */
/***************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#define NDmax 1000
#define NED 1000
#define NYR 14
int COMP_IS(const void *a,const void *b);
int COMP_DS(const void *a,const void *b);
void IRND(int nd,int idata[NDmax]);
void main_part(char fnameR[50],char fnameW[50]);
void LP3_COEF(int nd,double datat[NDmax],double *a_lp3,double *b_lp3,double *c_lp3);
double XP_LP3(double p,double a_lp3,double b_lp3,double c_lp3);
double XP_LP3_SP(int nd,double datat[NDmax],double p);
double GAMMAL(double x);
double TODA(double rx);
double GDISTPP(double eta,double p);
double LGAMMA1(double xx);
double LGAMMA2(double x);
double LINCGF(double y1,double z1);
/*---------------------------------------------------------------------------*/
int main(int argc,char *argv[])
{
    char fnameR[50]; /*Input file name*/
    char fnameW[50]; /*Output file name*/

    strcpy(fnameR,argv[1]);
    strcpy(fnameW,argv[2]);

    main_part(fnameR,fnameW);

    return 0;
}
/*---------------------------------------------------------------------------*/
void main_part(char fnameR[50],char fnameW[50])
{
    int i,j,iii,nd;
    int idata[NDmax+1];
    double datax[NDmax+1],datay[NDmax+1];
    double p,y1,y2;
    double a_lp3,b_lp3,c_lp3;
    int nnn=NYR;
    double Ty[]={0,2,5,10,20,30,50,100,200,300,500,1000,2000,3000,5000};
    char dat[500]="";
    char sv[50]="";

    FILE *fin,*fout,*fhis;

    srand((unsigned int)time(NULL));
    for(i=0;i<=NDmax;i++){datax[i]=0.0;datay[i]=0.0;idata[i]=0;}
    /*============*/
    /* Data input */
    /*============*/
    fin=fopen(fnameR,"r");
    fgets(dat,sizeof dat,fin);
    i=0;
    while(fgets(dat,sizeof dat,fin)!=NULL){i=i+1;datax[i]=atof(dat);}
    fclose(fin);
    nd=i;
    for(i=0;i<=nd;i++){idata[i]=i;}
    /*========================*/
    /* Sorting in small order */
    /*========================*/
    datax[0]=-99.9;
    qsort(datax,nd+1,sizeof(datax[0]),COMP_DS);
   /*========================*/
    /* Calc. of BTS           */
    /*========================*/
    fout=fopen(fnameW,"w");
    fprintf(fout,"i,iii,p,a_lp3,b_lp3,c_lp3,y1,y2,b_lp3,y2/y1\n");
    for(i=1;i<=NED;i++){
        IRND(nd,idata);
        for(j=1;j<=nd;j++){datay[j]=datax[idata[j]];}
            LP3_COEF(nd,datay,&a_lp3,&b_lp3,&c_lp3);
            for(iii=1;iii<=nnn;iii++){
                p=1.0-1.0/Ty[iii];
                y1=XP_LP3(p,a_lp3,b_lp3,c_lp3);
                y2=XP_LP3_SP(nd,datay,p);
                printf("%d,%d,%g,%g,%g,%g,%g,%g,%g,%g\n",i,iii,p,a_lp3,b_lp3,c_lp3,y1,y2,b_lp3,y2/y1);
                fprintf(fout,"%d,%d,%g,%g,%g,%g,%g,%g,%g,%g\n",i,iii,p,a_lp3,b_lp3,c_lp3,y1,y2,b_lp3,y2/y1);
        }
    }
    fclose(fout);

}
/*---------------------------------------------------------------------------*/
int COMP_IS(const void *a,const void *b)
{
    /*小さい順*/
    return *(int*)a-*(int*)b;
}
/*---------------------------------------------------------------------------*/
int COMP_DS(const void *a,const void *b)
{
    /*比較関数は戻り値が整数型なので*/
    /*整数値を返すようにする*/
    /*小さい順*/
    if(*(double*)a> *(double*)b)return  1;
    if(*(double*)a==*(double*)b)return  0;
    if(*(double*)a< *(double*)b)return -1;
}
/*---------------------------------------------------------------------------*/
void IRND(int nd,int idata[NDmax])
{
    int i;
    double r,a,b;

    idata[0]=0;
    a=1.0;
    b=(double)nd;
    for(i=1;i<=nd;i++){
        r=a+(b-a)*(double)rand()/(double)RAND_MAX;
        idata[i]=(int)(r+0.5);
    }
    qsort(idata,nd+1,sizeof(idata[0]),COMP_IS);
}
/*---------------------------------------------------------------------------*/
void LP3_COEF(int nd,double datat[NDmax],double *a_lp3,double *b_lp3,double *c_lp3)
{
    /* Log Pearson III distribution (Gamma distribution) */
    int i;
    double wmy,wsy,wry,ws,wc,w1,w2;
    
    wmy=0.0;
    for(i=1;i<=nd;i++){wmy=wmy+log(datat[i]);}
    wmy=wmy/(double)nd;
    ws=0.0;
    for(i=1;i<=nd;i++){ws=ws+(log(datat[i])-wmy)*(log(datat[i])-wmy);}
    ws=sqrt(ws/(double)nd);
    wsy=sqrt((double)nd/(double)(nd-1))*ws;
    wc=0.0;
    for(i=1;i<=nd;i++){wc=wc+pow((log(datat[i])-wmy)/ws,3.0);}
    wc=wc/(double)nd;
    w1=1.0+6.51/(double)nd+20.2/(double)nd/(double)nd;
    w2=1.48/(double)nd+6.77/(double)nd/(double)nd;
    wry=wc*(w1+w2*wc*wc);

    *b_lp3=4.0/(wry*wry);
    *a_lp3=wsy/sqrt(*b_lp3);
    if(wry<0.0)*a_lp3=-(*a_lp3);
    *c_lp3=wmy-(*a_lp3)*(*b_lp3);

if(*b_lp3<0.0)printf("%g\n",wry);

}
/*---------------------------------------------------------------------------*/
double XP_LP3(double p,double a_lp3,double b_lp3,double c_lp3)
{
    double q;

    q=p;
    if(a_lp3<0.0)q=1.0-p;
    return exp(c_lp3+a_lp3*GDISTPP(b_lp3,q));
}
/*---------------------------------------------------------------------------*/
double XP_LP3_SP(int nd,double datat[NDmax],double p)
{
    /* Log Pearson III distribution (Gamma distribution) */
    int i;
    double wmy,wsy,wry,ws,wc,w1,w2;
    double q,kp,ry;

    wmy=0.0;
    for(i=1;i<=nd;i++){wmy=wmy+log(datat[i]);}
    wmy=wmy/(double)nd;
    ws=0.0;
    for(i=1;i<=nd;i++){ws=ws+(log(datat[i])-wmy)*(log(datat[i])-wmy);}
    ws=sqrt(ws/(double)nd);
    wsy=sqrt((double)nd/(double)(nd-1))*ws;
    wc=0.0;
    for(i=1;i<=nd;i++){wc=wc+pow((log(datat[i])-wmy)/ws,3.0);}
    wc=wc/(double)nd;
    w1=1.0+6.51/(double)nd+20.2/(double)nd/(double)nd;
    w2=1.48/(double)nd+6.77/(double)nd/(double)nd;
    wry=wc*(w1+w2*wc*wc);

    q=p;
/*    if(wry<0.0)q=1.0-p;*/
/*    ry=fabs(wry);*/
    ry=wry;
    kp=2.0/ry*pow((1.0+ry*TODA(q)/6.0-ry*ry/36.0),3.0)-2.0/ry;

    return exp(wmy+wsy*kp);
}
/*---------------------------------------------------------------------------*/
double GAMMAL(double x)
{
    return exp(LGAMMA1(x));
}
/*---------------------------------------------------------------------------*/
double TODA(double rx)
{
    /* %-point of Non-exceeding probability */
    int i;
    double ry,ay,uay;
    double b[11],sum;

    if(rx==0.5){
        uay=0.0;
    }else{
        if(rx<0.5){
            ry=rx;
        }else{
            ry=1-rx;
        }
        ay=-log(4.0*ry*(1.0-ry));
        b[0] = 1.570796288;
        b[1] = 0.03706987906;
        b[2] =-0.0008364353589;
        b[3] =-0.0002250947176;
        b[4] = 0.000006841218299;
        b[5] = 0.000005824238515;
        b[6] =-0.00000104527497;
        b[7] = 0.00000008360937017;
        b[8] =-0.000000003231081277;
        b[9] = 0.00000000003657763036;
        b[10]= 0.0000000000006936233982;
        sum=b[0];
        for(i=1;i<=10;i++){sum=sum+b[i]*pow(ay,(double)i);}
        uay=sqrt(ay*sum);
        if(rx<0.5)uay=-uay;
    }
    return uay;
}
/*-----------------------------------------------------------------------*/
double GDISTPP(double eta,double p)
{
    /* eta=b_lp3,p=probability */
    int k;
    double etad,bd,et,g1,x,dk,xx,a,b,expa,q,qd,er;

    etad=eta;
    bd=p;
    et=log(etad);
    g1=log(bd)+LGAMMA1(etad);
    x=(et+g1)*1.0/etad;
    k=0;
    do{
        k=k+1;
        dk=(double)k;
        xx=x+log(dk);
    }while(LINCGF(xx,et)<g1);
    a=0.0;
    if(2<=k)a=x+log(dk-1.0);
    b=x+log(dk);
    expa=exp(a);
    if(k<=1)expa=0.0;
    k=0;
    do{
        k=k+1;
        q=0.5*(expa+exp(b));
        qd=log(q);
        er=LINCGF(qd,et)-g1;
        if(er<0.0){a=qd;expa=exp(a);}
        if(fabs(er)<1e-10)break;
        if(0.0<er)b=qd;
    }while(k<100000);

    return q;
}
/*-----------------------------------------------------------------------*/
double LGAMMA1(double xx)
{
    /* ln(GAMMA(xx))  */
    int i,m;
    double pi=M_PI;
    double x,gs,y,z,y1;

    x=xx;
    if(30.0<=x){
        gs=(x-0.5)*log(x)-x+0.5*log(2.0*pi)+1.0/(12.0*x)-1.0/(360.0*pow(x,3.0))+1.0/(1260.0*pow(x,5.0))-1.0/(1680.0*pow(x,7.0));
        if(10000.0<x)return gs;
        gs=gs+5.0/(5940.0*pow(x,9.0))-691.0/(360360.0*pow(x,11.0));
        return gs;
    }else{
        m=(int)x;
        y=x-(double)m;
        gs=LGAMMA2(y);
        if(x<=2.0){
            if(m<1)gs=gs-log(y);
            return gs;
        }else{
            z=0;
            for(i=1;i<=m-1;i++){y1=y+(double)i;z=z+log(y1);}
            gs=gs+z;
            return gs;
        }
    }
}
/*-----------------------------------------------------------------------*/
double LGAMMA2(double x)
{
    int i;
    double c[9],g;

    c[1]=-0.577191652;
    c[2]= 0.988205891;
    c[3]=-0.897056937;
    c[4]= 0.918206857;
    c[5]=-0.756704078;
    c[6]= 0.482199394;
    c[7]=-0.193527818;
    c[8]= 0.035868343;
    g=1.0;
    for(i=1;i<=8;i++){g=g+c[i]*pow(x,(double)i);}
    g=log(g);
    return g;
}
/*-----------------------------------------------------------------------*/
double LINCGF(double y1,double z1)
{
    int j;
    double z,y,f,d1,d2;

    z=exp(z1);
    y=exp(y1);
    f=1.0/z;
    d1=1.0/z;
    j=0;
    do{
        j=j+1;
        d2=d1*y/(z+(double)j);
        f=f+d2;
        if(y/(z+(double)j+1.0-y)*d2<f*1e-10&&y/(z+(double)j+1.0)<1.0)break;
        d1=d2;
    }while(j<100000);

    return z*y1-y+log(f);
}
/*-----------------------------------------------------------------------*/