各種基底によるπ計算I/O用Cプログラム
2003/02/15 日立製作所 & 早稲田大学 後 保範 ( 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
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
LP = 60/L10 ;
for (k=1; k<=NC; k=k+LP)
{ KEND = MIN(LP, NC-k+1) ; // 60 Characters in a Line
KETA = (k-1)*L10 + 1 ;
fprintf(FT1,"\n %8d ",KETA) ;
for (i=1; i<=KEND; i=i++)
{ NPI = Z[i+k] ;
CHNGBS(NPI,L10,LBASE,ISB) ; } }
}