DRM法によるラマヌジャン型公式π計算Cプログラム
2003/01/10 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------
1. 概要
DRM法でのラマヌジャン型π計算ルーチンである。桁数nに対して、演算量はn･log(n)2〜3
に比例する。πを求める公式の2.2 ラマヌジャン型公式 | を参照。
本計算には多数桁計算部品としてCの
乗除算 | 、
加減算 | 及び
Cのπ用I/O |
のルーチンが必要である。
2. プログラム
//=================================================================C
// DRM公式によるπ計算プログラム ( Ramanujan 版 , C )
// 2003/01/25 Yasunori Ushiro ( 後 保範 )
//=================================================================C
// Pi Value Computation Program by Ramanujan with DRM Method C
// Pi = K0*K1*SQRT(K2)/(K3*(Y+K4)) C
// Y = 0 C
// do i=1,ML C
// Y = (C + A*Y)/B with A,B,C by DRM Method C
// end do C
// Using FMT for High-Precision Multiplication C
//-----------------------------------------------------------------C
// Written by Yasunori Ushiro, 2003/01/25 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], C[ND], X[ND], PI[ND] ;
double W[4*ND], WK[6*ND] ;
int NPI[ND+3] ;
FILE *FT1 ;
double BASE=1.0e5 ;
int NORID=-1 ;
double VS ; // Number of Current Term in Ramanujan
// Pi Ramanujan Constant
int KD[3][6] = {
// K0 K1 K2 K3 K4 K5
128,10005,10005, 3, 13591409, 5, // ID=1
1, 9801, 2, 4, 1103, 3, // ID=2
1, 3528, 1, 1, 1123, 3 } ; // ID=3
// Function definition in the Program
double DRMRMA(int ID, int NN) ;
double DRMRM1(int N, int M) ;
void VRMSET(double *A, double VS, int *IB, int N, int NO) ;
// Function Definition in the PIOUT
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 MTCFM(double *A, double *B, double *C, int N, double *WK) ;
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 ADD(double *A, double *B, double *C, int N, int ID) ;
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 CMLT(double *A, double X, int N) ;
void NORM(double *A, int N) ;
void FNORM(double *A, int N) ;
void FLOAT(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, NN, NC, L10, MLP ;
int K2, T1, T2, T3 ;
double XX, ALG, ERR, ER1, ER2, VV, TT ;
char FNAME[] = "RamnDRM?.txt" ;
char FID[] = "12345" ;
// Parameter Read and Set
ID = PIIN (3, ND, &N, &NC, &NN, &L10) ;
// Set MLP
XX = 53360.0*53360.0*53360.0 ;
if(ID+1==2) XX = 99.0*99.0*99.0*99.0 ;
if(ID+1==3) XX = 882.0*882.0 ;
ALG = log(10.0)/log(XX) ;
MLP = 8.0*L10*NC*ALG/N + 1 ;
printf(" DRM Steps= %d, FMT Size= %d \n",MLP,N) ;
fprintf(FT1," DRM Steps= %d, FMT Size= %d \n",MLP,N) ;
T1 = clock() ; // Start Time
ERR = 0.0 ;
// Main Loop
VV = N/8 ;
for (k=1; k<=MLP;k++) {
VS = (MLP-k)*VV + 1.0 ; // Start Value
T2 = clock();
// Computing A,B,C by DRM Method
ER1 = DRMRMA(ID,N) ;
// Y = (C + A*Y)/B or Y = C/B (k=1)
if (k == 1)
{ FVSET(0, W, N) ;
W[0] = -N ; } // W=Floating(0)
else
{ ER2 = FMTCFM(A,PI,W,N,WK); } // W=A*Y, (Y:PI)
ERR = fmax3(ERR,ER1,ER2) ;
FADD(C,W,X,N,1) ; // X=C+W
ER1 = FDIV(X,B,PI,N,W,WK) ; // Y=X/B
ERR = fmax(ERR,ER1) ;
T3 = clock();
TT = (double)(T3 - T2)/CLOCKS_PER_SEC ;
printf(" DRM Step=%3d, Time=%8.2f (s)\n",k,TT) ;
}
// Pi=F=K0*K1*SQRT(K2)/(K3*(K4+K5*Y))
K2 = KD[ID][2] ;
if (K2==1)
{ FVSET(1,X,N) ; } // X=Floating(1)
else
{ ER1 = FCSQRT(K2,X,N,1,W,WK) ; } // X=SQRT(K2)
FCMT(X, KD[ID][0], A, N) ; // A=K0*X
FCMT(A, KD[ID][1], B, N) ; // B=K1*A
FVSET(KD[ID][4], A, N) ; // A=K4
FCMT(PI, KD[ID][5], C, N) ; // C=K5*Y, Y:PI
FADD(A, C, X, N, 1) ; // X=A+C
FCMT(X, KD[ID][3], C, N) ; // C=K3*X
ER2 = FDIV(B, C, PI, N, W, WK) ; // Pi=B/C
ERR = fmax3(ERR, ER1, ER2) ;
// Output Pi
PIOUT(PI, NC, L10, T1, ERR) ;
}
//==================================================================C
double DRMRMA(int ID, int NN)
//==================================================================C
// Outer Routine of DRM Method for Ramanujan C
// Set Initial Value to A,B,C C
// Compute A,B,C by DRM and Floating C
//------------------------------------------------------------------C
// ARGUMENTS : C
// ID I*4, In, Type of Ramanujan (1,2,3) C
// NN I*4, In, Number of Data Gloup C
//------------------------------------------------------------------C
// Written by Yasunori Ushiro, 2003/01/25 C
// (Hitachi Ltd. and Waseda University) C
//==================================================================C
{ int k, kk, N, M, MEP ;
double ERR, ER1, SG, VS1 ;
// A-Data
int IKA[3][6] = { 2,-1,6,-5,6,-1, 2,-1,4,-3,4,-1, 2,-1,4,-3,4,-1 } ;
// B-Data
int IKB[3][12] = { 3,0, 3,0, 1,0, 0,106720, 0,106720, 0,106720,
2,0, 1,0, 1,0, 0,39204, 0,198, 0,198,
2,0, 1,0, 1,0, 0,3528, 0,3528, 0,1 } ;
// C-Data
double CK[3][2] = { 545140134.0, 13591409.0, 26390.0, 1103.0,
21460.0, 1123.0 } ;
// Set Initial Value to A,B,C
N = 8 ; // Number of Elements for A,B,C
kk = 1 ; // Start Address of A,B,C
SG = 1 ; // Sign of C
for (k=1; k<=NN/N; k++) {
if(ID+1 != 2) SG = -SG ; // Change Sign of C
// Set A-Data
VS1 = VS + 1.0 ;
VRMSET(&A[kk], VS1, &IKA[ID][0], N,3) ; // Set A
// Set B-Data
VRMSET(&B[kk], VS, &IKB[ID][0], N,6) ; // Set B
// Set C-Data
SETV(0, &C[kk], N) ; // Clear C
C[kk+N-1] = VS ; // C=VS
NORM(&C[kk], N) ; // Normalize C
CMLT(&C[kk], CK[ID][0],N) ; // C=C*IKC(1)
C[kk+N-1] = C[kk+N-1] + CK[ID][1] ; // C=C+IKC(2)
NORM(&C[kk], N) ; // Normalize C
CMLT(&C[kk], SG, N) ; // C=C*IS
NORM(&C[kk], N) ; // Normalize C
// Count Up
VS = VS + 1.0 ;
kk = kk + N ;
}
// Compute A,B,X by DRM
M = NN/(2*N) ;
MEP = LOG2F(2*M) ; // MEP=log2(M)
ERR = 0.0 ;
for (k=1; k<=MEP; k++) {
ER1 = DRMRM1(N, M) ; // Compute A,B,C by DRM
ERR = fmax(ERR,ER1) ;
N = N*2 ;
M = M/2 ; }
// Change Floating A,B,C from Fixed
FLOAT(A, NN) ;
FLOAT(B, NN) ;
FLOAT(C, NN) ;
return (ERR) ; }
//==================================================================C
void VRMSET(double *A, double VV, int *IB, int N, int NO)
//==================================================================C
// Set Initial Value for A-Data and B-Data C
// ARGUMENTS : C
// A(N) R*8, I/O, Input and Output Data C
// VV R*8, I/O, Start Value of Computation Terms C
// IB(2,NO) I*4, In, Constant Table C
// N I*4, In, Number of Elements in A C
// NO I*4, In, Number of Constant Table C
//------------------------------------------------------------------C
// Written by Yasunori Ushiro, 2003/01/25 C
// (Hitachi Ltd. and Waseda University) C
//==================================================================C
{ int k ; double X ;
SETV(0, A, N) ; // Clear A
A[N-1] = 1.0 ; // A=1 by fixed
for (k=0; k<NO; k++) {
X = IB[2*k]*VV + IB[2*k+1] ; // X=c1*VS+c2
CMLT(A, X ,N) ; } // A=A*X by fixed
}
//==================================================================C
double DRMRM1(int N, int M)
//==================================================================C
// Inner Routine of DRM Method for Ramanujan C
// A(1-2,k) = A(1,k)*A(2,k) C
// B(1-2,k) = B(1,k)*B(2,k) C
// C(1-2,k) = A(1,k)*C(2,k) + B(2,k)*C(1,k) , for k=1,M C
//------------------------------------------------------------------C
// ARGUMENTS : C
// N I*4, In, Number of Computation Elements C
// M I*4, In, Number of Data Gloup C
//------------------------------------------------------------------C
// Written by Yasunori Ushiro, 2003/01/25 C
// (Hitachi Ltd. and Waseda University) C
//==================================================================C
{ int k, NW2, NW3,NK1, NK2 ;
double ERR, ER1, ER2 ;
ERR = 0.0 ;
NW2 = N + 1;
NW3 = 2*N + 1;
NK1 = 1 ;
for (k=1; k<=M; k++) {
NK2 = NK1 + N ;
ER1 = MTCFM(&A[NK1], &C[NK2], &W[1], N, WK) ; // W1=A1*C2
ER2 = MTCFM(&B[NK2], &C[NK1], &W[NW3], N, WK) ; // W3=B2*C1
ERR = fmax3(ERR, ER1, ER2) ;
ADD(&W[1], &W[NW3], &C[NK1], 2*N, 1) ; // C=W1+W3
// A(1-2,k) = A(1,k)*A(2,k)
ER1 = MTCFM(&A[NK1], &A[NK2], &W[1], N, WK) ; // W1=A1*A2
COPY(&W[1], &A[NK1], 2*N) ; // A=W1
// B(1-2,k) = B(1,k)*B(2,k)
ER2 = MTCFM(&B[NK1], &B[NK2], &W[1], N, WK) ; // W1=B1*B2
COPY(&W[1], &B[NK1], 2*N) ; // B=W1
ERR = fmax3(ERR, ER1, ER2) ;
NK1 = NK1 + 2*N ;
}
return (ERR); }