DRM法によるATAN公式π計算Cプログラム
2003/01/10 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------
1. 概要
DRM法によるATAN公式π計算ルーチンである。桁数nに対して、演算量はn･log(n)2〜3
に比例する。πを求める公式の2.1 Arctan公式 | を参照。
本計算には多数桁計算部品としてCの
乗除算 | 、
加減算 | 及び
Cのπ用I/O |
のルーチンが必要である。
2. プログラム
//=================================================================C
// DRM公式によるπ計算プログラム ( ATAN 版 , C )
// 2003/01/25 Yasunori Ushiro ( 後 保範 )
//=================================================================C
// Pi Value Computation Program by Atan with DRM Method C
// Pi = M1*Atan(1/X1) + M2*Atan(1/X2) + M3*Atan(1/X3) C
// Atan(1/X) = 1/X - 1/(3*X**3) + 1/(5*X**5) - 1/(7*X**7) +.. 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], E[ND], X[ND], Y[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 Atan
int INO[5] = {2,3,3,4,4}; // Number of Atan for Pi
// Constant of Atan formulation ; (M,X)
int IDATA[5][4][2] = {
16, 5, -4, 239, 0,0, 0,0, // ID=1
48,18, 32, 57, -20,239, 0,0, // ID=2
32,10, -4, 239, -16,515, 0,0, // ID=3
176,57, 28,239, -48,682, 96,12943, // ID=4
48,49, 128,57, -20,239, 48,110443 } ; // ID=5
// Function Definition in the Program
double DRMATN(int IX, int NN, int MLP) ;
double DRMAT2(int IX, int NN) ;
double DRMAT1(int N, int M) ;
// 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 the 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 the 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 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 IX, IB, T1, T2, T3 ;
double XX, ALG, ERR, ER1, TT ;
char FNAME[] = "AtanDRM?.txt" ;
char FID[] = "12345" ;
// Parameter Read and Set
ID = PIIN (2, ND, &N, &NC, &NN, &L10) ;
T1 = clock() ; // Start Time
ERR = 0.0 ;
SETV(0, PI, N+1) ; // PI=float(0)
// Main Loop
for (k=0; k<INO[ID]; k++) {
T2 = clock() ;
IX = IDATA[ID][k][1] ; // Get IX ; ATAN(1/IX)
IB = IDATA[ID][k][0] ; // Get IB ; IB*ATAN(1/IX)
// Set MLP
XX = IX ;
ALG = log(10.0)/log(XX) ;
MLP = 2.0*L10*NC*ALG/N + 1 ;
// Y=ATAN(1/IX)
ER1 = DRMATN(IX,N,MLP) ; // Compute Y=ATAN(1/IX)
ERR = fmax(ERR,ER1) ;
// Z=Z+IB*Y
FCMT(Y, IB, X, N) ; // X=IB*Y
FADD(PI, X, Y, N, 1) ; // Y=Pi+X
COPY(Y, PI, N+1) ; // Pi=Pi+IB*Y
T3 = clock() ;
TT = (double)(T3 - T2)/CLOCKS_PER_SEC ;
printf(" ATAN(1/%6d), DRM Step=%3d, Time=%8.2f (s)\n"
,IX,MLP,TT);
}
// Output Pi
PIOUT(PI, NC, L10, T1, ERR) ;
}
//==================================================================C
double DRMATN(int IX, int NN, int MLP)
//==================================================================C
// Compute Y=ATAN(1/IX) by DRM Method C
// Set Y,E,VS and Loop for MLP C
// Compute Y=Y+A/(B*E) C
// A/B=a1/b1+a2/(b2*X)+a3/(b3*X**2)+... by DRM C
// E = X**NN C
//------------------------------------------------------------------C
// ARGUMENTS : C
// IX I*4, In, Value IX for ATAN(1/IX) C
// NN I*4, In, Number of Data Gloup C
// MLP I*4, In, Number of Loops C
// ERR R*8, Out, Maximum Rounding Error to Convert Int. C
//------------------------------------------------------------------C
// Written by Yasunori Ushiro, 2003/01/25 C
// (Hitachi Ltd. and Waseda University) C
//==================================================================C
// Set Y,E
{ int k ;
double ERR, ER1 ;
SETV(0, Y, NN+1) ; // Y=0
FVSET(1, E, NN) ; // E=Float(1)
FCDIV(E, IX, E, NN) ; // E=Float(1/IX)
VS = 1.0 ;
// Main Loop
ERR = 0.0 ;
for (k=1; k<=MLP; k++) {
// Compute Y=Y+A/(B*E) with A,B,X,E by DRM
ER1 = DRMAT2(IX, NN) ;
ERR = fmax(ERR, ER1) ; }
return (ERR) ; }
//==================================================================C
double DRMAT2(int IX, int NN)
//==================================================================C
// Outer Routine of DRM Method for ATAN C
// Set Initial Value to A,B,X C
// Compute A,B,X by DRM and E=E*X C
// Y=Y+A/(B*E) by Floating C
//------------------------------------------------------------------C
// ARGUMENTS : C
// IX R*8, In, Value IX for ATAN(1/IX) 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
// Set Initial Value to A,B,X
{ int k, N, M, MEP ;
double AV, VX, ERR, ER1, ER2 ;
N = 4 ; // Element Size of Initial A,B
SETV(0, &A[1], NN) ; // Clear A
SETV(0, &B[1], NN) ; // Clear B
AV = 1.0 ;
for (k=1; k<=NN/N; k++) {
B[k*N] = VS ;
NORM(&B[k*(N-1)+1], N) ; // B(k)=VS+2*(k-1)
A[k*N] = AV ; // A(k)=(-1)**(k-1)
VS = VS + 2.0 ;
AV = -AV ;
/* fprintf(FT1,"A= %f %f %f %f \n",A[k*N-3],A[k*N-2],A[k*N-1],A[k*N]);
fprintf(FT1,"B= %f %f %f %f \n",B[k*N-3],B[k*N-2],B[k*N-1],B[k*N]); */
}
SETV(0, &X[1], N) ; // Clear X
VX = IX ; // Value of X in Atan(1/X)
X[N] = VX*VX ; // X=X**2
NORM(&X[1], N) ; // Fixed Normalize
// Compute A,B,X by DRM
M = NN/(2*N) ; // Half Number of Elements in A,B
MEP = LOG2F(2*M) ; // MEP is Step Number in DRM
ERR = 0.0 ;
for (k=1; k<=MEP; k++) {
ER1 = DRMAT1(N, M) ; // Compute A,B,X by DRM
ERR = fmax(ERR,ER1) ;
N = N*2 ;
M = M/2 ; }
// Change Floating A,B,X from Fixed
FLOAT(A, NN) ;
FLOAT(B, NN) ;
FLOAT(X, NN) ;
// E=E*X
ER1 = FMTCFM(X, E, W, NN, WK) ; // W=E*X
COPY(W, E, NN+1) ; // E=E*X
// Y=Y+A/(B*E)
ER2 = FMTCFM(B, E, W, NN, WK) ; // W=B*E
COPY(W, X, NN+1) ; // X=B*E
ERR = fmax3(ERR, ER1, ER2) ;
ER1 = FDIV(A, X, B, NN, W, WK) ; // B=A/X
FADD(Y, B, W, NN, 1) ; // W=Y+A/X
COPY(W, Y, NN+1) ; // Y=Y+A(B*E)
ERR = fmax(ERR, ER1) ;
return (ERR) ; }
//==================================================================C
double DRMAT1(int N, int M)
//==================================================================C
// Inner Routine of DRM Method for ATAN C
// W = X(1)*B(2,k) C
// A(1-2,k) = A(1,k)*W + A(2,k)*B(1,k) C
// B(1-2,k) = B(1,k)*B(2,k) for k=1,M C
// X(1-2) = X(1)*X(1) 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++) {
// W=X*B(2,k)
NK2 = NK1 + N ;
ER1 = MTCFM(&X[1], &B[NK2], &W[1], N, WK) ;
if(W[N] != 0.0)
{ printf(" ** Error Stop \n" ) ; }
// A(1-2,k)=A(1,k)*W+A(2,k)*B(1,k)
ER2 = MTCFM(&A[NK1], &W[NW2], &W[NW3],N,WK) ; // W(3)=A(1)*W(1)
ERR = fmax3(ERR, ER1, ER2) ;
ER1 = MTCFM(&A[NK2], &B[NK1], &W[1],N,WK) ; // W(1)=A(2)*B(1)
ADD(&W[1], &W[NW3], &A[NK1], 2*N, 1) ; // A=W(1)+W(3)
// B(1-2,k)=B(1,k)*B(2,k)
ER2 = MTCFM(&B[NK1], &B[NK2], &W[1],N,WK) ; // W=B(1)*B(2)
COPY(&W[1], &B[NK1], 2*N) ; // B=W
ERR = fmax3(ERR, ER1, ER2) ;
NK1 = NK1 + 2*N ; }
// X(1-2)=X*X
COPY(&X[1], &W[1], N) ; // W=X
ER1 = MTCFM(&W[1], &W[1], &X[1], N, WK) ; // X=W*W
ERR = fmax(ERR, ER1) ;
return (ERR); }