π計算I/O用Cプログラム
2003/01/30 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------
1. 概要
π計算用パラメータ入力と、πの値を出力するCのプログラム。
従来法方式によるAtan公式を利用する場合以外のπ計算ルーチンから使用される。
2. プログラム
#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]) ; } }
}