[ 前へ | 目次 | 次へ ]
Now Loading...
ダウンロード
PC (※スマートフォンでは動きません) でダウンロードし、ZIPファイルを右クリックメニューから展開して、できたフォルダ内の「 VCSSL.bat(バッチファイル) 」をダブルクリックすると起動します。 Linux等では「 VCSSL.jar 」をコマンド実行してください。
» 詳しい使用方法や、エラーで展開できない際の対応方法などはこちら

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

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

スポンサーリンク


使用方法

ダウンロードと展開(解凍)

まず、PC(スマホは未対応)で上の画面の「 ダウンロード 」ボタンを押してください。 するとZIP形式で圧縮されたファイルがダウンロードされます。

Microsoft® Windows® をご使用の方は、ここでまずZIPファイルを右クリックし、「プロパティ」を選んで開かれる画面で、 下の方にあるセキュリティ項目の「許可する」にチェックを入れて「OK」で閉じてください。 これを行わないと、ZIP展開やソフト起動時に、警告メッセージが出て展開完了/起動できない場合があります。

その後、ZIPファイルを右クリックして「すべて展開」や「ここに展開」などで展開(解凍)してください。 展開が成功すると、ZIPファイルと同じ名前のフォルダができ、その中にZIPファイルの中身が入っています。

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

プログラムの起動

Microsoft® Windows® をご使用の場合

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

※ もし、ここでセキュリティ警告メッセージが出て起動できない場合は、 ダウンロードと展開 の項目に書かれている手順通りに、ZIPファイルの展開をやり直してみてください。

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

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

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

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

java -jar VCSSL.jar
(プログラムの内容を書き変えながら使いたい場合は、代わりに VCSSL_Editor.jar を実行)

» 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で記述されており、このサイトではVCSSLのプログラミングガイドも無料公開しています。 コードを改造したい方や、新しいコードを書いてみたい方はぜひご活用ください!

» 公式プログラミングガイド一覧

ライセンス

このVCSSL/Vnanoコード( 拡張子が「.vcssl」や「.vnano」のファイル )は実質的な著作権フリー(パブリックドメイン) である CC0 の状態で公開しています。 記事中にC言語/C++/Java言語などでのサンプルコードが掲載されいてる場合は、それらについても同様です。 そのままでのご利用はもちろん、改造や流用などもご自由に行ってください。

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

※ Vnano の実行環境については、別途スクリプトエンジンのソースコードも一般公開しており、 何らかのソフトウェア内に組み込んでご利用いただく事も可能です。詳細はこちらをご参照ください。

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

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


スポンサーリンク



[ 前へ | 目次 | 次へ ]
Vnano版 | 積分値を求めるプログラム (数値積分)

矩形法/台形法/シンプソン法を用いて、積分の値を数値的に求めるコードです。
Vnano版 | 積分値のグラフ描画用データを出力するプログラム

数値的に積分を行い、結果の関数をグラフに描くためのデータを出力するコードです。
入力された数式を積分して値とグラフを表示するツール

画面上で数式を入力すると、それを数値的に積分し、値とグラフを表示してくれるGUIツールです。
シンプソン法による数値積分

積分の値を数値的に求めます。台形法よりも高精度な方法として、被積分関数を微小区間内で二次関数近似して求めた面積を足しあげる、シンプソン法を使用します。
台形法(台形近似)による数値積分

積分の値を数値的に求めます。長方形近似よりも高精度な方法として、台形で近似した微小領域を足しあげる方法を使用します。
矩形法(長方形近似)による数値積分

積分の値を数値的に求めます。長方形の短冊(矩形)で近似した微小領域を足しあげる、最も単純な方法を使用します。
小数(浮動小数点数)から分数へ近似的に変換するツール

小数(浮動小数点数)を、適当な誤差の範囲内で、近い分数に変換してくれるツールプログラムです。
円周率1万桁の計算(ガウス=ルジャンドル法)

ガウス=ルジャンドル法により、円周率を1万桁まで計算するプログラムです。
試し割り法による素数判定ツール

試し割り法を用いて、素数判定を行ってくれる簡易ツールです。
この階層の目次
お知らせ

Vnanoがオープンベータ版に移行、VCSSLの実行環境で標準で実行可能に
2021年04月07日 - ソフト内組み込み用スクリプトエンジン&言語「 Vnano 」がオープンベータ版に移行し、併せて、VCSSLの実行環境でもVnanoコードの実行が可能になりました。詳細をお知らせします。

リニアングラフの最新版をリリース、2D版でも描画エンジンの直接操作が可能に
2021年04月03日 - リニアングラフ2D/3Dの最新版をリリースしました。それぞれのアップデート内容をお知らせします。今回から、2D版でもJava言語APIによる描画エンジンの直接操作が可能になりました。

RINPn のオープンベータ版をリリース! 詳細な公式ガイドも同梱&公開
2021年03月08日 - 2019年より開発進行中のプログラム関数電卓「 RINPn(りんぷん)」が、正式リリースに向けた最終準備段階として、オープンベータ版へと移行しました。その詳細をお知らせします。

新着
[公式ガイドサンプル] 引き算の結果を画面に表示する

「VCSSLスタートアップガイド」内のサンプルコードです。引き算を行って、結果を画面に表示します。
2021年07月08日
[公式ガイドサンプル] 式を複数行にわたって書く

「VCSSLスタートアップガイド」内のサンプルコードです。足し算を行う式を、複数行にわたって記述します。
2021年07月07日
[公式ガイドサンプル] 足し算の結果を画面に表示する

「VCSSLスタートアップガイド」内のサンプルコードです。足し算の結果を求めて、画面に表示します。
2021年07月06日
Vnano版 | ローレンツ方程式を数値的に解くプログラム

ローレンツ方程式を4次ルンゲ=クッタ法によって解き、グラフ描画用のデータを出力するプログラムです。
2021年02月12日
Vnano版 | 積分値のグラフ描画用データを出力するプログラム

数値的に積分を行い、結果の関数をグラフに描くためのデータを出力するコードです。
2020年12月20日
開発元Twitterアカウント