#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#define Nmax 100
double TODA(double rx);
double SHENTON(double u);
/*-----------------------------------------------------------------------*/

int main(int argc, char* argv[])
{
    FILE *fin,*fout;
    char s[256];
    int i,ndata;
    double ppoint,prbs,tpoint;
    double xd[Nmax];
    double yd[Nmax];

    /* データ読み込み */
    fin=fopen(argv[1],"r");
    fgets(s,sizeof s,fin);
    fgets(s,sizeof s,fin);ndata=atoi(strtok(s,","));
    i=0;
    while(fgets(s,sizeof s,fin)!=NULL){
        xd[i]=atof(strtok(s,","));
        yd[i]=atof(strtok(NULL,","))+0.5;
        i=i+1;
    }
    fclose(fin);
    ndata=i-1;

    fout=fopen(argv[2],"w");
    fprintf(fout,"％点(入力値),非超過確率(table),非超過確率(Shenton),％点(戸田)\n");
    for(i=0;i<=ndata;i++){
        ppoint=xd[i];
        prbs=SHENTON(ppoint);
        tpoint=TODA(prbs);
        fprintf(fout,"%e,%e,%e,%e\n",xd[i],yd[i],prbs,tpoint);
    }
    fclose(fout);

    return 0;
}
/*-----------------------------------------------------------------------*/

double TODA(double rx)
{
    /* 非超過確率の％点 */
    int i;
    double ry,ay,uay,sum;
    double b[11];

    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=0.0;
        for(i=0;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 1.0-q;
}
/*-----------------------------------------------------------------------*/