/*
 calculation of PI.
 
 Example compilation:
 bcc32 gu015c.c fftsg_h.c
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

void cal(int);
void fmul1604(short int *,short int *,short int *,int,double *,double *);
void fmul1604_herf(short int *,short int *,short int *,int,double *,double *);
int  add1604(short int *,short int *,short int *,int);
int  sub1604(short int *,short int *,short int *,int);
void xmul1604(short int *,double,int);
int  xmul1604_add(short int *,double,short int *,int);
void div4(short int *,int);
void div2h(short int *,short int *,int);
void smul2x2(short int *,short int *,short int *);
void kuriage_h(short int *,int,int);
void kuriage(short int *,int);
void kurisage(short int *,int);
void result1604_file(short int *,int);
/* fft */
void fft_dataset(short int *,int,double *,int);
void fft_fuku_xx(double *,int);
void fft_fuku   (double *,double *,int);
void fft_dataput(short int *,int,double *,int);
int  fft_dataput_add(short int *,int,double *,int);
/* rrsqrt */
void rrsqrt(short int *,short int *,short int *,short int *,int,double *,double *);
void gksqrt_precalc(short int *,short int *,short int *,int,double *,double *);
void rrsqrt_maincalc(short int *,short int *,short int *,short int *,int,double *,double *);
/* gyaku */
void gyaku(short int *,short int *,short int *,int,double *,double *);
void gyaku_precalc(short int *,short int *,short int *,int,double *,double *);
void gyaku_maincalc(short int *,short int *,short int *,int,double *,double *);
/* sqrt & gyaku */
void sub9999(short int *,int);
/* fftsg_h.c */
void rdft(int, int, double *);

void main(void)
{
    int bx,snum;
    
    printf("length of pi =? ");
    scanf("%d",&snum);
    for(bx=1;bx<snum;bx*=2);
    cal(bx);
}

/*
    a=1.0;
    b=1.0/2.0;
    t=1.0/2.0;
    x=2;
    for(i=0;i<n;i++){
        y=(a+b)/4.0;
        b=sqrt(a*b);
        a=y+b/2.0;
        t=t-x*(a-b);
        x*=2;
    }
    pi=(a+b)/t;
*/

void cal(int bx)
{
    int i,x,kt,log2_bx,ret;
    short int *a,*b,*t,*y,*u,*z1;
    double *acp,*bcp;
    time_t stime,rtime;
    
    stime=clock();
    
    kt=bx/4;
    for(log2_bx = 1; (1 << log2_bx) < bx; log2_bx++);
    
    a =(short int *)malloc((kt+4   )*sizeof(short int));
    b =(short int *)malloc((kt+4   )*sizeof(short int));
    t =(short int *)malloc((kt+4   )*sizeof(short int));
    u =(short int *)malloc((kt+4   )*sizeof(short int));
    y =(short int *)malloc((kt+4   )*sizeof(short int));
    z1=(short int *)malloc((kt+4+10)*sizeof(short int));
    acp=(double   *)malloc((kt     )*sizeof(double));
    bcp=(double   *)malloc((kt     )*sizeof(double));
    if(bcp == NULL){
        printf("Memory Allocation Failure!\n");
        exit(1);
    }
    printf("bx = %d  ",bx);
    printf("use memory = %d\n",(kt+4)*sizeof(short int)*6+kt*sizeof(double)*2);
    
    rtime=clock();
    x=2;
    a[0]=10000;for(i=1;i<kt+4;i++) a[i]=0;
    b[0]= 5000;for(i=1;i<kt+4;i++) b[i]=0;
    t[0]= 5000;for(i=1;i<kt+4;i++) t[i]=0;
    printf("-- %dmsec\n",clock()-rtime);
    
    for(i=1;i<log2_bx;i++){
        rtime=clock();
        ret=add1604(a,b,y,kt); if(ret==1) y[0]+=10000;
        div4(y,kt);
        fmul1604_herf(a,b,u,kt,acp,bcp);
        rrsqrt(u,b,a,z1,kt,acp,bcp);
        div2h(b,a,kt);
        ret=add1604(a,y,a,kt); if(ret==1) a[0]+=10000;
        sub1604(a,b,y,kt);
        xmul1604(y,(double)x,kt);
        sub1604(t,y,t,kt);
        x*=2;
        printf("%2d %dmsec\n",i,clock()-rtime);
    }
    
    rtime=clock();
    ret=add1604(a,b,a,kt); if(ret==1) a[0]+=10000;
    gyaku(t,b,z1,kt,acp,bcp);
    fmul1604_herf(a,b,y,kt,acp,bcp);
    result1604_file(y+1,kt-1);
    printf("-- %dmsec\n",clock()-rtime);
    
    free(a);
    free(b);
    free(t);
    free(u);
    free(y);
    free(z1);
    free(acp);
    free(bcp);
    
    printf("Total time : %d msec.\n",clock()-stime);
}

void gksqrt_precalc(short int *a,short int *x,short int *t,int kt,double *acp,double *bcp)
{
    /* x=x+x*(1-a*x*x)/2 */
    
    int i,shr,ret;
    
    for(i=0;i<kt+6;i++) t[i]=0;
    shr=kt*2;
    
    /* step1 a*x */
    fft_dataset(a,kt,acp,shr);
    rdft(shr, 1, acp);
    fft_dataset(x,kt,bcp,shr);
    rdft(shr, 1, bcp);
    fft_fuku(acp,bcp,shr);
    rdft(shr,-1, acp);
    fft_dataput(t,kt+2,acp,shr);
    
    /* step2 a*x*x */
    t++;
    fft_dataset(t,kt,acp,shr);
    rdft(shr, 1, acp);
    fft_fuku(acp,bcp,shr);
    rdft(shr,-1, acp);
    fft_dataput(t,kt+2,acp,shr);
    
    if(t[0]==1){
        t++;
        /* step3 x*(1-a*x*x)/2 */
        fft_dataset(t,kt,acp,shr);
        rdft(shr, 1, acp);
        fft_fuku(acp,bcp,shr);
        rdft(shr,-1, acp);
        fft_dataput(t,kt+2,acp,shr);
        div2h(t,t,kt+2);
        
        /* step4 x+x*(1-a*x*x)/2 */
        ret=sub1604(x,t,x,kt+2);
        if(ret==-1){printf("error gsqpre s\n");exit(1);}
    }
    else{
        t++;
        /* step3 x*(1-a*x*x)/2 */
        sub9999(t,kt+2);
        fft_dataset(t,kt,acp,shr);
        rdft(shr, 1, acp);
        fft_fuku(acp,bcp,shr);
        rdft(shr,-1, acp);
        fft_dataput(t,kt+2,acp,shr);
        div2h(t,t,kt+2);
        
        /* step4 x+x*(1-a*x*x)/2 */
        ret=add1604(x,t,x,kt+2);
        if(ret==1){printf("error gsqpre a\n");exit(1);}
    }
}

void fmul1604(short int *a,short int *b,short int *out,int kt,double *acp,double *bcp)
{
    int shr;
    
    shr=kt*2;
    
    fft_dataset(a,kt,acp,shr);
    rdft(shr, 1, acp);
    fft_dataset(b,kt,bcp,shr);
    rdft(shr, 1, bcp);
    fft_fuku(acp,bcp,shr);
    rdft(shr,-1, acp);
    fft_dataput(out,kt,acp,shr);
}

void fmul1604_herf(short int *a,short int *b,short int *out,int kt,double *acp,double *bcp)
{
    int i,shr,ret;
    
    if( a==out || b==out ){printf("error m1604h\n");exit(1);}
    
    shr=kt;
    for(i=0;i<kt+2;i++) out[i]=0;
    
    smul2x2(a+kt/2,b+kt/2,out+kt);
    
    fft_dataset(a+kt/2,kt/2,acp,shr);
    rdft(shr, 1, acp);
    fft_dataset(b     ,kt/2,bcp,shr);
    rdft(shr, 1, bcp);
    fft_fuku(acp,bcp,shr);
    rdft(shr,-1, acp);
    ret=fft_dataput_add(out+kt/2,kt/2+2,acp,shr);
    if(ret==1) kuriage(out,kt/2);
    
    fft_dataset(a     ,kt/2,acp,shr);
    rdft(shr, 1, acp);
    fft_fuku(bcp,acp,shr);
    rdft(shr,-1, bcp);
    fft_dataput_add(out,kt,bcp,shr);
    
    fft_dataset(b+kt/2,kt/2,bcp,shr);
    rdft(shr, 1, bcp);
    fft_fuku(acp,bcp,shr);
    rdft(shr,-1, acp);
    ret=fft_dataput_add(out+kt/2,kt/2+2,acp,shr);
    if(ret==1) kuriage(out,kt/2);
    
    out[kt]=0;out[kt+1]=0;
}

void fft_dataput(short int *c,int out_len,double *acp,int shr)
{
    int i,deep;
    double scale,rgs;
    __int64 cy;
    
    scale=2.0/shr;
    
    deep=shr-out_len;
    if(deep>3) deep=3;
    
    cy=0;
    for(i=out_len-1+deep;i>out_len-1;i--){
        rgs=acp[i-1]*scale+0.5+cy;
        cy=rgs*0.0001;
    }
    for(i=out_len-1;i>0;i--){
        rgs=acp[i-1]*scale+0.5+cy;
        cy=rgs*0.0001;
        c[i]=rgs-(double)cy*10000.0;
    }
    c[0]=cy;
}

int fft_dataput_add(short int *c,int out_len,double *acp,int shr)
{
    int i,deep;
    double scale,rgs;
    __int64 cy;
    
    scale=2.0/shr;
    
    deep=shr-out_len;
    if(deep>3) deep=3;
    
    cy=0;
    for(i=out_len-1+deep;i>out_len-1;i--){
        rgs=acp[i-1]*scale+0.5+cy;
        cy=rgs*0.0001;
    }
    for(i=out_len-1;i>0;i--){
        rgs=acp[i-1]*scale+0.5+cy+c[i];
        cy=rgs*0.0001;
        c[i]=rgs-(double)cy*10000.0;
    }
    c[0]+=cy;
    
    if(c[0]>=10000){
        c[0]-=10000;
        return(1);
    }
    else{
        return(0);
    }
}

void fft_fuku(double *acp,double *bcp,int shr)
{
    int i;
    double ra,rb,rc,rd;
    
    acp[0]=acp[0]*bcp[0];
    acp[1]=acp[1]*bcp[1];
    for(i=2;i<shr;i+=2){
        ra=acp[i  ];
        rb=acp[i+1];
        rc=bcp[i  ];
        rd=bcp[i+1];
        acp[i  ]=ra*rc-rb*rd;
        acp[i+1]=rb*rc+ra*rd;
    }
}

void fft_dataset(short int *a,int kt,double *acp,int shr)
{
    int i;
    
    for(i=0;i<kt;i++) acp[i]=a[i];
    for(i=kt;i<shr;i++) acp[i]=0.0;
}

void fft_fuku_xx(double *acp,int shr)
{
    int i;
    double ra,rb;
    
    acp[0]=acp[0]*acp[0];
    acp[1]=acp[1]*acp[1];
    for(i=2;i<shr;i+=2){
        ra=acp[i  ];
        rb=acp[i+1];
        acp[i  ]=ra*ra-rb*rb;
        acp[i+1]=ra*rb*2;
    }
}

void div4(short int *a,int kt)
{
    int i,d,rem;
    
    rem=0;
    for(i=0;i<kt;i++){
        d=a[i]+rem;
        a[i]=d>>2;
        rem=(d&3)*10000;
    }
}

void div2h(short int *a,short int *b,int kt)
{
    int i,d,rem;
    
    rem=0;
    for(i=0;i<kt;i++){
        d=a[i]+rem;
        b[i]=d>>1;
        rem=(d&1)*10000;
    }
}

void result1604_file(short int *c,int out_len)
{
    FILE *fp;
    int i,fe,count,p01;
    char file_name[13];
    
    sprintf(file_name,"pi.txt");
    
    if((fp=fopen(file_name,"w"))==NULL){
        printf("file open error!\n");
        exit(1);
    }
    
    fprintf(fp,"PI=3.\n\n");
    
    fe=out_len/25*25;
    count=0;
    for(i=0;i<fe;i+=25){
        
        fprintf(fp,"%04d%04d%02d ",c[i  ],c[i+1],p01=c[i+2]/100);
        fprintf(fp,"%02d%04d%04d ",c[i+2]-p01*100,c[i+3],c[i+4]);
        fprintf(fp,"%04d%04d%02d ",c[i+5],c[i+6],p01=c[i+7]/100);
        fprintf(fp,"%02d%04d%04d ",c[i+7]-p01*100,c[i+8],c[i+9]);
        fprintf(fp,"%04d%04d%02d\n",c[i+10],c[i+11],p01=c[i+12]/100);
        
        fprintf(fp,"%02d%04d%04d ",c[i+12]-p01*100,c[i+13],c[i+14]);
        fprintf(fp,"%04d%04d%02d ",c[i+15],c[i+16],p01=c[i+17]/100);
        fprintf(fp,"%02d%04d%04d ",c[i+17]-p01*100,c[i+18],c[i+19]);
        fprintf(fp,"%04d%04d%02d ",c[i+20],c[i+21],p01=c[i+22]/100);
        fprintf(fp,"%02d%04d%04d\n",c[i+22]-p01*100,c[i+23],c[i+24]);
        
        count+=100;
        if(count==1000){
            count=0;
            fprintf(fp,"\n");
        }
    }
    
    fe=out_len;
    count=0;
    for(;i<fe;i++){
        fprintf(fp,"%02d",p01=c[i]/100);
        count+=2;
        if(count%50==0) fprintf(fp,"\n"); else if(count%10==0) fprintf(fp," ");
        
        fprintf(fp,"%02d",c[i]-p01*100);
        count+=2;
        if(count%50==0) fprintf(fp,"\n"); else if(count%10==0) fprintf(fp," ");
    }
    fprintf(fp,"\n");
    
    fclose(fp);
}

void gyaku(short int *a,short int *x,short int *t,int kt,double *acp,double *bcp)
{
    int i,bp;
    double gy,sm;
    
    a[kt]=0;a[kt+1]=0;
    
    for(i=0;i<kt+4;i++) x[i]=0;
    
    gy=0.0;
    sm=0.0001;
    for(i=0;i<5;i++){
        gy+=a[i]*sm;
        sm*=0.0001;
    }
    gy=1.0/gy;
    
    while((int)gy==0) gy*=10000.0;
    for(i=0;i<5;i++){
        x[i]=gy;
        gy-=x[i];
        gy*=10000.0;
    }
    
    gyaku_precalc(a,x,t,8,acp,bcp);
    for(i=6;i<20;i++) x[i]=0;
    
    bp=8;
    while(bp<=kt){
        /* x=2x-Axx */
        gyaku_maincalc(a,x,t,bp,acp,bcp);
        bp*=2;
    }
    x[kt]=0;x[kt+1]=0;
    t[kt]=0;t[kt+1]=0;
}

void gyaku_precalc(short int *a,short int *x,short int *t,int kt,double *acp,double *bcp)
{
    /* x=x+x*(1-ax) */
    
    int i,shr,ret;
    
    for(i=0;i<kt+6;i++) t[i]=0;
    shr=kt*2;
    
    /* step1 a*x */
    fft_dataset(a,kt,acp,shr);
    rdft(shr, 1, acp);
    fft_dataset(x,kt,bcp,shr);
    rdft(shr, 1, bcp);
    fft_fuku(acp,bcp,shr);
    rdft(shr,-1, acp);
    fft_dataput(t,kt+2,acp,shr);
    
    if(t[0]==1){
        t++;
        /* step2 x*(1-a*x) */
        fft_dataset(t,kt,acp,shr);
        rdft(shr, 1, acp);
        fft_fuku(acp,bcp,shr);
        rdft(shr,-1, acp);
        fft_dataput(t,kt+2,acp,shr);
        /* step3 x+x*(1-a*x) */
        ret=sub1604(x,t,x,kt+2);
        if(ret==-1){printf("error gakpre s\n");exit(1);}
    }
    else{
        t++;
        /* step2 x*(1-a*x) */
        sub9999(t,kt+2);
        fft_dataset(t,kt,acp,shr);
        rdft(shr, 1, acp);
        fft_fuku(acp,bcp,shr);
        rdft(shr,-1, acp);
        fft_dataput(t,kt+2,acp,shr);
        /* step3 x+x*(1-a*x) */
        ret=add1604(x,t,x,kt+2);
        if(ret==1){printf("error gakpre a\n");exit(1);}
    }
}

void gyaku_maincalc(short int *a,short int *x,short int *t,int kt,double *acp,double *bcp)
{
    /* x=x+x*(1-ax) */
    
    int i,shr,ret;
    short int h1[2],h2[2];
    
    for(i=kt;i<kt+6;i++) t[i]=0;
    shr=kt;
    
    /* step1 a*x */
    smul2x2(a+kt,x,t+kt);
    
    fft_dataset(a,kt/2,acp,shr);
    rdft(shr, 1, acp);
    fft_dataset(x,kt/2,bcp,shr);
    rdft(shr, 1, bcp);
    fft_fuku(acp,bcp,shr);
    rdft(shr,-1, acp);
    fft_dataput(t,kt,acp,shr);
    
    fft_dataset(a+kt/2,kt/2,acp,shr);
    rdft(shr, 1, acp);
    fft_fuku(acp,bcp,shr);
    rdft(shr,-1, acp);
    ret=fft_dataput_add(t+kt/2,kt/2+2,acp,shr);
    if(ret==1) kuriage(t,kt/2);
    
    ret=xmul1604_add(a,(double)(x[kt/2]*10000+x[kt/2+1]),t+kt/2+2,kt/2);
    if(ret>0) kuriage_h(t,kt/2+2,ret);
    
    smul2x2(a+kt/2,x+kt/2,h1);
    ret=add1604(t+kt,h1,t+kt,2);
    if(ret==1) kuriage(t,kt);
    
    
    if(t[0]==1){
        
        /* step2 (1-a*x)*x */
        t++;
        smul2x2(t+kt,x,h1);
        
        smul2x2(t+kt/2,x+kt/2,h2);
        
        fft_dataset(t+kt/2,kt/2,acp,shr);
        rdft(shr, 1, acp);
        fft_fuku(acp,bcp,shr);
        rdft(shr,-1, acp);
        fft_dataput(t+kt/2,kt/2+2,acp,shr);
        
        ret=add1604(t+kt,h2,t+kt,2);
        if(ret==1) kuriage(t,kt);
        
        ret=add1604(t+kt,h1,t+kt,2);
        if(ret==1) kuriage(t,kt);
        
        /* step3 x+x*(1-a*x) */
        ret=sub1604(x+kt/2,t+kt/2,x+kt/2,kt/2+2);
        if(ret==-1) kurisage(x,kt/2);
    }
    else{
        /* step2 (1-a*x)*x */
        t++;
        sub9999(t+kt/2,kt/2+2);
        smul2x2(t+kt,x,h1);
        
        smul2x2(t+kt/2,x+kt/2,h2);
        
        fft_dataset(t+kt/2,kt/2,acp,shr);
        rdft(shr, 1, acp);
        fft_fuku(acp,bcp,shr);
        rdft(shr,-1, acp);
        fft_dataput(t+kt/2,kt/2+2,acp,shr);
        
        ret=add1604(t+kt,h2,t+kt,2);
        if(ret==1) kuriage(t,kt);
        
        ret=add1604(t+kt,h1,t+kt,2);
        if(ret==1) kuriage(t,kt);
        
        /* step3 x+x*(1-a*x) */
        ret=add1604(x+kt/2,t+kt/2,x+kt/2,kt/2+2);
        if(ret==1) kuriage(x,kt/2);
    }
}

void kuriage_h(short int *a,int n,int cy)
{
    int i,rgs;
    
    for(i=n-1;i>=0;i--){
        rgs=a[i]+cy;
        cy=rgs*0.0001;
        a[i]=rgs-cy*10000;
        if(cy==0) break;
    }
    if(cy==1){printf("error kuriage_h");exit(1);}
}

void kuriage(short int *a,int n)
{
    int i,cy;
    
    cy=1;
    for(i=n-1;i>=0;i--){
        a[i]+=cy;
        if(a[i]>=10000){
            cy=1;
            a[i]-=10000;
        }
        else{
            cy=0;
            break;
        }
    }
    if(cy==1){printf("error kuriage");exit(1);}
}

void kurisage(short int *a,int n)
{
    int i,cy;
    
    cy=-1;
    for(i=n-1;i>=0;i--){
        a[i]+=cy;
        if(a[i]<0){
            cy=-1;
            a[i]+=10000;
        }
        else{
            cy=0;
            break;
        }
    }
    if(cy==-1){printf("error kurisage");exit(1);}
}

void sub9999(short int *a,int kt)
{
    int i;
    
    for(i=0;i<kt;i++){
        a[i]=9999-a[i];
    }
}

void smul2x2(short int *a,short int *b,short int *c)
{
    int i,j,cy,t[4];
    
    for(i=0;i<4;i++) t[i]=0;
    
    for(i=0;i<2;i++){
        for(j=0;j<2;j++){
            t[i+j+1]+=a[i]*b[j];
        }
    }
    
    cy=0;
    for(i=3;i>=0;i--){
        t[i]+=cy;
        cy=t[i]*0.0001;
        t[i]-=cy*10000;
    }
    
    for(i=0;i<2;i++){
        c[i]=t[i];
    }
}

void rrsqrt(short int *a,short int *y,short int *x,short int *t1,int kt,double *acp,double *bcp)
{
    int i,bp;
    double gy,sm;
    
    a[kt]=0;a[kt+1]=0;
    
    for(i=0;i<kt+4;i++) x[i]=0;
    for(i=0;i<kt+4;i++) y[i]=0;
    
    gy=0.0;
    sm=0.0001;
    for(i=0;i<5;i++){
        gy+=a[i]*sm;
        sm*=0.0001;
    }
    gy=1.0/sqrt(gy);
    
    while((int)gy==0) gy*=10000.0;
    for(i=0;i<5;i++){
        x[i]=gy;
        gy-=x[i];
        gy*=10000.0;
    }
    
    gy=0.0;
    sm=0.0001;
    for(i=0;i<5;i++){
        gy+=a[i]*sm;
        sm*=0.0001;
    }
    gy=sqrt(gy);
    
    while((int)gy==0) gy*=10000.0;
    for(i=0;i<5;i++){
        y[i]=gy;
        gy-=y[i];
        gy*=10000.0;
    }
    
    gksqrt_precalc(a,x,t1,8,acp,bcp);
    gksqrt_precalc(a,x,t1,16,acp,bcp);
    fmul1604(a,x,y,16,acp,bcp);
    for(i=0;i<15;i++) y[i]=y[i+1];
    
    for(i= 6;i<20;i++) x[i]=0;
    for(i=10;i<20;i++) y[i]=0;
    
    bp=8;
    while(bp<kt){
        /* x=x+x(1-yx)   */
        /* y=y+x(a-yy)/2 */
        gyaku_maincalc(y,x,t1,bp,acp,bcp);
        rrsqrt_maincalc(a,x,y,t1,bp*2,acp,bcp);
        bp*=2;
    }
    
    y[kt]=0;y[kt+1]=0;
    x[kt]=0;x[kt+1]=0;
    t1[kt]=0;t1[kt+1]=0;
}

void rrsqrt_maincalc(short int *a,short int *x,short int *y,short int *t1,int kt,double *acp,double *bcp)
{
    /* x=x+x*(1-a*x*x)/2 */
    /* y=y+x*(a-y*y)/2   */
    
    int i,shr,ret;
    short int h1[2],h2[2];
    
    for(i=kt;i<kt+6;i++) t1[i]=0;
    shr=kt;
    
    /* step1 y*y */
    fft_dataset(y,kt/2,acp,shr);
    rdft(shr, 1, acp);
    fft_fuku_xx(acp,shr);
    rdft(shr,-1, acp);
    fft_dataput(t1,kt,acp,shr);
    
    ret=xmul1604_add(y,(y[kt/2]*10000+y[kt/2+1])*2.0,t1+kt/2+2,kt/2);
    if(ret>0) kuriage_h(t1,kt/2+2,ret);
    
    smul2x2(y+kt/2,y+kt/2,h1);
    ret=add1604(t1+kt,h1,t1+kt,2);
    if(ret==1) kuriage(t1,kt);
    
    /* step2 a-y*y */
    ret=sub1604(a,t1,t1,kt+2);
    
    
    if(ret==-1){
        /* step3 x*(a-y*y) */
        t1++;
        sub9999(t1+kt/2,kt/2+2);
        smul2x2(t1+kt,x,h1);
        
        smul2x2(t1+kt/2,x+kt/2,h2);
        
        fft_dataset(t1+kt/2,kt/2,acp,shr);
        rdft(shr, 1, acp);
        fft_dataset(x      ,kt/2,bcp,shr);
        rdft(shr, 1, bcp);
        fft_fuku(acp,bcp,shr);
        rdft(shr,-1, acp);
        fft_dataput(t1+kt/2,kt/2+2,acp,shr);
        
        ret=add1604(t1+kt,h2,t1+kt,2);
        if(ret==1) kuriage(t1,kt);
        
        ret=add1604(t1+kt,h1,t1+kt,2);
        if(ret==1) kuriage(t1,kt);
        
        /* step4 x*(a-y*y)/2 */
        div2h(t1+kt/2,t1+kt/2,kt/2+2);
        
        /* step5 y+x*(a-y*y)/2 */
        ret=sub1604(y+kt/2,t1+kt/2,y+kt/2,kt/2+2);
        if(ret==1) kurisage(y,kt/2);
    }
    else{
        /* step3 x*(a-y*y) */
        t1++;
        smul2x2(t1+kt,x,h1);
        
        smul2x2(t1+kt/2,x+kt/2,h2);
        
        fft_dataset(t1+kt/2,kt/2,acp,shr);
        rdft(shr, 1, acp);
        fft_dataset(x      ,kt/2,bcp,shr);
        rdft(shr, 1, bcp);
        fft_fuku(acp,bcp,shr);
        rdft(shr,-1, acp);
        fft_dataput(t1+kt/2,kt/2+2,acp,shr);
        
        ret=add1604(t1+kt,h2,t1+kt,2);
        if(ret==1) kuriage(t1,kt);
        
        ret=add1604(t1+kt,h1,t1+kt,2);
        if(ret==1) kuriage(t1,kt);
        
        /* step4 x*(a-y*y)/2 */
        div2h(t1+kt/2,t1+kt/2,kt/2+2);
        
        /* step5 y+x*(a-y*y)/2 */
        ret=add1604(y+kt/2,t1+kt/2,y+kt/2,kt/2+2);
        if(ret==-1) kuriage(y,kt/2);
    }
}

int add1604(short int *a,short int *b,short int *c,int kt)
{
    int i,carry;
    
    carry=0;
    for(i=kt-1;i>=0;i--){
        c[i]=a[i]+b[i]+carry;
        if(c[i]>=10000){
            carry=1;
            c[i]-=10000;
        }
        else{
            carry=0;
        }
    }
    return(carry);
}

int sub1604(short int *a,short int *b,short int *c,int kt)
{
    int i,carry;
    
    carry=0;
    for(i=kt-1;i>=0;i--){
        c[i]=a[i]-b[i]+carry;
        if(c[i]<0){
            carry=-1;
            c[i]+=10000;
        }
        else{
            carry=0;
        }
    }
    return(carry);
}

void xmul1604(short int *a,double x,int kt)
{
    int i;
    double d;
    __int64 cy;
    
    cy=0;
    for(i=kt-1;i>=0;i--){
        d=a[i]*x+cy;
        cy=d*0.0001;
        a[i]=d-(double)cy*10000.0;
    }
}

int xmul1604_add(short int *a,double x,short int *c,int kt)
{
    int i;
    double d;
    __int64 cy;
    
    cy=0;
    for(i=kt-1;i>=0;i--){
        d=a[i]*x+cy+c[i];
        cy=d*0.0001;
        c[i]=d-(double)cy*10000.0;
    }
    return(cy);
}