//=================================================================C
//  DRM公式によるπ計算プログラム ( ATAN 版 , C )
//                   2003/01/25 Yasunori Ushiro ( 後　保範 )
//=================================================================C
//  Pi Value Computation Program by Atan with DRM 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
//   Using FMT for High-Precision Multiplication                   C
//-----------------------------------------------------------------C
//    Written by Yasunori Ushiro,  2003/01/25                      C
//        (Hitachi Ltd. and Waseda University)                     C
//=================================================================C
#include <stdio.h>
#include <math.h>
#include <time.h>
// Global Value 
#define ND  270001
#define L10V 5
double A[ND], B[ND], E[ND], X[ND], Y[ND], PI[ND] ;
double W[4*ND], WK[6*ND] ;
int    NPI[ND+3] ; 
FILE   *FT1 ;
double BASE=1.0e5 ;
int    NORID=-1 ;
double VS ;                  // Number of Current Term in Atan
int INO[5] = {2,3,3,4,4};    // Number of Atan for Pi 
// Constant of Atan formulation ; (M,X)
int IDATA[5][4][2] = { 
	   16, 5,  -4, 239,   0,0,     0,0,           // ID=1
       48,18,   32, 57,  -20,239,  0,0,           // ID=2
       32,10,  -4, 239,  -16,515,  0,0,           // ID=3
       176,57,  28,239,  -48,682,  96,12943,      // ID=4
        48,49,  128,57,  -20,239,  48,110443 } ;  // ID=5
// Function Definition in the Program
double DRMATN(int IX, int NN, int MLP) ;
double DRMAT2(int IX, int NN) ;
double DRMAT1(int N, int M) ;
// Function Definition in the PIOUT
int PIIN(int MID, int NX, int *N, int *NC, int *NN, int *L10) ;
void PIOUT(double *PI, int NC, int L10, int T1, double ERR) ;
// Function Definition in the MULT
double MTCFM(double *A, double *B, double *C, int N, double *WK) ;
double FMTCFM(double *A, double *B, double *C,int N, double *WK) ;
double FDIV(double *A, double *B, double *C, int N, double *FW,
				 double *WK) ;
double FCSQRT(int IA, double *C, int N, int ID, double *FW,
			  double *WK) ;
double FSQRT(double *A, double *C, int N, int ID, double *FW,
			  double *WK) ; 
void CFMT4(double *A, int N, int ID, double *WK) ;
void CFMT(double *A, int N, int ID, double *WK) ;
void CFMTTB(int N2, int ID, double *TBL) ;
void CFMTSB(double *A, int M, int L, double *TB, double *WK) ;
// Function definition in the ADD
void ADD(double *A, double *B, double *C, int N, int ID) ;
void FADD(double *A, double *B, double *C, int N, int ID) ;
void FCADD(int IA, double *B, double *C, int N, int ID) ;
void FCMT(double *A, int IB, double *C, int N) ;
void FCDIV(double *A, int IB, double *C, int N) ;
void NORM(double *A, int N) ;
void FNORM(double *A, int N) ;
void FLOAT(double *A, int N) ;
void COPY(double *A, double *C, int N) ;
void SETV(int IV, double *C, int N) ;
void FVSET(int IV, double *C, int N) ;
double fmax(double A1, double A2) ;
double fmax3(double A1, double A2, double A3) ;
int MIN(int N1, int N2) ;
int LOG2F(int N) ;
int power(int N1, int N2) ;
//=================================================================C
//  Main Program                                                   C
//-----------------------------------------------------------------C
void main () 
{  int    k, ID, N, NN, NC, L10, MLP ;
   int    IX, IB, T1, T2, T3 ;
   double XX, ALG, ERR, ER1, TT ;
   char   FNAME[] = "AtanDRM?.txt" ;
   char   FID[]   = "12345" ;
//  Parameter Read and Set
   ID = PIIN (2, ND, &N, &NC, &NN, &L10) ;
   T1 = clock() ;                     // Start Time
   ERR = 0.0 ;
   SETV(0, PI, N+1) ;                 //  PI=float(0)
//  Main Loop
   for (k=0; k<INO[ID]; k++) {
      T2 = clock() ; 
      IX  = IDATA[ID][k][1] ;         //  Get IX ; ATAN(1/IX)
      IB  = IDATA[ID][k][0] ;         //  Get IB ; IB*ATAN(1/IX)
//    Set MLP
      XX  = IX ;
      ALG = log(10.0)/log(XX) ;
      MLP = 2.0*L10*NC*ALG/N + 1 ;
//    Y=ATAN(1/IX)
      ER1 = DRMATN(IX,N,MLP) ;        //  Compute Y=ATAN(1/IX)
      ERR = fmax(ERR,ER1) ;
//    Z=Z+IB*Y
      FCMT(Y, IB, X, N) ;             //  X=IB*Y
      FADD(PI, X, Y, N, 1) ;          //  Y=Pi+X
      COPY(Y, PI, N+1) ;              //  Pi=Pi+IB*Y
      T3 = clock() ;
      TT = (double)(T3 - T2)/CLOCKS_PER_SEC ; 
      printf("  ATAN(1/%6d),  DRM Step=%3d,  Time=%8.2f (s)\n"
              ,IX,MLP,TT);
     }
//  Output Pi
    PIOUT(PI, NC, L10, T1, ERR) ;
  }
//==================================================================C
double DRMATN(int IX, int NN, int MLP)
//==================================================================C
//   Compute Y=ATAN(1/IX) by DRM Method                             C
//      Set Y,E,VS and Loop for MLP                                 C
//        Compute Y=Y+A/(B*E)                                       C
//          A/B=a1/b1+a2/(b2*X)+a3/(b3*X**2)+...   by DRM           C
//          E  = X**NN                                              C
//------------------------------------------------------------------C
//   ARGUMENTS :                                                    C
//      IX        I*4, In,   Value IX for ATAN(1/IX)                C
//      NN        I*4, In,   Number of Data Gloup                   C
//      MLP       I*4, In,   Number of Loops                        C
//      ERR       R*8, Out,  Maximum Rounding Error to Convert Int. C
//------------------------------------------------------------------C
//    Written by Yasunori Ushiro,  2003/01/25                       C
//        (Hitachi Ltd. and Waseda University)                      C
//==================================================================C
//  Set Y,E
{  int    k ;
   double ERR, ER1 ;
   SETV(0, Y, NN+1) ;             //  Y=0
   FVSET(1, E, NN) ;              //  E=Float(1)
   FCDIV(E, IX, E, NN) ;          //  E=Float(1/IX)
   VS = 1.0 ;
//  Main Loop
   ERR = 0.0 ;
   for (k=1; k<=MLP; k++) {
//   Compute Y=Y+A/(B*E) with A,B,X,E by DRM
      ER1 = DRMAT2(IX, NN) ;
      ERR = fmax(ERR, ER1) ;  }
   return (ERR) ;  }
//==================================================================C
double DRMAT2(int IX, int NN)
//==================================================================C
//   Outer Routine of DRM Method for ATAN                           C
//      Set Initial Value to A,B,X                                  C
//      Compute A,B,X by DRM and E=E*X                              C
//      Y=Y+A/(B*E) by Floating                                     C
//------------------------------------------------------------------C
//   ARGUMENTS :                                                    C
//      IX        R*8, In,   Value IX for ATAN(1/IX)                C
//      NN        I*4, In,   Number of Data Gloup                   C
//------------------------------------------------------------------C
//    Written by Yasunori Ushiro,  2003/01/25                       C
//        (Hitachi Ltd. and Waseda University)                      C
//==================================================================C
//  Set Initial Value to A,B,X
{  int k, N, M, MEP ;
   double AV, VX, ERR, ER1, ER2 ;
   N = 4 ;                        //  Element Size of Initial A,B
   SETV(0, &A[1], NN) ;           //  Clear A
   SETV(0, &B[1], NN) ;           //  Clear B
   AV = 1.0 ;
   for (k=1; k<=NN/N; k++) {
      B[k*N] = VS ;
      NORM(&B[k*(N-1)+1], N) ;    //  B(k)=VS+2*(k-1)
      A[k*N] = AV ;               //  A(k)=(-1)**(k-1)
      VS     = VS + 2.0 ;
      AV     = -AV ; 
/*   fprintf(FT1,"A= %f %f %f %f \n",A[k*N-3],A[k*N-2],A[k*N-1],A[k*N]);
   fprintf(FT1,"B= %f %f %f %f \n",B[k*N-3],B[k*N-2],B[k*N-1],B[k*N]);  */
     }
   SETV(0, &X[1], N) ;            //  Clear X
   VX   = IX ;                    //  Value of X in Atan(1/X)
   X[N] = VX*VX ;                 //  X=X**2
   NORM(&X[1], N) ;               //  Fixed Normalize
//  Compute A,B,X by DRM
   M = NN/(2*N) ;                 //  Half Number of Elements in A,B
   MEP = LOG2F(2*M) ;             //  MEP is Step Number in DRM
   ERR = 0.0 ;
   for (k=1; k<=MEP; k++) {
      ER1 = DRMAT1(N, M) ;        //  Compute A,B,X by DRM
      ERR = fmax(ERR,ER1) ;
      N = N*2 ;
      M = M/2 ;  }
//  Change Floating A,B,X from Fixed
   FLOAT(A, NN) ;
   FLOAT(B, NN) ;
   FLOAT(X, NN) ;
//  E=E*X
   ER1 = FMTCFM(X, E, W, NN, WK) ;     // W=E*X
   COPY(W, E, NN+1) ;                  // E=E*X
//  Y=Y+A/(B*E)
   ER2 = FMTCFM(B, E, W, NN, WK) ;     // W=B*E
   COPY(W, X, NN+1) ;                  // X=B*E
   ERR = fmax3(ERR, ER1, ER2) ;
   ER1 = FDIV(A, X, B, NN, W, WK) ;    // B=A/X
   FADD(Y, B, W, NN, 1) ;              // W=Y+A/X
   COPY(W, Y, NN+1) ;                  // Y=Y+A(B*E)
   ERR = fmax(ERR, ER1) ;
   return (ERR) ;  }
//==================================================================C
double DRMAT1(int N, int M)
//==================================================================C
//   Inner Routine of DRM Method for ATAN                           C
//      W        = X(1)*B(2,k)                                      C
//      A(1-2,k) = A(1,k)*W + A(2,k)*B(1,k)                         C
//      B(1-2,k) = B(1,k)*B(2,k)       for k=1,M                    C
//        X(1-2)   = X(1)*X(1)                                      C
//------------------------------------------------------------------C
//   ARGUMENTS :                                                    C
//      N         I*4, In,   Number of Computation Elements         C
//      M         I*4, In,   Number of Data Gloup                   C
//------------------------------------------------------------------C
//    Written by Yasunori Ushiro,  2003/01/25                       C
//        (Hitachi Ltd. and Waseda University)                      C
//==================================================================C
{  int    k, NW2, NW3,NK1, NK2 ;
   double ERR, ER1, ER2 ; 
   ERR = 0.0 ;
   NW2 = N + 1;
   NW3 = 2*N + 1;
   NK1 = 1 ;
   for (k=1; k<=M; k++) {
//   W=X*B(2,k)
      NK2 = NK1 + N ;
      ER1 = MTCFM(&X[1], &B[NK2], &W[1], N, WK) ;
      if(W[N] != 0.0) 
	    { printf(" ** Error Stop \n" ) ; }
//   A(1-2,k)=A(1,k)*W+A(2,k)*B(1,k)
      ER2 = MTCFM(&A[NK1], &W[NW2], &W[NW3],N,WK) ;  // W(3)=A(1)*W(1)
      ERR = fmax3(ERR, ER1, ER2) ;
      ER1 = MTCFM(&A[NK2], &B[NK1], &W[1],N,WK) ;    // W(1)=A(2)*B(1)
      ADD(&W[1], &W[NW3], &A[NK1], 2*N, 1) ;         // A=W(1)+W(3)
//   B(1-2,k)=B(1,k)*B(2,k)
      ER2 = MTCFM(&B[NK1], &B[NK2], &W[1],N,WK) ;    // W=B(1)*B(2)
      COPY(&W[1], &B[NK1], 2*N) ;                    // B=W
      ERR = fmax3(ERR, ER1, ER2) ;
      NK1 = NK1 + 2*N ;   }
//  X(1-2)=X*X
   COPY(&X[1], &W[1], N) ;                           // W=X
   ER1 = MTCFM(&W[1], &W[1], &X[1], N, WK) ;         // X=W*W
   ERR = fmax(ERR, ER1) ;
   return (ERR);  }