//  Atan公式によるπ計算プログラム(N**2版、Cプログラム）
//　　　　　　　　　　　　　　　　　　　　　　2002/9/26  後　保範
//=================================================================C
//  Pi Value Computation Program by Atan Normal Method             C
//    Pi = M1*Atan(1/X1) + M2*Atan(1/X2) + M3*Atan(1/X3)           C
//    Atan(1/X) = 1/X - 1/(3*X**3) + 1/(5*X**5) - 1/(7*X**7) +..   C
//-----------------------------------------------------------------C
//    Written by Yasunori Ushiro,  2002/09/09                      C
//        (Hitachi Ltd. and Waseda University)                     C
//=================================================================C
#include <stdio.h> 
#include <math.h>
#include <time.h>
// Global Define
#define ND   1000010
#define L10V 5
double X, P10, D10;
int M, N2, L10;
int IV[ND+1], IXD[ND+1];
int INO[5] = {2,3,3,4,4};    // π計算公式のAtanの個数 
// π計算用各公式の定数　(M,X)の組、最大4個(Atanの個数）
int IDATA[5][4][2] = { 16, 5, -4,239, 0,0, 0,0,
   48,18,  32, 57, -20,239,  0,0,      32,10, -4,239, -16,515,  0,0,
  176,57, 28,239, -48,682, 96,12943,  48,49, 128,57, -20,239, 48,110443};
//=================================================================C
//   IV = M*Atan(1/X) Computation by Normal Method                 C
//      IXD = 1/X**K                                               C
//-----------------------------------------------------------------C
//    Written by Yasunori Ushiro,  2002/09/08                      C
//        (Hitachi Ltd. and Waseda University)                     C
//=================================================================C
  void ATANV()     // M*Atan(1/X)の多数桁の値の計算
  { int    i, k, KOSU, LST, ISG;
    double X2, DX2, W, W1, W2, VK, DK, WK, WK1, WK2;
//  Initial Constant
      X2   = X*X ;
      DX2  = 1.0/X2 ;
      KOSU = N2*L10/(log(X2)/log(10.0));   // 多数桁に必要なAtanの項数
//  Initial Value Set
      for(i=2; i<=N2; i++)
       { IXD[i] = 0 ; }
      IXD[1] = M ;
//  Initial Computation, IV=IXD=M/X
      W2 = 0.0;
      for(i=1; i<=N2; i++)
       { W  = IXD[i] + W2*P10 ;
         W1 = floor(W/X) ;
         W2 = W - W1*X ;
         IXD[i] = W1 ;
         IV[i]  = IV[i] + IXD[i] ; }    // IV = M/X + 前のAtanの計算値
//  Main Computation
      LST =  1 ;
      ISG = -1 ;
      for (k=1; k<=KOSU; k++)           //  計算項数だけ反復計算
       { W2  = 0.0 ;
         WK2 = 0.0 ;
         VK  = 2*k + 1 ;
         DK  = 1.0/VK ;
         for(i=LST; i<=N2; i++)          // 多数桁は1要素(int)
          { W      = IXD[i] + W2*P10 ;   // に10進5桁ごとに分けて記憶
            W1     = floor(W*DX2) ;      // floor: 切り捨てて整数値にする
            W2     = W - W1*X2 ;
            IXD[i] = W1 ;                // W1 = 多数桁のM/X**(2k+1)の値
            WK     = W1 + WK2*P10 ;      // W1 = 旧W1/X**2 = 旧W1*DX2
            WK1    = floor(WK*DK) ;
            WK2    = WK - WK1*VK ;       // WK2 =　多数桁のW1/(2k+1)
            IV[i]  = IV[i] + ISG*WK1 ; }
         if(IXD[LST] == 0)
	      { LST = LST + 1;               // 無駄な計算を省くための処理
            if(IXD[LST] == 0)  LST = LST + 1;  }
         ISG = -ISG ;   }
    }
//=================================================================C
  void NORMAL()       // 多数桁の各要素を10進5桁に収める。
  { int    i;
    double W, W1; 
// First Loop
      for(i=N2; i>=2; i--)
       { W       = IV[i] ;           // 10進5桁を超えたものは、一つ前の
         W1      = floor(W*D10) ;    // 要素に桁上げ処理をする。
         IV[i-1] = IV[i-1] + W1 ;
         IV[i]   = W - W1*P10 ; }
// Second  Loop
      for(i=N2; i>=2; i--)           // 再度正確な桁揃え　
       { if(IV[i] >= P10)               // 10**5以上の要素の処理
          { IV[i]   = IV[i] - P10 ;
            IV[i-1] = IV[i-1] + 1 ; }
         if(IV[i] < 0)                  // 負要素の処理
          { IV[i]   = IV[i] + P10 ;
            IV[i-1] = IV[i-1] - 1 ;}
        }
   }
//===================================================================C
// Main Program                                                      C
//===================================================================C
 void main()
 { int    i, k, ID, NN, N, IX, KETA, T1, T2, T3, T4, KEND;
   double TT3, TT4;
   char   FNAME[] = "AtanOrg?.txt" ;
   char   FID[]   = "12345"	;
   FILE   *FT1;
//   Constant Set
      L10 = L10V ;                 // L10 = 5
      P10 = pow(10.0, L10) ;       // P10 = 10**5
      D10 = 1.0/P10 ;              // D10 = 1/10**5
//   Parameter Read
      printf(" Type In ID, (ID=1,2,3,4,5) \n") ;
      printf(" 1;Machin, 2;Gauss, 3;Klingen., 4;Stoemer, 5;Takano \n") ;
      scanf("%d",&ID) ;           // 計算公式選択
      if(ID <= 0 || ID >= 6) ID = 1 ;
      printf(" Type In N (N=Number of Pi Decimal-Digits) \n") ;
      scanf("%d",&NN) ;           // 計算桁数入力　
      N = NN/L10 ;
      if(N >= ND-3) N = ND - 3 ;
      NN = N*L10 ;
      N2 = N + 3 ;                //  計算要素数を求める(誤差を考慮して）
//  Open Output File
      ID = ID - 1;
      FNAME[7] = FID[ID] ;
      FT1 = fopen(FNAME,"w") ;
      printf(" ** Pi Computation Start ** \n") ;
      printf("  Number of Pi Decimal-Digits= %d \n",NN) ;
	  fprintf(FT1," ** Pi Computation with the original arctan method ** \n") ;
      fprintf(FT1,"  Number of Pi Decimal-Digits= %d, ID= %d \n",NN,ID+1) ;
      T1 = clock(); 
//  Clear Pi Area
      for(i=1; i <= N2; i++) IV[i] = 0 ;
//  Pi Computation Start
      for(k=0; k < INO[ID]; k++)     //  Atanの個数だけ反復計算
       { T2 = clock() ;
         M  = IDATA[ID][k][0] ;      //  Set M
         X  = IDATA[ID][k][1] ;      //  Set X
         ATANV() ;                   //  IV = IV + M*Atan(1/X) ---> π
         T3  = clock() ;
         TT3 = (double)(T3 - T2)/CLOCKS_PER_SEC ; 
         IX = X ;
         printf(" Atan(1/ %d ) Computation,  Time(s)= %f \n",IX,TT3) ;
       }
//  Output 
      NORMAL() ;             //  各要素を10進5桁の値にそろえる
      T4  = clock() ;
      TT4 = (double)(T4-T1)/CLOCKS_PER_SEC ;     // πの計算時間
      printf(" Pi Computation,  Time(s)= %f \n",TT4) ;
	  fprintf(FT1," Pi Computation,  Time(s)= %f \n\n",TT4) ;
      fprintf(FT1," Pi = 3.") ;
      for (k=2; k<=N2-2; k=k+16)      //  Output Pi Value
       { KEND = k+15 ;                //  80 Decimal in a Line
         if(KEND > N2-2) KEND = N2-2 ;
         KETA = (k-2)*L10 + 1 ;
         fprintf(FT1,"\n %8d",KETA) ;
         for (i=k; i<=KEND; i=i+2)
          { fprintf(FT1," %5.5d%5.5d",IV[i],IV[i+1]) ;  }
       }
  }