π計算の原理プログラム(atan公式)
2011/04/11 早稲田大学 後 保範 ( Ushiro Yasunori )
1. 概要
多数桁のπ(円周率)の値を計算するCのソースプログラムを掲載する。
Arctan公式を使用し、「定義式」通りと「DRM法」による2種類を掲載する。
計算桁数をnとすると、計算量は「定義式」がO(n2)に、「DRM法」はO(n・log(n)2〜3)
となる。
2. 10進多数桁π計算の原理プログラム
(1) 従来のArctan公式計算プログラム
(a) 計算方法 : 4バイト整数に10進5桁を入れて計算。演算量は桁数の2乗に比例する。
(b) 実行方法 : 使用公式と桁数を画面上から応答する。
画面上に利用したArctanとπの計算時間を表示する。πの値はAtanOrg?.txtに
出力する。?は1〜5の利用公式の番号である。
(c) プログラム : piatan.c
(2) DRM法によるAtan公式計算プログラム
(a) 計算方法 : 倍精度実数に10進5桁を入れて計算。演算量は桁数nのn・log(n)2〜3に比例する。
多数桁乗算は複素数のFMT(高速剰余法)を使用している。
(b) 実行方法 : 使用公式と桁数を画面上から応答する。
画面上に利用したArctanとπの計算時間を表示する。πの値はAtanDRM?.txtに
出力する。?は1〜5の利用公式の番号である。
πの計算時間と共にRounding Errorが出力される。これは、FMTによる多数桁
乗算時に丸め誤差のため実数を整数にする過程で発生した最大誤差である。
これが0.1以上の値となる場合は1要素に詰める桁数を10進5桁から4桁にする
必要がある。
(c) プログラム :drmatan.c, mult.c,
add.c, piout.c
3. Arctan公式によるπ計算
arctan(1/p) = 1/p - 1/(3p3) + 1/(5p5) - 1/(7p7) + …… なる展開式で計算。
多数あるが下記に有名な公式及び日本人の発見した公式を示す。
(1) マチンの公式 : π/4 = 4arctan(1/5) - arctan(1/239)
(2) ガウスの公式 : π/4 = 12arctan(1/18) + 8arctan(1/57) - 5arctan(1/239)
(3) クリンゲンシュティルナの公式:π/4 = 8arctan(1/10) - arctan(1/239) - 4arctan(1/515)
(4) シュテルマーの公式 : π/4 = 44arctan(1/57) + 7arctan(1/239) - 12arctan(1/682)
+ 24arctan(1/12943)
(5) 高野喜久雄の公式: π/4 = 12arctan(1/49) + 32arctan(1/57) - 5arctan(1/239)
+ 12arctan(1/110443)