#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
void CHNGBS(int IN, int NCR, int NBASE, int IS)
//==========================================================C
{  int  i, k, IV, IDV, IC, IOUT[20] ;
   char *CTB1[63] = {
         "0","1","2","3","4","5","6","7","8","9"
        ,"a","b","c","d","e","f","g","h","i","j"
        ,"k","l","m","n","o","p","q","r","s","t"
        ,"u","v","w","x","y","z","A","B","C","D"
        ,"E","F","G","H","I","J","K","K","M","N"
        ,"O","P","Q","R","S","T","U","V","W","X"
        ,"Y","Z"," " } ; 
//  Change to Charcter
   IV = IN ;
   for (k=NCR-1; k>=0; k=k-1) {
     IDV     = IV/NBASE ;
     IC      = IV - IDV*NBASE ;
     IOUT[k] = IC + IS ;
     IV      = IDV ; }
//  Write
     for (i=0; i<NCR; i++)
      { fprintf(FT1,"%s",CTB1[IOUT[i]]) ; }
  }
//==================================================================C
int XPIIN(int MID, int ND, int *NW, int *NCW, int *NNW, int *L10W,
         int *LBASEW, int *ISBW )
//------------------------------------------------------------------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
//       LBASE  I*4, Out, Base Value (3=< LBASE =<62)               C
//       ISB    I*4, Out, 0: Char. 0 Start,  10; Char. a Start      C
//------------------------------------------------------------------C
//  Written by Y. Ushiro ( Waseda Uni. and Hitachi )  2003.02.15    C
//==================================================================C
{  int ID, N, NC, NN, L10, LBASE, ISB, NB, NDMY ;
   double B2L ; 
   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" } ;
//  Base Read
   printf(" Type In Base Value ( from 3 to 62 ) \n") ;
   scanf("%d",&LBASE) ;
//  ISB Read
   printf(" Type In 0 or 1 (0; 0 start, 1; a start) \n") ;
   scanf("%d",&ISB) ;
   if(ISB != 0) ISB = 1 ;
   ISB = ISB*10 ;
   if(LBASE < 3 || LBASE > 62-ISB) LBASE = 26 ;
//  Constant Set
   B2L   = log(LBASE*1.0)/log(10.0) ;
   L10   = 5.7/B2L ;
   BASE  = power(LBASE, L10) ;      
   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 %d-Base) \n",LBASE) ;
   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) ;
   printf("  Base=%d,  Digit-Number of a Elements=%d \n",LBASE,L10) ; 
   fprintf(FT1,"  Base=%d,  Digit-Number of a Elements=%d \n",LBASE,L10) ;    
   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
   *LBASEW=LBASE;    *ISBW=ISB; 
   return (ID) ;
 }
//==================================================================C
void XPIOUT(double *Z, int NC, int L10, int LBASE, int ISB,
           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
//       LBASE       I*4, In,  Base Value (3=< LBASE =<62)          C
//       ISB         I*4, In,  0: Char. 0 Start,  10; Char. a Start 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.02.15    C
//==================================================================C
{  int    i, k, T4, KEND, KETA, LP, NPI ; 
   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.") ;
   LP = 60/L10 ;
   for (k=1; k<=NC; k=k+LP)      
     { KEND = MIN(LP, NC-k+1) ;     //  80 Characters in a Line
       KETA = (k-1)*L10 + 1 ;
       fprintf(FT1,"\n %6d ",KETA) ;
       for (i=1; i<=KEND; i=i++)
        { NPI = Z[i+k] ;
          CHNGBS(NPI,L10,LBASE,ISB) ;  } }
    }