円周率1万桁の計算（ガウス=ルジャンドル法）

このプログラムは、 ガウス=ルジャンドル法 を用いて、円周率を1万桁まで計算するVCSSLプログラムです。

使用方法

ダウンロードと展開（解凍）

まず、PC（スマホは未対応）で上の画面の「 ダウンロード 」ボタンを押してください。 するとZIP形式で圧縮されたファイルがダウンロードされるので、そのZIPファイルを右クリックして「すべて展開」や「ここに展開」などで展開（解凍）してください。 展開が成功すると、ZIPファイルと同じ名前のフォルダができ、その中にZIPファイルの中身が入っています。

» 展開がエラーで止まってしまう場合は…

なお、Linux® 等をご使用で、右クリックメニューから展開するとファイル名が文字化けしてしまう場合は、 コマンドライン端末でZIPファイルのある場所まで cd した上で「 unzip -O cp932 ZIPファイル名 」で展開してみてください。

プログラムの起動

Microsoft® Windows® をご使用の場合

上記の通りにZIPファイルを展開したフォルダ内にある、 「 VCSSL.bat（種類はバッチファイル） 」をダブルクリック実行してください。 もしプログラムの内容を書き変えながら使いたい場合は、代わりに「 VCSSL_Editor.bat 」を実行してください。

実行すると、最初にメモリー使用量や、（必要な場合のみ）Java®実行環境を自動で入手するか※ 等を尋ねられるので、適時答えると、プログラムが起動します。2回目以降はすぐに起動します。

※ ここで入手したJava®実行環境は、ZIPファイルを展開した中の「 jre 」フォルダ内にダウンロードされ、このプログラムの実行のみに使用されます。PC全体に影響する形でインストールされる事はありません。

Linux® 等やその他のOSをご使用の場合

ZIPファイルを展開したフォルダ内へコマンドライン端末で cd して、以下の通り入力して実行してください：

» javaコマンドが使用できない等のエラーが表示される場合は…

操作方法

プログラムを実行すると、負荷に関する注意メッセージが表示され、実際に円周率を計算するかどうかを尋ねられます。そこで「はい」を選択すると、円周率の計算がスタートします。

計算時間は、最近のパソコンであれば概ね10秒を切りますが、古いパソコンや、低性能の省電力パソコンなどでは、より長い時間を要します。 計算途中で強制終了したい場合はウィンドウを閉じてください。

計算が完了すると、画面に円周率1万桁と、計算に要した時間が表示されます。

最後に、計算結果をファイルに保存するか尋ねられます。「はい」を選択した場合は「pi.txt」に結果が保存されます。

題材解説

ガウス=ルジャンドル法

ガウス=ルジャンドル法 は、円周率を計算するための強力なアルゴリズムです。 同じ変数に対して計算を繰り返す事によって、真の値に収束させていくタイプのアルゴリズムで、特筆すべきはその収束の速さです。 具体的には、ループを一周するごとに、値の精度（正しい桁数）が毎回倍増していくという、驚異的な収束特性を誇ります。 つまり、ループ回数に対して、精度が指数関数的に向上します。

円周率を求めるアルゴリズムとしては、他にも arctan 関数や arcsin 関数のマクローリン展開（0近傍でのテイラー展開） を用いたものが有名です。 しかしガウス=ルジャンドル法は、それらに比べて遥かに高速である上に、アルゴリズムの複雑さも大差ありません。

ガウス=ルジャンドル法による円周率計算の処理

ガウス=ルジャンドル法では、以下の漸化式を繰り返し処理する事によって、円周率計算に必要な値の組を収束させ、最後にそれらの値から円周率を求めます。

\[ a_{n+1} = \frac{a_n + b_n}{2} \ \ … (1) \]
\[ b_{n+1} = \sqrt{a_n b_n} \ \ … (2) \]
\[ c_{n+1} = c_n - d_n(a_n-a_{n+1})^2 \ \ … (3) \]
\[ d_{n+1} = 2 d_n \ \ … (4) \]

なお、各変数の初期値は以下のように指定します。

\[ a_{0} = 1 \ \ … (5) \]
\[ b_{0} = \frac{1}{\sqrt{2}} \ \ … (6) \]
\[ c_{0} = \frac{1}{4} \ \ … (7) \]
\[ d_{0} = 1 \ \ … (8) \]

これらの\( a_n \) 〜 \( d_n \)の（\(a_nとb_nが収束した時点における\)）収束値の組を用いて、円周率\( \pi \)の近似値は以下のように求まります。

\[ \pi = \frac{ ( a_n + b_n )^2 }{ 4 c } \ \ … (9) \]

こうして得られる\( \pi \)における信頼できる桁数は、\( 2^n \) のオーダーで極めて俊敏に向上します。

ニュートン法による平方根計算の処理

さて、上で述べたガウス=ルジャンドル法の処理は極めて単純ですが、一つ飛ばしている点があります。それは(2)式と(6)式で平方根を使用している事についてです。 つまり、ガウス=ルジャンドル法によって円周率を計算するためには、まず平方根を計算する処理を用意する必要があります。

このプログラムでは、平方根の計算にはニュートン法を採用しています。 ニュートン法もまた、繰り返し計算によって値を収束させていくタイプのアルゴリズムです。ここではニュートン法そのものについて詳しく述べる事はせず、平方根を求める処理だけを簡潔に解説します。

まず、変数\(x\)が、与えられた\(a\) の平方根 \(\sqrt{a}\) の値に等しい時に \( f(x) = 0 \) となるような、関数\( f(x) \)を定義します。具体的な定義は以下のようになります。

\[ f(x) = x^2 - a \]

この\( f(x) \)は、確かに\( x = √a \)の時に\( 0 \)となります。この\( f(x) \)を基にして、以下のような漸化式を作ります。

\[ x_{n+1} = x_{n} - \frac{f(x)}{f'(x)} \]

これがニュートン法の一般的な漸化式で、\( x_{n} \)の収束値として\(f(x_{n}) = 0\)を満たす値が得られます（収束するためにいくつかの条件はありますが）。ここに上で定義した\( f(x) \)を具体的に代入して：

\[ x_{n+1} = x_{n} - \frac{ x_{n}^2 - a }{ 2 x_{n} } \]

このままでも良いのですが、少しでも計算回数を節約したいので、もう少しまとめると：

\[ x_{n+1} = \frac{1}{2} \{ x_n + \frac{a}{x_n} \} \ \ … (10) \]

という漸化式が得られます。

(10)の漸化式は、繰り返し処理する度に、\(x_{n+1}\)の値が\(\sqrt{a}\)の値に近づいて行きます。つまり適当な \( x_0 \) からスタートして、\(x_{n+1}\)と\(x_n\)の値がほぼ同じ値に収束するまで回すと、\( x_n \)に\( \sqrt{a} \)の値が得られます。

コード解説

ここからは、実際に円周率を求めるコード内容について解説します。

まずは、コード全体をそのまま掲載します。

以上の通りです。それでは、細部を見ていきましょう。

ライブラリの読み込み

まずは、プログラムの先頭領域で、必要なライブラリを読み込んでいます。

円周率を計算する肝心な部分に関しては、ゼロから組んでいるため、特別なライブラリは必要ありません。 しかし、計算時間を計測したり、結果を読みやすい形に加工して表示したりするために、標準ライブラリである「 Time 」ライブラリと「 Text 」ライブラリを import で読み込んでいます。

ニュートン法による平方根計算の実装

続いて、後のガウス=ルジャンドル法の実装で必要になる、平方根を求める処理を実装しています。

具体的には、(10)式をほぼそのまま、収束するまで繰り返しています。 収束の判定には、一回前のループ時点での値を変数 valueStock に保持しておき、1回処理後に値が変わっていなければ、収束したとみなす方法を用いています（末尾桁が振動するようなケースでは、ある程度絶対値が接近した時点で収束と見なす必要がありますが、今回の場合は必要ありません）。

なお、2個目の引数 start には、ニュートン法の初期値を指定するようにしています。ニュートン法では、初期値が収束値に近いほど、高速な動作が期待できます。

今回は、ガウス=ルジャンドル法の収束部で平方根を求めるので、非常に近い平方根の値を、収束ループ毎に何度も求める事になります。そこで、ガウス=ルジャンドル法の前回ループ時の結果を保持しておき、次回ループ時の初期値に与える事で、平方根の収束を速めています。

ガウス=ルジャンドル法による円周率計算の実装

続いていよいよ、ガウス=ルジャンドル法で円周率を求める部分の実装です。

これは素直に、初期値を設定して、(1)〜(4)式を繰り返し処理しています。平方根と同様、変数 a と b の値を、変数 aStock と bStock に保持しておき、収束ループを一周しても値が変わっていなければ、収束したと見なしています（これも場合によっては、絶対値がある程度接近した時点で収束と見なす必要がありますが、今回は必要ありません）。

引数 printProgress には、計算の進行状況を表示するかどうかを指定できるようにしています。true の場合は、収束ループが一周するごとにタイムを表示します。

その他、プログラムの制御的な部分など

円周率の計算に関する実装は以上です。これだけでも、pi 関数をコールすると、円周率が計算されて返ってきます。

実際のプログラムには、この他にも pi 関数をコールして、結果を見やすい形にして表示するような、プログラムの制御的な部分もあります。しかしそれらは本題材にあまり関係が無いので、詳しい解説は省略します。 おおまかな流れとしては、まず process 関数内でvarfloat型（多倍長精度の浮動小数点を扱う型）の演算桁数を設定して、続いて円周率を計算し、最後に円周率を getOutputText 関数内で見やすい形にまとめ、画面に表示して終了といった具合です。

ライセンス

このVCSSLコード（ 拡張子が「.vcssl」のファイル ）は実質的な著作権フリー（パブリックドメイン) である CC0 の状態で公開しています。 そのままでのご利用はもちろん、言語の種類を問わず、改造や流用などもご自由に行ってください。

※ ただし、このVCSSLコードの配布フォルダ内には、ダウンロード後すぐに実行できるように、 VCSSLの実行環境も同梱されており、そのライセンス文書は「 License 」フォルダ内に同梱されています （要約すると、商用・非商用問わず自由に使用できますが、使用の結果に対して開発元は一切の責任を負いません、といった具合の内容です）。 配布フォルダ内の各構成物の一覧やライセンスについては「 ReadMe_使用方法_必ずお読みください.txt 」をご参照ください。

この記事中の商標などについて

• OracleとJavaは、Oracle Corporation 及びその子会社、関連会社の米国及びその他の国における登録商標です。文中の社名、商品名等は各社の商標または登録商標である場合があります。
• Windows は、米国 Microsoft Corporation の米国およびその他の国における登録商標です。この記事は独立著作物であり、Microsoft Corporation と関連のある、もしくはスポンサーを受けるものではありません。
• Linux は、Linus Torvalds 氏の米国およびその他の国における商標または登録商標です。
• その他、文中に使用されている商標は、その商標を保持する各社の各国における商標または登録商標です。