// Atan公式によるπ計算プログラム(N**2版、Cプログラム) //                      2002/9/26 後 保範 //=================================================================C // Pi Value Computation Program by Atan Normal 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 //-----------------------------------------------------------------C // Written by Yasunori Ushiro, 2002/09/09 C // (Hitachi Ltd. and Waseda University) C //=================================================================C #include #include #include // Global Define #define ND 1000010 #define L10V 5 double X, P10, D10; int M, N2, L10; int IV[ND+1], IXD[ND+1]; int INO[5] = {2,3,3,4,4}; // π計算公式のAtanの個数 // π計算用各公式の定数 (M,X)の組、最大4個(Atanの個数) int IDATA[5][4][2] = { 16, 5, -4,239, 0,0, 0,0, 48,18, 32, 57, -20,239, 0,0, 32,10, -4,239, -16,515, 0,0, 176,57, 28,239, -48,682, 96,12943, 48,49, 128,57, -20,239, 48,110443}; //=================================================================C // IV = M*Atan(1/X) Computation by Normal Method C // IXD = 1/X**K C //-----------------------------------------------------------------C // Written by Yasunori Ushiro, 2002/09/08 C // (Hitachi Ltd. and Waseda University) C //=================================================================C void ATANV() // M*Atan(1/X)の多数桁の値の計算 { int i, k, KOSU, LST, ISG; double X2, DX2, W, W1, W2, VK, DK, WK, WK1, WK2; // Initial Constant X2 = X*X ; DX2 = 1.0/X2 ; KOSU = N2*L10/(log(X2)/log(10.0)); // 多数桁に必要なAtanの項数 // Initial Value Set for(i=2; i<=N2; i++) { IXD[i] = 0 ; } IXD[1] = M ; // Initial Computation, IV=IXD=M/X W2 = 0.0; for(i=1; i<=N2; i++) { W = IXD[i] + W2*P10 ; W1 = floor(W/X) ; W2 = W - W1*X ; IXD[i] = W1 ; IV[i] = IV[i] + IXD[i] ; } // IV = M/X + 前のAtanの計算値 // Main Computation LST = 1 ; ISG = -1 ; for (k=1; k<=KOSU; k++) // 計算項数だけ反復計算 { W2 = 0.0 ; WK2 = 0.0 ; VK = 2*k + 1 ; DK = 1.0/VK ; for(i=LST; i<=N2; i++) // 多数桁は1要素(int) { W = IXD[i] + W2*P10 ; // に10進5桁ごとに分けて記憶 W1 = floor(W*DX2) ; // floor: 切り捨てて整数値にする W2 = W - W1*X2 ; IXD[i] = W1 ; // W1 = 多数桁のM/X**(2k+1)の値 WK = W1 + WK2*P10 ; // W1 = 旧W1/X**2 = 旧W1*DX2 WK1 = floor(WK*DK) ; WK2 = WK - WK1*VK ; // WK2 = 多数桁のW1/(2k+1) IV[i] = IV[i] + ISG*WK1 ; } if(IXD[LST] == 0) { LST = LST + 1; // 無駄な計算を省くための処理 if(IXD[LST] == 0) LST = LST + 1; } ISG = -ISG ; } } //=================================================================C void NORMAL() // 多数桁の各要素を10進5桁に収める。 { int i; double W, W1; // First Loop for(i=N2; i>=2; i--) { W = IV[i] ; // 10進5桁を超えたものは、一つ前の W1 = floor(W*D10) ; // 要素に桁上げ処理をする。 IV[i-1] = IV[i-1] + W1 ; IV[i] = W - W1*P10 ; } // Second Loop for(i=N2; i>=2; i--) // 再度正確な桁揃え  { if(IV[i] >= P10) // 10**5以上の要素の処理 { IV[i] = IV[i] - P10 ; IV[i-1] = IV[i-1] + 1 ; } if(IV[i] < 0) // 負要素の処理 { IV[i] = IV[i] + P10 ; IV[i-1] = IV[i-1] - 1 ;} } } //===================================================================C // Main Program C //===================================================================C void main() { int i, k, ID, NN, N, IX, KETA, T1, T2, T3, T4, KEND; double TT3, TT4; char FNAME[] = "AtanOrg?.txt" ; char FID[] = "12345" ; FILE *FT1; // Constant Set L10 = L10V ; // L10 = 5 P10 = pow(10.0, L10) ; // P10 = 10**5 D10 = 1.0/P10 ; // D10 = 1/10**5 // Parameter Read printf(" Type In ID, (ID=1,2,3,4,5) \n") ; printf(" 1;Machin, 2;Gauss, 3;Klingen., 4;Stoemer, 5;Takano \n") ; scanf("%d",&ID) ; // 計算公式選択 if(ID <= 0 || ID >= 6) ID = 1 ; printf(" Type In N (N=Number of Pi Decimal-Digits) \n") ; scanf("%d",&NN) ; // 計算桁数入力  N = NN/L10 ; if(N >= ND-3) N = ND - 3 ; NN = N*L10 ; N2 = N + 3 ; // 計算要素数を求める(誤差を考慮して) // Open Output File ID = ID - 1; FNAME[7] = FID[ID] ; FT1 = fopen(FNAME,"w") ; printf(" ** Pi Computation Start ** \n") ; printf(" Number of Pi Decimal-Digits= %d \n",NN) ; fprintf(FT1," ** Pi Computation with the original arctan method ** \n") ; fprintf(FT1," Number of Pi Decimal-Digits= %d, ID= %d \n",NN,ID+1) ; T1 = clock(); // Clear Pi Area for(i=1; i <= N2; i++) IV[i] = 0 ; // Pi Computation Start for(k=0; k < INO[ID]; k++) // Atanの個数だけ反復計算 { T2 = clock() ; M = IDATA[ID][k][0] ; // Set M X = IDATA[ID][k][1] ; // Set X ATANV() ; // IV = IV + M*Atan(1/X) ---> π T3 = clock() ; TT3 = (double)(T3 - T2)/CLOCKS_PER_SEC ; IX = X ; printf(" Atan(1/ %d ) Computation, Time(s)= %f \n",IX,TT3) ; } // Output NORMAL() ; // 各要素を10進5桁の値にそろえる T4 = clock() ; TT4 = (double)(T4-T1)/CLOCKS_PER_SEC ; // πの計算時間 printf(" Pi Computation, Time(s)= %f \n",TT4) ; fprintf(FT1," Pi Computation, Time(s)= %f \n\n",TT4) ; fprintf(FT1," Pi = 3.") ; for (k=2; k<=N2-2; k=k+16) // Output Pi Value { KEND = k+15 ; // 80 Decimal in a Line if(KEND > N2-2) KEND = N2-2 ; KETA = (k-2)*L10 + 1 ; fprintf(FT1,"\n %8d",KETA) ; for (i=k; i<=KEND; i=i+2) { fprintf(FT1," %5.5d%5.5d",IV[i],IV[i+1]) ; } } }