/*★ 円周率(pai)を計算する                    前田　稔   ★*/
/*   pai= 4*(4*a_tan(1/5) - a_tan(1/239))                  */
/*   a_tan(x)= x - x**3/3 + x**5/5 - x**7/7 + x**9/9 - ... */
/*   1989年12月   PC-9801UV(10MHZ) では、１７分間以上      */
#include <stdio.h>
#define     S       1000        //配列のサイズ（Ｓ＊４桁）
#define     C1      3000        //第一項のループ回数
#define     C2      900         //第二項のループ回数

short   pai[S],t[S];
short   a[S],w[S];

void    a_tan(short[],short,short);
void    a_tan2(short[],short,short);
void    add(short[],short[],short[]);
void    sub(short[],short[],short[]);
void    mlt(short[],short[],short);
void    dev(short[],short[],short);
void    disp(char*,short[]);

//★ MAIN PROGRAM
void    main()
{   short   i;

    //※ a_tan(1/5) の計算
    // pai[]= 1/5
    for(i=0; i<S; i++)  pai[i]= 0;
    pai[1]=2000;
    // a_tan(pai[1/5], 25, 3000)
    a_tan(pai,5*5,C1);
    // pai[]= pai[]*4
    mlt(pai,pai,4);

    //※ a_tan(1/239) の計算
    // t[]= 1/239
    for(i=0; i<S; i++)  t[i]= 0;
    t[0]=1;
    dev(t,t,239);
    // a_tan2(t[1/239], 239, 900)
    a_tan2(t,239,C2);

    //※ (a_tan(1/5) - a_tan(1/239)) * 4
    sub(pai,pai,t);
    mlt(pai,pai,4);
    disp("pai=",pai);
}

//★ a_tan(x)= x - x**3/3 + x**5/5 - x**7/7 + x**9/9 - ...
//   a_tan(pai[1/5], 25, 3000)
//   x==1/5 のときは、毎回 25 で割ればよい
void    a_tan(short x[], short k, short n)
{   short   i,j;
    for(i=0; i<S; i++)  {  a[i]= x[i];  x[i]= 0;  }
    for(i=0; i<n; i++)
    {   j= i+i+1;
        dev(w,a,j);
        if (i & 1)  sub(x,x,w);
        else        add(x,x,w);
        dev(a,a,k);
    }
}

//★ a_tan2(x)= x - x**3/3 + x**5/5 - x**7/7 + x**9/9 - ...
//   a_tan2(t[1/239], 239, 900)
//   x==1/239 のときは、毎回 239 で二回割る
void    a_tan2(short x[], short k, short n)
{   short   i,j;
    for(i=0; i<S; i++)  {  a[i]= x[i];  x[i]= 0;  }
    for(i=0; i<n; i++)
    {   j= i+i+1;
        dev(w,a,j);
        if (i & 1)  sub(x,x,w);
        else        add(x,x,w);
        dev(a,a,k);
        dev(a,a,k);
    }
}

//★ x[]= y[] + z[]
void    add(short x[], short y[], short z[])
{   short   i;
    for(i=0; i<S; i++)  x[i]= y[i]+z[i];
    for(i=S-1; i>0; i--)
        if (x[i]>=10000) {  x[i]= x[i]-10000;  x[i-1]++;  }
}

//★ x[]= y[] - z[]
void    sub(short x[], short y[], short z[])
{   short   i;
    for(i=0; i<S; i++)  x[i]= y[i]-z[i];
    for(i=S-1; i>0; i--)
        if (x[i]<0)     {  x[i]= x[i]+10000;  x[i-1]--;  }
}

//★ x[]= y[] * n
void    mlt(short x[], short y[], short n)
{   short   i,v;
    long    w;
    v= 0;
    for(i=S-1; i>=0; i--)
    {   w= y[i];  w*= n;  w+= v;
        x[i]= w%10000;  v= w/10000;
    }
}

//★ x[]= y[] / n
void    dev(short x[], short y[], short n)
{   short   i;
    long    w;
    w= 0;
    for(i=0; i<S; i++)
    {   w*= 10000;  w+= y[i];
        x[i]= w/n;  w%= n;
    }
}

//★ 配列の表示
void    disp(char s[], short x[])
{   short   i,j,k;
    printf("\n%s\n    1: %d.",s,x[0]);
    for(k=i=1; i<S; k++)
    {   if (k>1)    printf("   ");
        for(j=0; i<S && j<10; j++,i++)  printf("%04d ",x[i]);
        if (i<S)    printf("\n%5d:",k*40+1);
    }
}