多数桁π計算の原理プログラム
2002/12/20 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
2006/09/24 Windows版実行プログラム及び実行結果を追加
1. 概要
多数桁のπ(円周率)の値を計算するFORTRAN及びCのソースプログラムを掲載する。
本プログラムは内容を理解してもらうことを目的とし、計算時間及び使用メモリ量の効率化
は目的としていない。本格的なプログラムは金田研究室から他のプログラムと合わせて
公開する予定。
2. 10進多数桁π計算の原理プログラム
(1) 従来のArctan公式計算プログラム
(a) 計算公式 ： πを求める公式の2.1 Arctan公式 | を参照。
(b) 計算方法 ： 4バイト整数に10進5桁を入れて計算。演算量は桁数の２乗に比例する。
(c) 実行方法 ： 使用公式(
2.1 Arctan公式 | )と桁数を画面上から応答する。
画面上に利用したArctanとπの計算時間を表示する。πの値はAtanOrg?.txtに
出力する。?は1〜5の利用公式の番号である。
(d) プログラム ： FORTRAN |
とC | のソースプログラム
Cソースプログラム： orgatan.c
Windows版実行プログラム： orgatan.exe
実行結果(πの10進20万桁)： orgatan.txt
(2) AGM法(算術幾何平均、相加相乗平均）での計算プログラム
(a) 計算公式 ： πを求める公式の3.1 ガウス・ルジャンドル公式 | を参照。
(b) 計算方法 ： 倍精度実数に10進5桁を入れて計算。演算量は桁数nのn･log2に比例する。
多数桁乗算は複素数のFMT(高速剰余法)を使用している。
(c) 実行方法 ： 桁数を画面上から応答する。計算桁数は2のべき乗の5倍である。
現在の利用公式はウス・ルジャンドル公式であるが、利用公式の応答要。
画面上にAGMの各stepとπの計算時間を表示する。πの値はAGM-2nd?.txtに出力
する。πの計算時間と共にRounding Errorが出力される。これは、FMTによる
多数桁乗算時に丸め誤差のため実数を整数にする過程で発生した最大誤差
である。これが0.1以上の値となる場合は１要素に詰める桁数を10進5桁から
4桁にする必要がある。
(d) プログラム ： FORTRAN |
とC | のソースプログラム
Cソースプログラム： agm.c, pibase.c,
mult.c, add.c
Windows版実行プログラム： agm.exe
実行結果(πの10進20万桁)： agm.txt
(3) DRM法によるAtan公式計算プログラム
(a) 計算公式 ： πを求める公式の2.1 Arctan公式 | を参照。
(b) 計算方法 ： 倍精度実数に10進5桁を入れて計算。演算量は桁数nのn･log2〜3に比例する。
多数桁乗算は複素数のFMT(高速剰余法)を使用している。
(c) 実行方法 ： 使用公式(
2.1 Arctan公式 | )と桁数を画面上から応答する。
画面上に利用したArctanとπの計算時間を表示する。πの値はAtanDRM?.txtに
出力する。?は1〜5の利用公式の番号である。
πの計算時間と共にRounding Errorが出力される。これは、FMTによる多数桁
乗算時に丸め誤差のため実数を整数にする過程で発生した最大誤差である。
これが0.1以上の値となる場合は１要素に詰める桁数を10進5桁から4桁にする
必要がある。
(d) プログラム ： FORTRAN |
とC | のソースプログラム
Cソースプログラム： drmatan.c, pibase.c,
mult.c, add.c
Windows版実行プログラム： drmatan.exe
実行結果(πの10進20万桁)： drmatan.txt
(4) DRM法によるラマヌジャン公式計算プログラム
(a) 計算公式 ： πを求める公式の2.2 ラマヌジャン型公式 | を参照。
(b) 計算方法 ： 倍精度実数に10進5桁を入れて計算。演算量は桁数nのn･log2〜3に比例する。
多数桁乗算は複素数のFMT(高速剰余法)を使用している。
(c) 実行方法 ： 使用公式(
2.2 ラマヌジャン型公式 | )と桁数を画面上から応答する。
画面上に利用したDRMの各step単位とπの計算時間を表示する。πの値は
RamnDRM?.txtに出力する。?は1〜3の利用公式の番号である。
πの計算時間と共にRounding Errorが出力される。これは、FMTによる多数桁
乗算時に丸め誤差のため実数を整数にする過程で発生した最大誤差である。
これが0.1以上の値となる場合は１要素に詰める桁数を10進5桁から4桁にする
必要がある。
(d) プログラム ： FORTRAN |
とC | のソースプログラム
Cソースプログラム： drmrmnj.c, pibase.c,
mult.c, add.c
Windows版実行プログラム： drmrmnj.exe
実行結果(πの10進20万桁)： drmrmnj.txt
(5) DRM法による連分数公式計算プログラム
(a) 計算公式 ： πを求める公式の4. 連分数公式 | を参照。
(b) 計算方法 ： 倍精度実数に10進5桁を入れて計算。演算量は桁数nのn･log2〜3に比例する。
多数桁乗算は複素数のFMT(高速剰余法)を使用している。
(c) 実行方法 ：
(d) プログラム ： FORTRANとCのソースプログラム
3. 任意基底多数桁π計算の原理プログラム
DRM法によるAtan公式計算プログラムとマヌジャン公式計算プログラムを載せる。
3〜62進数までの任意の基底による円周率の多数桁計算を行う。
任意基底への変換は小数点以下に対してで、最初の値は常に3と表示する。
結果の表示方法は下記の2種類がある。
(a) 表示方法に'0'と応答
0〜 9 ： そのまま0〜9の文字で表示
10〜35 ： a=10,b=11,c=12,...,y=34,z=35と英小文字で表示
36〜61 ： A=36,B=37,C=38,...,Y=60,Z=61と英大文字で表示
(b) 表示方法に'1'と応答
英小文字(26進数)や英大小文字(52進数)で表示するときに都合が良い。
0〜25 ： a=0,b=1,c=2,...,y=24,z=25と英小文字で表示
26〜51 ： A=26,B=27,C=28,...,Y=50,Z=51と英大文字で表示
52〜61 ： 表示不可
(1) DRM法によるAtan公式計算プログラム
(a) 計算公式 ： πを求める公式の2.1 Arctan公式 | を参照。
(b) 計算方法 ： 倍精度実数に指定した基底数(3〜62)を3〜11桁入れて計算。倍精度実数に
入れる桁数は自動調整する。演算量は桁数nのn･log2〜3に比例する。
多数桁乗算は複素数のFMT(高速剰余法)を使用している。
(c) 実行方法 ： 最初に、使用する基底数(3〜62)と表示方法(0 or 1)を応答する。
次に、使用公式(
2.1 Arctan公式 | )と桁数を画面上から応答する。
画面上に利用したArctanとπの計算時間を表示する。πの値はAtanDRM?.txtに
出力する。?は1〜5の利用公式の番号である。
(d) プログラム ： FORTRAN |
とC | のソースプログラム
Cソースプログラム： xdrmatan.c, xpibase.c,
mult.c, add.c
Windows版実行プログラム： xdrmatan.exe
実行結果(πの16進(0〜e)20万桁)： xdrmatan.txt
(2) DRM法によるラマヌジャン公式計算プログラム
(a) 計算公式 ： πを求める公式の2.2 ラマヌジャン型公式 | を参照。
(b) 計算方法 ： 倍精度実数に指定した基底数(3〜62)を3〜11桁入れて計算。倍精度実数に
入れる桁数は自動調整する。演算量は桁数nのn･log2〜3に比例する。
多数桁乗算は複素数のFMT(高速剰余法)を使用している。
(c) 実行方法 ： 最初に、使用する基底数(3〜62)と表示方法(0 or 1)を応答する。
次に、使用公式(
2.2 ラマヌジャン型公式 | )と桁数を画面上から応答する。
画面上に利用したDRMの各step単位とπの計算時間を表示する。πの値は
RamnDRM?.txtに出力する。?は1〜3の利用公式の番号である。
(d) プログラム ： FORTRAN |
とC | のソースプログラム
Cソースプログラム： xdrmrmnj.c, xpibase.c,
mult.c, add.c
Windows版実行プログラム： xdrmrmn.exe
実行結果(πの26進(a〜z)20万桁)： xdrmrmnj.txt