5.プログラムの例
・具体的なプログラム示します。Webページ上では空白が消されてインデントがバラバラになってます。お使いの際は、CalcPi.cをダウンロードしてください。
・このプログラムには計算時間がわかるように、_dos_gettime()
を組み込みました。
・計算桁数は１０００桁としてありますが、３０００桁までは大丈夫です。
もちろんコンパイルオプションには注意してください。スタック不足ではコンパイルできません。
・このプログラムはあまり工夫のないものです。それだからこそ、工夫のしがいがあると思います。
・１０万桁、１００万桁の計算は、最初に述べた通りDOSでは無理なため、Windows用プログラムとし、Delphi4
を使って開発しました。プログラムは載せませんでしたが、実行ファイルのダウンロードはできるようにしておきました。自己解凍型の圧縮ファイル
PI.EXE （226KB）にしてあります。
ただ単に円周率を計算するだけ、しかも遅いですから、何の役にもたちませんが・・・(^^;
/**************************************************************************
* Stormer の公式による円周率の計算
Copyright (C) 1998-1999 S.H.
powebee@pluto.dti.ne.jp
**************************************************************************/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <dos.h>
#define BASE 10000
#define TARGET 1000
#define B8 1.80618 /* 2*log10(8 ) */
#define B57 3.51175 /* 2*log10(57) */
#define B239 4.75680 /* 2*log10(239) */
#define BU8 64 /* 8^2 */
#define BU57 3249 /* 57^2 */
#define BU239 57121 /* 239^2 */
/**************************************************************************/
/* Global変数 */
/**************************************************************************/
unsigned long gMaxarray; /* 配列の要素数
*/
unsigned long gNofTerm; /* 項数 */
unsigned long gAc[1100], gDc[1100];
long gPc[1000];
struct dostime_t s,e;
/**************************************************************************/
/* 関数の宣言 */
/**************************************************************************/
void compute(unsigned long nof, unsigned
long bun);
void divide_b(unsigned long bun);
void divide_add(unsigned long nodd);
void divide_sub(unsigned long nodd);
void outprint(void);
main()
{
long i;
/* 10,000進法ゆえ，1配列要素には 9999までの数が収まる。
*/
/* よって，n桁を求める時は，その4分の1の配列要素があればよい
*/
/* 2を加算するのは，余裕を持たせるため */
gMaxarray = TARGET/4 + 2;
/* 初期化 i=0 の要素は使用しない */
for (i=1; i<=gMaxarray; i++)
*(gAc+i) = *(gDc+i) = *(gPc+i) = 0;
/* 計算開始時間の取得 */
_dos_gettime(&s);
/* 第1級数 */
*(gAc + 1) = 192;
gNofTerm = TARGET / B8 + 2; /* 計算項数を求める
*/
/* 計算を行う */
compute(gNofTerm, BU8);
/* 第2級数 */
*(gAc + 1) = 456;
gNofTerm = TARGET / B57 + 2;
/* 計算を行う */
compute(gNofTerm, BU57);
/* 第3級数 */
*(gAc + 1) = 956;
gNofTerm = TARGET / B239 + 2;
/* 計算を行う */
compute(gNofTerm, BU239);
/* 計算終了時間の取得 */
_dos_gettime(&e);
/* 出力する */
outprint();
}
/* 計算ルーチン */
/* function :級数を計算する。 */
void compute(unsigned long nof, unsigned
long bun)
{
unsigned long cnttrm;
for(cnttrm=1; cnttrm<=nof; cnttrm++){
/* 分母で除算 */
divide_b(bun);
/* 2n-1 での除算。nに応じて符号が変わる。*/
if ( cnttrm % 2 > 0 )
divide_add(cnttrm);
else
divide_sub(cnttrm);
}
}
/* bunbo で除算する。 */
void divide_b(unsigned long bun)
{
unsigned long r=0, x=0, cntd;
for(cntd=1; cntd<=gMaxarray; cntd++){
x = *(gAc + cntd) + r * BASE;
*(gAc + cntd) = x / bun;
r = x - *(gAc + cntd) * bun;
}
}
/* function : 2n-1 で除算する。 */
/* 引数 : *p :ap,bp,cpのいずれか。 */
/* nodd :項の番号 2n-1 の n のこと。 */
/* これから，符号も分かる。 */
void divide_sub(unsigned long nodd)
{
unsigned long r=0,x=0,cntd,n2;
n2 = 2 * nodd -1 ;
for(cntd=1; cntd<=gMaxarray; cntd++){
x = *(gAc + cntd) + r * BASE;
*(gDc + cntd) = x / n2;
r = x - *(gDc + cntd) * n2;
/* 桁下がり調整 */
*(gPc + cntd) -= *(gDc + cntd);
if( *(gPc + cntd) < 0){
*(gPc + cntd) += BASE;
*(gPc + cntd - 1) -= 1;
}
}
}
void divide_add( unsigned long nodd)
{
unsigned long r=0,x=0,cntd, n2;
n2 = 2 * nodd -1 ;
for(cntd=1; cntd<=gMaxarray; cntd++){
x = *(gAc + cntd) + r * BASE;
*(gDc + cntd) = x / n2;
r = x - *(gDc + cntd) * n2;
/* 桁上がりの調整 */
*(gPc + cntd) += *(gDc + cntd);
if(*(gPc + cntd) >= BASE ){
*(gPc + cntd) -= BASE;
*(gPc + cntd - 1) += 1;
}
}
}
/*==================================================
function : 出力成形処理
結果を格納する配列には、先頭が0のものも
あります。そのまま出力すると、0 は空白に
なるので、これを補正します。
==================================================*/
void outprint()
{
char outpt[TARGET+100];
char sc,hsc;
unsigned long cnt1,cnt2,cnt3;
unsigned long temp1,temp2;
FILE *fp;
/* 出力用配列の初期化*/
for (cnt1=1; cnt1<=TARGET; cnt1++)
*(outpt + cnt1)=0;
/* 10000進法ゆえ、4桁分 */
/* これを１つずつ取り出します */
cnt3 = 4;
for(cnt1=2; cnt1<=gMaxarray; cnt1++){
temp1 = *(gPc + cnt1);
for(cnt2=1;cnt2<=4;cnt2++){
temp2 = temp1 / 10;
*(outpt + (cnt3-cnt2+1)) = temp1 - temp2
* 10;
temp1 = temp2;
}
cnt3 += 4;
}
fp = fopen("PI_TABLE.txt","w");
/* open new file for writing */
/* 桁数の出力 */
fprintf(fp,"n = %d\n",TARGET);
/* 時間の出力 */
sc = e.second-s.second;
hsc = e.hsecond-s.hsecond;
if (e.hsecond > s.hsecond)
hsc = e.hsecond - s.hsecond;
else{
hsc = 100 + e.hsecond - s.hsecond;
sc -= 1;
}
fprintf(fp,"計算開始時間 = %d 時 %d
分 %d 秒 %d\n",s.hour, s.minute, s.second,
s.hsecond);
fprintf(fp,"計算終了時間 = %d 時 %d
分 %d 秒 %d\n",e.hour, e.minute, e.second,
e.hsecond);
fprintf(fp,"計算時間 = %d 秒 %d\n",sc,
hsc);
fprintf(fp,"PI=%1ld.\n",*(gPc +
1));
fprintf(fp," 1 ");
for(cnt1=1; cnt1<=TARGET; cnt1++){
fprintf(fp,"%1d",*(outpt + cnt1));
if(cnt1 % 10 == 0)
if(cnt1 % 50 == 0){
fprintf(fp,"\n");
fprintf(fp,"%10ld ",cnt1+1);
}else fprintf(fp," ");
}
}
