数値計算プログラム課題

数値計算プログラム課題：円周率πを計算する

マチンの公式
練習課題その５：マチンの公式によりπを求めよ
数値計算の精度を上げるには
より高精度にπを求めよ

1. マチンの公式 　計算機でπを求める時に良く使われるのが、マチン(Machin)の公式です。

マチンの公式

　　π　＝　16 × arctan(1/5) − 4 × arctan(1/239)　　

　また、arctan()はテーラー展開により、次式で求められます。

arctan(x)

　　arctan(x)　＝　x　−　(x^3)/3　＋　(x^5)/5　−　(x^7)/7　…　　（但し |x|<=1）　　
　　　＝　Σ_{n=1}^{∞} (-1)^(n+1) (x^(2n-1)) / (2n-1)　　

2. 練習課題その５：マチンの公式によりπを求めよ

　次のプログラムは、マチンの公式によりπを求めるプログラムの一部です。このプログラムを参考にしてπを求めるプログラムを完成させましょう。

　少し難しくなってきました。このプログラムについて、注意すべき点をいくつか上げます。

1. forループで、arctan()の第n項を1から30まで変数PIに足していくと良いでしょう。

       for ( n=1; n<=30; n++) は、次のようにして{}の中を30回ループします。

          (A)  n = 1:    ループの前にnに1を代入する 
          (B)  n<=30:    n<=30の間だけループする
          (C)   n++ :    毎ループ終了時にnを1増やす

2. CalcATAN_Term()は、arctan(x)のn番目の項を計算する関数です。

3. int main(void)の前にある"double CalcATAN_Term(int n, double x);"は関数のプロトタイプ宣言と呼ばれるものです。
main関数の後ろでmain関数の中で使われている関数が定義されている場合に必要になります。
具体的には、コンパイラがプログラムをコンパイルする場合、最初の行から順に解釈していきます。main関数よりも後ろにmain関数の中で使う関数が定義されていた場合、コンパイラは「こんな関数は定義されていない」と怒ります。
そこで、あらかじめ、定義する関数の名前や関数の型、引数の型のみをmain関数よりも前に書いておくことにより、コンパイラはコンパイル時にまず「このプログラムではこういう関数が使われるんだ」と覚えてくれます。それにより、関数がどこに定義されていても怒らなくなります。

4. "printf("PI[%d] = %1.111lf\n", n, PI);"の中の"%1.111lf"は、桁数を指定して浮動小数点の数値を表示しています。この例では、整数部1桁、小数部111桁で表示されます。

★"%1.10lf"に変えて表示の変化を確認せよ

5. 関数pow(x, y)はべき関数(Power)で、xのy乗をdouble型で返します。x、yはdouble型で与えます。

/* 課題 --- πの数値計算プログラム */ #include <stdio.h> #include <math.h> /* 関数CalcATAN_Term()のプロトタイプ宣言 */ double CalcATAN_Term(int n, double x); int main(void) { int n; double PI; PI = 0.0; /* * n=30 までΣを計算する * n=1,2,…について、arctan()の項をPIに加算していく */ for( n = 1; n <= 30; n++) { （ここに入るプログラムを考えよう） printf("PI[%d] = %1.111lf\n", n, PI); } return 0; } double CalcATAN_Term(int n, double x) { double val; val = pow(x, (double)(2*n-1)) / (double)(2*n - 1); if ( (n%2==0) ) val = -val; return val; }

　プログラムを打ち込み終ったら、コンパイルし、実行して下さい。
結果は次のようになります。　注: 式の形によっては全く同じにならないことがあります。

PI[1] = 3.183263598326360188650596683146432042121887207031250000000000000000000000000000000000000000000000000000000000000 PI[2] = 3.140597029326060773968265493749640882015228271484375000000000000000000000000000000000000000000000000000000000000 PI[3] = 3.141621029325035063806126345298252999782562255859375000000000000000000000000000000000000000000000000000000000000 PI[4] = 3.141591772182177777494871406815946102142333984375000000000000000000000000000000000000000000000000000000000000000 PI[5] = 3.141592682404399816675777401542291045188903808593750000000000000000000000000000000000000000000000000000000000000 PI[6] = 3.141592652615309066987947517191059887409210205078125000000000000000000000000000000000000000000000000000000000000 PI[7] = 3.141592653623555442266024329001083970069885253906250000000000000000000000000000000000000000000000000000000000000 PI[8] = 3.141592653588602956915565300732851028442382812500000000000000000000000000000000000000000000000000000000000000000 PI[9] = 3.141592653589836636740528774680569767951965332031250000000000000000000000000000000000000000000000000000000000000 PI[10] = 3.141592653589792671908753618481568992137908935546875000000000000000000000000000000000000000000000000000000000000 PI[11] = 3.141592653589794448265593018732033669948577880859375000000000000000000000000000000000000000000000000000000000000 PI[12] = 3.141592653589794448265593018732033669948577880859375000000000000000000000000000000000000000000000000000000000000 PI[13] = 3.141592653589794448265593018732033669948577880859375000000000000000000000000000000000000000000000000000000000000 PI[14] = 3.141592653589794448265593018732033669948577880859375000000000000000000000000000000000000000000000000000000000000 PI[15] = 3.141592653589794448265593018732033669948577880859375000000000000000000000000000000000000000000000000000000000000 PI[16] = 3.141592653589794448265593018732033669948577880859375000000000000000000000000000000000000000000000000000000000000 PI[17] = 3.141592653589794448265593018732033669948577880859375000000000000000000000000000000000000000000000000000000000000 PI[18] = 3.141592653589794448265593018732033669948577880859375000000000000000000000000000000000000000000000000000000000000 PI[19] = 3.141592653589794448265593018732033669948577880859375000000000000000000000000000000000000000000000000000000000000 PI[20] = 3.141592653589794448265593018732033669948577880859375000000000000000000000000000000000000000000000000000000000000 PI[21] = 3.141592653589794448265593018732033669948577880859375000000000000000000000000000000000000000000000000000000000000 PI[22] = 3.141592653589794448265593018732033669948577880859375000000000000000000000000000000000000000000000000000000000000 PI[23] = 3.141592653589794448265593018732033669948577880859375000000000000000000000000000000000000000000000000000000000000 PI[24] = 3.141592653589794448265593018732033669948577880859375000000000000000000000000000000000000000000000000000000000000 PI[25] = 3.141592653589794448265593018732033669948577880859375000000000000000000000000000000000000000000000000000000000000 PI[26] = 3.141592653589794448265593018732033669948577880859375000000000000000000000000000000000000000000000000000000000000 PI[27] = 3.141592653589794448265593018732033669948577880859375000000000000000000000000000000000000000000000000000000000000 PI[28] = 3.141592653589794448265593018732033669948577880859375000000000000000000000000000000000000000000000000000000000000 PI[29] = 3.141592653589794448265593018732033669948577880859375000000000000000000000000000000000000000000000000000000000000

　実行結果を見ると、nが増えるにつれ、PIは正しい値に収束していくことが分かります。
しかし、n=11以降ではPIはそれ以上変化しません。ある桁以下は0が続いているので、正しい円周率に収束したわけではではありません。

　これは、計算に使用したdouble型の精度が原因です。n=11以降はdouble型で表現できるよりも小さな幅で収束していきます。このため、計算結果の丸めにより、PIはある値から変化しないのです。

円周率計算プログラムの解答例

/* 課題7-1 "exercise7-1.c" --- πの数値計算プログラム */ /* マチンの公式で円周率を求める */ #include <stdio.h> #include <math.h> /* 関数CalcATAN_Term()のプロトタイプ宣言 */ double CalcATAN_Term(int n, double x); int main(void) { int n; double PI; PI = 0.0; /* * arctan(x)をテーラー展開した第1項から第30項まで足していく * 第n項はCalcATAN_Term(n,x)で表される。 * for文で n=30 まで足しあわせる。 * PIはマチンの公式であたえる。 */ for(n=1; n<=30; n++) { PI = PI + (16.0 * CalcATAN_Term(n, 1.0/5.0) - 4.0 * CalcATAN_Term(n, 1.0/239.0)); printf("PI[%d] = %1.111lf\n", n, PI); } return 0; } double CalcATAN_Term(int n, double x) /* arctan(x)をテーラー展開した第n項を返す */ { double val; val = pow(x, (double)(2*n-1)) / (double)(2*n - 1); if ( (n%2==0) ) val = -val; return val; }

3. 数値計算の精度を上げるには

　より精度の高い型の変数を使用することにより、より正しい円周率を求めることができます。
実習で使用するシステムでは、long doubleという型により倍の精度での計算を行なうことができます。

　プログラム中の"double"を全て"long double"に置き換えましょう。
ただし、printf文の中の"%1.111lf"を"%1.111Lf"に変える必要があります。
すると、実行結果は以下のようになり、より正確な円周率が求まります。

PI[1] = 3.183263598326360011708802133512108412105590105056762695312500000000000000000000000000000000000000000000000000000 PI[2] = 3.140597029326060483237282419834287758297897141801513498649001121520996093750000000000000000000000000000000000000 PI[3] = 3.141621029325034594236822505700749617445715368673865545726698594113080430234585804782909690402448177337646484375 PI[4] = 3.141591772182177464197909768805470665390451975891436757351950060079262203038830136847536778077483177185058593750 PI[5] = 3.141592682404399686420390462942965624972967358560440429222170886989877469641641027919831685721874237060546875000 PI[6] = 3.141592652615308777329463926117373603725933150205709520895824509483982255042544551315586431883275508880615234375 PI[7] = 3.141592653623554931175618439188173690342857441668227947223737950481026008631157964146041194908320903778076171875 PI[8] = 3.141592653588602397842285075566152872293070827983711315150961999507171762413548776748939417302608489990234375000 PI[9] = 3.141592653589836016665814488559529373726337088092911901387159060480402245074671441216196399182081222534179687500 PI[10] = 3.141592653589791866097393435881213288321901066944055130665093259108040096849734368333884049206972122192382812500 PI[11] = 3.141592653589793463927488673978407108119363493676226541070115918431629280227479483755814726464450359344482421875 PI[12] = 3.141592653589793405571954760934853968181392388764760576338692105581411329051633174458402208983898162841796875000 PI[13] = 3.141592653589793407719438408934857132336406726621461108197088800712479250698550004017306491732597351074218750000 PI[14] = 3.141592653589793407639901977527449778931845696677362841622594930664127271802144036882964428514242172241210937500 PI[15] = 3.141592653589793407642864023938484125793989076941919612900143935093049546436461127996153663843870162963867187500 PI[16] = 3.141592653589793407642753186072780816107157893599060663518131979391324708206223448314631241373717784881591796875 PI[17] = 3.141592653589793407642757350889552502419112350609077804536181412130529770621034657551717828027904033660888671875 PI[18] = 3.141592653589793407642757193816462842799929623328301125379895659918069728000489249097881838679313659667968750000 PI[19] = 3.141592653589793407642757199759769021312220618943718753498368854214857666395932511704813805408775806427001953125 PI[20] = 3.141592653589793407642757199534227987756897957558380510989310219766425073817384117091933148913085460662841796875 PI[21] = 3.141592653589793407642757199542809546403098603072427337505880138912367147141502243812283268198370933532714843750 PI[22] = 3.141592653589793407642757199542482138312552772977407821255640701268215669861305627819092478603124618530273437500 PI[23] = 3.141592653589793407642757199542494464264196851286867379514473244803054313711854206303542014211416244506835937500 PI[24] = 3.141592653589793407642757199542494079078207973839696768318884727817590606091524563225902966223657131195068359375 PI[25] = 3.141592653589793407642757199542494079078207973839696768318884727817590606091524563225902966223657131195068359375 PI[26] = 3.141592653589793407642757199542494079078207973839696768318884727817590606091524563225902966223657131195068359375 PI[27] = 3.141592653589793407642757199542494079078207973839696768318884727817590606091524563225902966223657131195068359375 PI[28] = 3.141592653589793407642757199542494079078207973839696768318884727817590606091524563225902966223657131195068359375 PI[29] = 3.141592653589793407642757199542494079078207973839696768318884727817590606091524563225902966223657131195068359375

4. アドバンスト課題：より高精度にπを求めよ

　上記の方法では精度が上がるといっても、n=24以降はこれ以上収束していません。では、円周率を何千万桁も
求めるにはどうしたら良いでしょうか。ここまで早く終った人は、考えてみて下さい。
（ヒント：「多倍長演算」）
この課題ができた人は、担当教員宛にソースファイルを送ってください。（学籍番号、氏名を記入のこと）

戻る