//=================================================================C
//  AGM公式によるπ計算プログラム(O(N*Log2(N)**2)版 , FORTRAN
//　　　　　　　　　　　　 2003/01/06 Yasunori Ushiro ( 後　保範 )
//=================================================================C
//  Pi Value Computation Program by AGM ( Gauss-Legendre Method )  C
//     A=1, B=1/SQRT(2), T=1/4, S=1                                C
//     do k=1,log2(N)                                              C
//        W=A,  A=(W+B)/2,  B=SQRT(W*B)                            C
//        T=T-S*(A-W)**2,   S=2*S                                  C
//     end do                                                      C
//     Pi = (A+B)**2/(4*T)                                         C
//   Using FMT for High-Precision Multiplication                   C
//-----------------------------------------------------------------C
//    Written by Yasunori Ushiro,  2002/09/08                      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], S[ND], PI[ND], W[ND] ;
double C[ND], FW[4*ND], WK[6*ND] ;
double BASE=1.0e5;
int    NORID=-1;
FILE   *FT1;
// Function Definition in Pi 
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 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 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 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, NC, NN, L10, LOOP ;
   int    T1, T2, T3 ;
   double ERR, ERR1, ERR2, ERR3, ERR4, TT ;
//  Parameter Read and Set
   ID = PIIN (1, ND, &N, &NC, &NN, &L10) ;
//  Initial Value Set (A=1, B=1/SQRT(2), T=1/4, S=1)
   T1 = clock() ;                      // Start Time
   FVSET(1, A, N) ;                    // A=1
   ERR = FCSQRT(2, B, N,-1,FW,WK) ;    // B=1/SQRT(2)
   FCDIV(A, 4, PI, N) ;                // T=1/4, (T:PI)
   FVSET(1, S, N) ;                    // S=1
//  Iteration until Conversion
   LOOP = LOG2F(N*L10) ;
   for (k=1; k<=LOOP; k++) {
//    W=A,  A=(W+B)/2,  B=SQRT(W*B)
      T2 = clock() ;                    // Time set
      COPY(A, W, N+1) ;                 // W=A
      FADD(W, B, FW, N, 1) ;            // FW=W+B
      FCDIV(FW, 2, A, N) ;              // A=(W+B)/2
      ERR1 = FMTCFM(W, B, FW ,N,WK) ;   // FW=W*B
      COPY(FW, C, N+1) ;                // C=W*B
      ERR2 = FSQRT(C, B, N,1,FW,WK) ;   // B=SQRT(W*B)
      ERR = fmax3(ERR, ERR1, ERR2) ;
//    T=T-S*(A-W)**2,   S=2*S  
      FADD(A, W, C, N, -1) ;            // C=A-W
      ERR3 = FMTCFM(C,C, FW, N,WK) ;    // FW=C**2
      COPY(FW, C, N+1) ;                // C=(A-W)**2
      ERR4 = FMTCFM(S,C, FW, N,WK) ;    // FW=S*(A-W)**2
      COPY(PI, C, N+1) ; 	              
      FADD(C, FW, PI, N, -1) ;          // T=T-S*(A-W)**2
      FCMT(S, 2, C, N) ;                // C=2*S
      COPY(C, S, N+1) ;                 // S=2*S
      ERR = fmax3(ERR, ERR3, ERR4) ;
      T3 = clock() ;                    // Time set
      TT = (double)(T3 - T2)/CLOCKS_PER_SEC ; 
      printf("  AGM Step= %d  Time= %f (s) \n",k,TT) ;
     }
//  Pi = (A+B)**2/(4*T)
   FADD(A, B, W, N, 1) ;               // W=(A+B)
   ERR1 = FMTCFM(W, W, FW, N,WK) ;     // FW=W**2
   COPY(FW, W, N+1) ;                  // W=(A+B)**2
   FCMT(PI, 4, C, N) ;                 // C=4*T, (T:PI) 
   ERR2 = FDIV(W,C, PI, N, FW,WK) ;    // Pi=(A+B)**2/(4*T)
   ERR = fmax3(ERR,ERR1,ERR2) ;
//  Output Pi
   PIOUT(PI, NC, L10, T1, ERR) ;
  }