/***************************************************/
/* gccSFT_LP3.c (F-test) */
/***************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#define NDmax 1000
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 CALPP_LP3(double a_lp3,double b_lp3,double c_lp3,double xx);
double CALPP_LP3SP(int nd,double datat[NDmax],double xx);
double GAMMAL(double x);
double TODA(double rx);
double SHENTON(double u);
double GDISTPP(double eta,double p);
double LGAMMA1(double xx);
double LGAMMA2(double x);
double LINCGF(double y1,double z1);
double FUP(int n1,int n2,double u);
/*---------------------------------------------------------------------------*/
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,nd,md;
    double datax[NDmax+1],datay[NDmax+1];
    double beta0,xeps,ueps,eps0,eps,px,fv;
    double a_lp3,b_lp3,c_lp3;
    char dat[500]="";
    char sv[50]="";


    FILE *fin,*fout;

    beta0=0.05;
    for(i=0;i<=NDmax;i++){datax[i]=0.0;datay[i]=0.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;
    /*========================*/
    /* Sorting in small order */
    /*========================*/
    datax[0]=-99.9;
    qsort(datax,nd+1,sizeof(datax[0]),COMP_DS);
    eps0=1.0-pow(1.0-beta0,1.0/(double)nd);
    md=nd-1;
    xeps=datax[nd];
    LP3_COEF(md,datax,&a_lp3,&b_lp3,&c_lp3);
    if(b_lp3<1e4){
        px=CALPP_LP3(a_lp3,b_lp3,c_lp3,xeps);
    }else{
        px=CALPP_LP3SP(nd,datay,xeps);
    }
    ueps=TODA(px);
    fv=(double)(md-1)/(double)(md+1)*ueps*ueps;
    eps=0.5*FUP(1,md-1,fv);

    printf("md=%d\n",md);
    printf("a=%g b=%g c=%g\n",a_lp3,b_lp3,c_lp3);
    printf("xeps=%g prob_x=%g T=%.0f\n",xeps,px,1.0/(1.0-px));
    printf("ueps=%g fv=%g 2eps=%g\n",ueps,fv,2.0*eps);
    if(eps>eps0)printf("eps=%g eps0=%g Not critical region\n",eps,eps0);
    if(eps<=eps0)printf("eps=%g eps0=%g Critical region\n",eps,eps0);

    fout=fopen(fnameW,"w");
        fprintf(fout,"md=%d\n",md);
        fprintf(fout,"a=%g b=%g c=%g\n",a_lp3,b_lp3,c_lp3);
        fprintf(fout,"xeps=%g prob_x=%g T=%.0f\n",xeps,px,1.0/(1.0-px));
        fprintf(fout,"ueps=%g fv=%g 2eps=%g\n",ueps,fv,2.0*eps);
        if(eps>eps0)fprintf(fout,"eps=%g eps0=%g Not critical region\n",eps,eps0);
        if(eps<=eps0)fprintf(fout,"eps=%g eps0=%g Critical region\n",eps,eps0);
    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 CALPP_LP3(double a_lp3,double b_lp3,double c_lp3,double xx)
{
    double wp,etad,qd,et,g1,px;

    wp=(log(xx)-c_lp3)/a_lp3;
    etad=b_lp3;
    qd=log(wp);
    et=log(etad);
    g1=LINCGF(qd,et);
    px=exp(g1-LGAMMA1(etad));

    if(a_lp3<0.0)px=1.0-px;

    return px;
}
/*---------------------------------------------------------------------------*/
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 CALPP_LP3SP(int nd,double datat[NDmax],double xx)
{
    int i;
    double wmy,wsy,wry,ws,wc,w1,w2;
    double kp,ry,z;

    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);

    kp=(log(xx)-wmy)/wsy;
    ry=wry;
    z=6.0/ry*(pow(ry/2.0*(kp+2.0/ry),1.0/3.0)+ry*ry/36.0-1.0);

    return 1.0-SHENTON(z);
}
/*---------------------------------------------------------------------------*/
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 SHENTON(double u)
{
    /* 標準正規分布％点に対する上側確率 */
    int nt=30;
    int ii;
    double ue,af,qh,ag,ab,q;

    ue=0.0;
    if(u>=0)ue=u;
    if(u<0)ue=-u;
    af=ue*ue;
    qh=0.0;
    ag=exp(-0.5*af)*0.3989422804014327;
    ab=-1.0;
    for(ii=nt;ii>=1;ii--){
        qh=(double)ii*af/(2.0*(double)ii+1.0+ab*qh);
        ab=-ab;
    }
    qh=0.5-ag*ue/(1.0-qh);
    if(u<0)qh=1.0-qh;
    q=qh;
    return q;
}
/*-----------------------------------------------------------------------*/
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);
}
/*-----------------------------------------------------------------------*/
double FUP(int n1,int n2,double u)
{
    /* F分布の上側確率 2*epsilon */
    int i,j,mdf,ndf;
    double p,w1,w2,w3;
    double pi=M_PI;

    w1=(double)n2/((double)n2+(double)n1*u);
    w2=sqrt(1.0-w1);
    mdf=2*(int)(n1/2)-n1+2;
    ndf=2*(int)(n2/2)-n2+2;
    if(1<mdf&&1<ndf){
        p=w1;
        w3=p*(1.0-w1);
    }
    if(1<mdf&&ndf<=1){
        p=sqrt(w1);
        w3=p*(1.0-w1)*0.5;
    }
    if(mdf<=1&&1<ndf){
        p=1.0-w2;
        w3=w1*w2*0.5;
    }
    if(mdf<=1&&ndf<=1){
        p=1.0-2.0/pi*atan(w2/sqrt(w1));
        w3=w2*sqrt(w1)/pi;
    }
    for(i=ndf;i<=n2;i=i+2){
        j=i;
        if(n2<=i)break;
        p=p-2.0/(double)i*w3;
        w3=w3*w1*(double)(mdf+i)/(double)i;
    }
    for(i=mdf;i<=n1;i=i+2){
        if(n1<=i)break;
        p=p+2.0/(double)i*w3;
        w3=w3*(1.0-w1)*(double)(i+j)/(double)i;
    }
    return p;
}
/*---------------------------------------------------------------------------*/