#include <stdio.h>
#include <math.h>
#include <time.h>
#include <stdlib.h>
#include <string.h>
// Global Value 
extern double BASE ;
extern int NORID ;
extern FILE *FT1 ;
// Function Definition
void NORM(double *A, int N) ;
int LOG2F(int N) ;
int MIN(int N1, int N2) ;
int power(int N1, int N2) ;
//==================================================================C
int PIIN(int MID, int ND, int *NW, int *NCW, int *NNW, int *L10W)
//------------------------------------------------------------------C
//   Read and Set for Pi Computation Parameter                      C
//       MID    I*4, In,  Computing Method                          C
//                        MID = 0 ; Original Atan(1/X)              C
//                              1 ; AGM Method                      C
//                              1 ; Atan(1/X) with the DRM          C
//                              2 ; Ramanujan with the DRM          C
//       ND     I*4, In,  Size of Work Dimension (ND >= N)          C
//       ID     I*4, Out, ID in the Computing Method                C
//       N      I*4, Out, Number of Computing Elements for Pi       C
//       NC     I*4, Out, Number of Output Elements for Pi          C
//       NN     I*4, Out, Computating Decimal-Digits of Pi          C
//       L10    I*4, Out, Number of Decimal-Digits in a Element     C
//------------------------------------------------------------------C
//  Written by Y. Ushiro ( Waseda Uni. and Hitachi )  2003.01.25    C
//==================================================================C
{  int ID, N, NC, NN, L10, NB, NDMY ;
   int NSIZE[4] = {5,1,5,3} ;           // Number of Case in MID
   char *MSG1[4] = { " Type In ID, (ID=1,2,3,4,5) ",
         " Type In ID, (ID=1)         ", " Type In ID, (ID=1,2,3,4,5) ",
         " Type In ID, (ID=1,2,3)     " } ;
   char *MSG2[4] = {
         " 1;Machin, 2;Gauss, 3;Klingen., 4;Stoemer, 5;Takano",
         " 1;Gauss-Legendre (Now Only)                       ",
         " 1;Machin, 2;Gauss, 3;Klingen., 4;Stoemer, 5;Takano",
         " 1;Chudnovsky, 2;Fast-Ramanujan 3;Free-Sqrt        " } ;
   char FNAME[13] ;    
   char *FILET[4] = { "AtanORG?.txt", "AGM-2nd?.txt",
                      "AtanDRM?.txt", "RamnDRM?.txt" } ;
   char *NAME[4] = { "Atan Original", "AGM Gauss-Legendre",
                     "DRM with Atan", "DRM with Ramanujan" } ;
//  Constant Set
   L10  = 5 ;                       //  L10=5
   BASE  = power(10, L10) ;          //  BASE=10**5
   NORID = -1 ;                      //  Nearest Rounding
   NDMY  = 4 ;
//  Check MID
   if(MID < 0 || MID > 4) MID = 0 ;
//  ID Read 
   printf("%s \n",MSG1[MID]) ; 
   printf("%s \n",MSG2[MID]) ;
   scanf("%d",&ID) ;                 // Read ID (Pi Formulation)
   if(ID <= 0 || ID > NSIZE[MID]) ID = 1 ;
//  Output File Open
   strcpy(FNAME,FILET[MID]) ;
   FNAME[7]   = '0'+ ID ;
   FT1 = fopen(FNAME,"w") ;
//  Pi Decimal-Digits Read
   printf(" Type In N (N=Number of Pi Decimal-Digits) \n") ;
   scanf("%d",&NN) ;                 //  Input Pi Decimal-Digits
   N  = (NN+L10-1)/L10 + NDMY ;
   NC = N - NDMY ;                   //  Computation Decimal-Digits
   NN = NC*L10 ;                     //  Reset Input NN 
   NB = LOG2F(NC) ;                  //  NB=log2(N)
   N  = power(2, NB) ;               //  N  : Number of Elements
   if(N-NDMY < NC) N = 2*N ;         //  N is Number of Elements
//  Parameter Write
   printf(" ** Pi Computation by %s \n",NAME[MID]) ;
   fprintf(FT1," ** Pi Computation by %s ** \n",NAME[MID]) ;
   printf("  ID=%d,  Pi Decimal=%d,  FMT Size=%d\n",ID,NN,N) ; 
   fprintf(FT1,"  ID=%d,  Pi Decimal=%d,  FMT Size=%d\n",ID,NN,N) ; 
   ID = ID - 1 ;
//  Dimension Size Check
   if(N > ND) 
     { printf(" ** Error Stop ** N > ND \n") ;
	   exit(1) ;   }                         
   *NW=N;  *NCW=NC;  *NNW=NN; *L10W=L10;    // Set Return Value 
   return (ID) ;
 }
//==================================================================C
void PIOUT(double *Z, int NC, int L10, int T1, double ERR)
//------------------------------------------------------------------C
//   Output Pi Value on File-1                                      C
//       PI(0:NC+3)  R*8, In,  Pi Decimal-Value with Floating       C
//       NC          I*4, In,  Number of Output Elements for PI     C
//       L10         I*4, In,  Number of Decimals in a Element      C
//       T1          I*4, In,  Start Time                           C
//       ERR         R*8, In,  Rounding Error                       C
//------------------------------------------------------------------C
//  Written by Y. Ushiro ( Waseda Uni. and Hitachi )  2003.01.25    C
//==================================================================C
{  int    i, k, T4, KEND, KETA, NPI[18] ; 
   double TT ;
     NORID = 1 ;                       //  Positive Rounding
     NORM(&Z[1],NC+3) ;                   //  Positive Normalize 
     T4 = clock() ;                    //  End Time
     TT = (double)(T4 - T1)/CLOCKS_PER_SEC ; 
     printf("  Total Time= %f (s)  Rounding Error= %e \n",TT,ERR) ;
     fprintf(FT1,"  Total Time= %f (s)  Rounding Error= %e \n",TT,ERR) ;
//  Output Pi Value in the File 
     fprintf(FT1,"  Pi = 3. ") ;
     for (k=1; k<=NC; k=k+16)      //  Output Pi Value
      { KEND = MIN(16, NC-k+1) ;     //  80 Decimal in a Line
        KETA = (k-1)*L10 + 1 ;
        fprintf(FT1,"\n %8d",KETA) ;
        for (i=1; i<=KEND; i=i+2)
         { NPI[i] = Z[i+k] ;  NPI[i+1] = Z[i+k+1] ;
           fprintf(FT1," %5.5d%5.5d",NPI[i],NPI[i+1]) ;  } }
    }