AGM法のπ計算Cプログラム
2003/01/10 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------
1. 概要
AGM法を利用したπ計算ルーチンである。桁数nに対して、演算量ははn･log(n)2に
比例する。πを求める公式の3.1 ガウス・ルジャンドル公式 | を参照。
本計算には多数桁計算部品として
Cの乗除算 | 、
Cの加減算 | 及び
Cのπ用I/O |
のルーチンが必要である。
2. プログラム
//=================================================================C
// AGM公式によるπ計算プログラム(O(N*Log2(N)**2)版 , FORTRAN
// 2003/01/06 Yasunori Ushiro ( 後 保範 )
//=================================================================C
// Pi Value Computation Program by AGM ( Gauss-Legendre Method ) C
// A=1, B=1/SQRT(2), T=1/4, S=1 C
// do k=1,log2(N) C
// W=A, A=(W+B)/2, B=SQRT(W*B) C
// T=T-S*(A-W)**2, S=2*S C
// end do C
// Pi = (A+B)**2/(4*T) C
// Using FMT for High-Precision Multiplication C
//-----------------------------------------------------------------C
// Written by Yasunori Ushiro, 2002/09/08 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], S[ND], PI[ND], W[ND] ;
double C[ND], FW[4*ND], WK[6*ND] ;
double BASE=1.0e5;
int NORID=-1;
FILE *FT1;
// Function Definition in Pi
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 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 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 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, NC, NN, L10, LOOP ;
int T1, T2, T3 ;
double ERR, ERR1, ERR2, ERR3, ERR4, TT ;
// Parameter Read and Set
ID = PIIN (1, ND, &N, &NC, &NN, &L10) ;
// Initial Value Set (A=1, B=1/SQRT(2), T=1/4, S=1)
T1 = clock() ; // Start Time
FVSET(1, A, N) ; // A=1
ERR = FCSQRT(2, B, N,-1,FW,WK) ; // B=1/SQRT(2)
FCDIV(A, 4, PI, N) ; // T=1/4, (T:PI)
FVSET(1, S, N) ; // S=1
// Iteration until Conversion
LOOP = LOG2F(N*L10) ;
for (k=1; k<=LOOP; k++) {
// W=A, A=(W+B)/2, B=SQRT(W*B)
T2 = clock() ; // Time set
COPY(A, W, N+1) ; // W=A
FADD(W, B, FW, N, 1) ; // FW=W+B
FCDIV(FW, 2, A, N) ; // A=(W+B)/2
ERR1 = FMTCFM(W, B, FW ,N,WK) ; // FW=W*B
COPY(FW, C, N+1) ; // C=W*B
ERR2 = FSQRT(C, B, N,1,FW,WK) ; // B=SQRT(W*B)
ERR = fmax3(ERR, ERR1, ERR2) ;
// T=T-S*(A-W)**2, S=2*S
FADD(A, W, C, N, -1) ; // C=A-W
ERR3 = FMTCFM(C,C, FW, N,WK) ; // FW=C**2
COPY(FW, C, N+1) ; // C=(A-W)**2
ERR4 = FMTCFM(S,C, FW, N,WK) ; // FW=S*(A-W)**2
COPY(PI, C, N+1) ;
FADD(C, FW, PI, N, -1) ; // T=T-S*(A-W)**2
FCMT(S, 2, C, N) ; // C=2*S
COPY(C, S, N+1) ; // S=2*S
ERR = fmax3(ERR, ERR3, ERR4) ;
T3 = clock() ; // Time set
TT = (double)(T3 - T2)/CLOCKS_PER_SEC ;
printf(" AGM Step= %d Time= %f (s) \n",k,TT) ;
}
// Pi = (A+B)**2/(4*T)
FADD(A, B, W, N, 1) ; // W=(A+B)
ERR1 = FMTCFM(W, W, FW, N,WK) ; // FW=W**2
COPY(FW, W, N+1) ; // W=(A+B)**2
FCMT(PI, 4, C, N) ; // C=4*T, (T:PI)
ERR2 = FDIV(W,C, PI, N, FW,WK) ; // Pi=(A+B)**2/(4*T)
ERR = fmax3(ERR,ERR1,ERR2) ;
// Output Pi
PIOUT(PI, NC, L10, T1, ERR) ;
}