//=================================================================C
//  DRM公式によるπ計算プログラム ( Ramanujan 版 , C )
//                      2003/01/25 Yasunori Ushiro ( 後　保範 )
//=================================================================C
//    Pi Value Computation Program by Ramanujan with DRM Method    C
//        Pi = K0*K1*SQRT(K2)/(K3*(Y+K4))                          C
//          Y = 0                                                  C
//          do i=1,ML                                              C
//            Y = (C + A*Y)/B  with A,B,C by DRM Method            C
//          end do                                                 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], C[ND], X[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 Ramanujan
//  Pi Ramanujan Constant
int KD[3][6] = {
//       K0   K1    K2  K3     K4    K5
        128,10005,10005, 3, 13591409, 5,        // ID=1
          1, 9801,    2, 4,     1103, 3,        // ID=2
          1, 3528,    1, 1,     1123, 3 } ;     // ID=3
// Function definition in the Program
double DRMRMA(int ID, int NN) ;
double DRMRM1(int N, int M) ;
void VRMSET(double *A, double VS, int *IB, int N, int NO) ;
// 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 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 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 CMLT(double *A, double X, 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    K2, T1, T2, T3 ;
   double XX, ALG, ERR, ER1, ER2, VV, TT ;
   char   FNAME[] = "RamnDRM?.txt" ;
   char   FID[]   = "12345" ; 
//  Parameter Read and Set
   ID = PIIN (3, ND, &N, &NC, &NN, &L10) ;
//  Set MLP
   XX = 53360.0*53360.0*53360.0 ;
   if(ID+1==2) XX = 99.0*99.0*99.0*99.0 ;
   if(ID+1==3) XX = 882.0*882.0 ;
   ALG = log(10.0)/log(XX) ;
   MLP = 8.0*L10*NC*ALG/N + 1 ;
   printf("  DRM Steps= %d,  FMT Size= %d \n",MLP,N) ;
   fprintf(FT1,"  DRM Steps= %d,  FMT Size= %d \n",MLP,N) ;
   T1 = clock() ;                       // Start Time
   ERR = 0.0 ;
   ER2 = 0.0 ;
//  Main Loop
   VV = N/8 ;
   for (k=1; k<=MLP;k++) {
      VS = (MLP-k)*VV + 1.0 ;           // Start Value
      T2 = clock();
//   Computing A,B,C by DRM Method 
      ER1 = DRMRMA(ID,N) ;
	  ERR = fmax(ERR,ER1) ;
//   Y = (C + A*Y)/B or Y = C/B (k=1)
      if (k == 1) 
        { FVSET(0, W, N) ;
          W[0] = -N ; }                 // W=Floating(0)
      else 
	  { ER2 = FMTCFM(A,PI,W,N,WK); }    // W=A*Y, (Y:PI)
      FADD(C,W,X,N,1) ;                 // X=C+W
      ER1 = FDIV(X,B,PI,N,W,WK) ;       // Y=X/B
      ERR = fmax3(ERR,ER1,ER2) ;
      T3 = clock();
      TT = (double)(T3 - T2)/CLOCKS_PER_SEC ; 
      printf("  DRM Step=%3d,  Time=%8.2f (s)\n",k,TT) ;
	   }
//  Pi=F=K0*K1*SQRT(K2)/(K3*(K4+K5*Y))
   ER1 = 0.0 ;
   K2 = KD[ID][2] ;
   if (K2==1) 
     { FVSET(1,X,N) ; }                  // X=Floating(1)
   else
     { ER1 = FCSQRT(K2,X,N,1,W,WK) ; }   // X=SQRT(K2)
   FCMT(X, KD[ID][0], A, N) ;            // A=K0*X
   FCMT(A, KD[ID][1], B, N) ;            // B=K1*A
   FVSET(KD[ID][4], A, N) ;              // A=K4
   FCMT(PI, KD[ID][5], C, N) ;           // C=K5*Y, Y:PI
   FADD(A, C, X, N, 1) ;                 // X=A+C
   FCMT(X, KD[ID][3], C, N) ;            // C=K3*X
   ER2 = FDIV(B, C, PI, N, W, WK) ;      // Pi=B/C
   ERR = fmax3(ERR, ER1, ER2) ;
//  Output Pi
   PIOUT(PI, NC, L10, T1, ERR) ;
  }
//==================================================================C
double DRMRMA(int ID, int NN)
//==================================================================C
//   Outer Routine of DRM Method for Ramanujan                      C
//      Set Initial Value to A,B,C                                  C
//      Compute A,B,C by DRM and Floating                           C
//------------------------------------------------------------------C
//   ARGUMENTS :                                                    C
//      ID        I*4, In,   Type of Ramanujan (1,2,3)              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
{  int    k, kk, N, M, MEP ;
   double ERR, ER1, SG, VS1 ;
//  A-Data
   int IKA[3][6] = { 2,-1,6,-5,6,-1, 2,-1,4,-3,4,-1, 2,-1,4,-3,4,-1 } ;
//  B-Data
   int IKB[3][12] = { 3,0, 3,0, 1,0, 0,106720, 0,106720, 0,106720,
                     2,0, 1,0, 1,0, 0,39204,  0,198,    0,198,
                     2,0, 1,0, 1,0, 0,3528,   0,3528,   0,1 } ;
//  C-Data
   double CK[3][2] = { 545140134.0, 13591409.0, 26390.0, 1103.0,
                      21460.0, 1123.0 } ;
//  Set Initial Value to A,B,C
   N  = 8 ;                     // Number of Elements for A,B,C
   kk = 1 ;                     // Start Address of A,B,C
   SG = 1 ;                       // Sign of C
   for (k=1; k<=NN/N; k++) {
     if(ID+1 != 2) SG = -SG ;     // Change Sign of C 
//   Set A-Data
      VS1 = VS + 1.0 ;
	  VRMSET(&A[kk], VS1, &IKA[ID][0], N,3) ;    // Set A
//   Set B-Data
      VRMSET(&B[kk], VS, &IKB[ID][0], N,6) ;     // Set B
//   Set C-Data 
      SETV(0, &C[kk], N) ;                   //  Clear C
      C[kk+N-1] = VS ;                       //  C=VS
      NORM(&C[kk], N) ;                      //  Normalize C
      CMLT(&C[kk], CK[ID][0],N) ;            //  C=C*IKC(1)
      C[kk+N-1] = C[kk+N-1] + CK[ID][1] ;    //  C=C+IKC(2)
      NORM(&C[kk], N) ;	  	  	             //  Normalize C
      CMLT(&C[kk], SG, N) ;                  //  C=C*IS
      NORM(&C[kk], N) ;	  	  	             //  Normalize C
//   Count Up
      VS = VS + 1.0 ;
      kk = kk + N ;
     }
//  Compute A,B,X by DRM
   M = NN/(2*N) ;
   MEP = LOG2F(2*M) ;                   // MEP=log2(M)
   ERR = 0.0 ;
   for (k=1; k<=MEP; k++) {
      ER1 = DRMRM1(N, M) ;              // Compute A,B,C by DRM
      ERR = fmax(ERR,ER1) ;
      N = N*2 ;
      M = M/2 ;  }
//  Change Floating A,B,C from Fixed
   FLOAT(A, NN) ;
   FLOAT(B, NN) ;
   FLOAT(C, NN) ;
   return (ERR) ;  } 
//==================================================================C
void VRMSET(double *A, double VV, int *IB, int N, int NO)
//==================================================================C
//   Set Initial Value for A-Data and B-Data                        C
//   ARGUMENTS :                                                    C
//      A(N)      R*8, I/O,  Input and Output Data                  C
//      VV        R*8, I/O,  Start Value of Computation Terms       C
//      IB(2,NO)  I*4, In,   Constant Table                         C
//      N         I*4, In,   Number of Elements in A                C
//      NO        I*4, In,   Number of Constant Table               C
//------------------------------------------------------------------C
//    Written by Yasunori Ushiro,  2003/01/25                       C
//        (Hitachi Ltd. and Waseda University)                      C
//==================================================================C
{  int k ;   double X ;
   SETV(0, A, N) ;                     //  Clear A
   A[N-1] = 1.0 ;                      //  A=1 by fixed
   for (k=0; k<NO; k++) {
      X = IB[2*k]*VV + IB[2*k+1] ;     //  X=c1*VS+c2
      CMLT(A, X ,N) ;    }             //  A=A*X by fixed
  }
//==================================================================C
double DRMRM1(int N, int M)
//==================================================================C
//   Inner Routine of DRM Method for Ramanujan                      C
//      A(1-2,k) = A(1,k)*A(2,k)                                    C
//      B(1-2,k) = B(1,k)*B(2,k)                                    C
//      C(1-2,k) = A(1,k)*C(2,k) + B(2,k)*C(1,k) , for k=1,M        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++) {
      NK2 = NK1 + N ;
      ER1 = MTCFM(&A[NK1], &C[NK2], &W[1], N, WK) ;     // W1=A1*C2
      ER2 = MTCFM(&B[NK2], &C[NK1], &W[NW3], N, WK) ;   // W3=B2*C1
      ERR = fmax3(ERR, ER1, ER2) ;  
      ADD(&W[1], &W[NW3], &C[NK1], 2*N, 1) ;            // C=W1+W3
//   A(1-2,k) = A(1,k)*A(2,k) 
      ER1 = MTCFM(&A[NK1], &A[NK2], &W[1], N, WK) ;     // W1=A1*A2
      COPY(&W[1], &A[NK1], 2*N) ;                       // A=W1 
//   B(1-2,k) = B(1,k)*B(2,k)
      ER2 = MTCFM(&B[NK1], &B[NK2], &W[1], N, WK) ;     // W1=B1*B2
      COPY(&W[1], &B[NK1], 2*N) ;                       // B=W1
      ERR = fmax3(ERR, ER1, ER2) ;
	  NK1 = NK1 + 2*N ;
     }
   return (ERR);  }