円周率を１億桁まで求めてみる

円周率ってどうやって求めるのか気になったので実装してみた。
求め方はいろいろあるらしくウィキペディアにあったこのアルゴリズムを利用した。

ガウス＝ルジャンドルのアルゴリズム（英語: Gauss–Legendre algorithm）は、円周率を計算する際に用いられる数学の反復計算アルゴリズムである。円周率を計算するものの中では非常に収束が速く、2009年にこの式を用いて2,576,980,370,000桁（約2兆6000億桁）の計算がなされた。

初期値の設定

{\displaystyle a_{0}=1\qquad b_{0}={\frac {1}{\sqrt {2}}}\qquad t_{0}={\frac {1}{4}}\qquad p_{0}=1}

反復式
a, b が希望する精度（桁数）になるまで以下の計算を繰り返す。小数第n位まで求めるとき　log2 n 回程度の反復でよい。

{\displaystyle {\begin{aligned}a_{n+1}&={\frac {a_{n}+b_{n}}{2}}\\b_{n+1}&={\sqrt {a_{n}b_{n}}}\\t_{n+1}&=t_{n}-p_{n}(a_{n}-a_{n+1})^{2}\\p_{n+1}&=2p_{n}\end{aligned}}}

π の算出
円周率 π は、a, b, t を用いて以下のように近似される。

{\displaystyle \pi \approx {\frac {(a+b)^{2}}{4t}}}

上記を素直にコードにすると以下になって、実行してみると　pi = 3.140579250522169 になった

#include <math.h>
#include <iostream>
#include <iomanip>

void gauss_Legendre_algorithm()
{
   std::cout << "Start.." << std::endl;

   // 初期値
   double a0 = 1.0;
   double b0 = 1.0 / sqrt(2.0);
   double t0 = 1.0 / 4.0;
   double p0 = 1.0;

   // １回目
   double a1 = (a0 + b0) / 2.0;
   double b1 = sqrt(a0 * b0);
   double t1 = t0 - p0 * pow((a0 - a1), 2.0);
   double p1 = 2 * p0;

   // PI計算
   double pi = 0.0;
   pi = pow((a1 + b1), 2.0) / (4.0 * t1);
   std::cout << " pi = "
             << std::fixed
             << std::setprecision(15)
             << pi
             << std::endl;

   std::cout << "End.." << std::endl;
}

int main(int argc, char *argv[])
{
   gauss_Legendre_algorithm();
   return 0;
}

繰り返すことで精度があがるので６回くらい回せるコードに書き換えてみた。

#include <math.h>
#include <iostream>
#include <iomanip>

struct Parameter
{
   double a;
   double b;
   double t;
   double p;

   Parameter() : a(1.0), b(1.0 / sqrt(2.0)), t(1.0 / 4.0), p(1.0)
   {
   }
};

double gauss_Legendre_algorithm(Parameter &param)
{
   double a = (param.a + param.b) / 2.0;
   double b = sqrt(param.a * param.b);
   double t = param.t - param.p * pow((param.a - a), 2.0);
   double p = 2.0 * param.p;

   param.a = a;
   param.b = b;
   param.t = t;
   param.p = p;

   // pi計算
   return pow((a + b), 2.0) / (4.0 * t);
}

int main(int argc, char *argv[])
{
   std::cout << "Start.." << std::endl;
   Parameter param;
   double pi;

   for (int i = 1; i < 6; i++)
   {
      pi = gauss_Legendre_algorithm(param);
      std::cout << "pi = "
                << std::fixed
                << std::setprecision(20)
                << pi
                << std::endl;
   }

   std::cout << "End.." << std::endl;
   return 0;
}

結果は以下

pi = 3.14057925052216857509
pi = 3.14159264621354283875
pi = 3.14159265358979400418
pi = 3.14159265358979400418
pi = 3.14159265358979400418

どうも double では精度の限界があるようだ、それもそのはず　double (倍精度浮動小数点数) では少数は１５桁までしか表現できないらしい。

では、どうすれば１００万桁まで求められるか？
任意精度を保持できる「多倍長浮動小数点数」を使うと良さそうだ。
PythonやJavaなどは標準で実装されているようだが、 C++にはないので探してみた。

以下を使えば１億桁（それ以上も）表現できそうだ。

まずは１００万桁で実装してみた。

#include <math.h>
#include <iostream>
#include <iomanip>
#include <gmpxx.h> //brew install gmp
#include <string>

mpf_class power(mpf_class x, int n)
{
   mpf_class powval(1);
   for (int num = n; num != 0; --num)
   {
      powval *= x;
   }
   return powval;
}

struct Parameter
{
   mpf_class a;
   mpf_class b;
   mpf_class t;
   mpf_class p;

   Parameter() : a(1.0), b(1.0 / sqrt(mpf_class(2))), t(1.0 / 4.0), p(1.0)
   {
   }
};

mpf_class gauss_Legendre_algorithm(Parameter &param)
{
   mpf_class r;
   mpf_class a = (param.a + param.b) / 2;
   mpf_class b = sqrt(param.a * param.b);
   mpf_class t = param.t - param.p * power((param.a - a), 2);
   mpf_class p = 2 * param.p;

   r = power((a + b), 2) / (4 * t);

   param.a = a;
   param.b = b;
   param.t = t;
   param.p = p;

   return r;
}

void displayWrappedText(const std::string &text, int width)
{
   for (size_t i = 0; i < text.size(); i += width)
   {
      std::cout << text.substr(i, width) << std::endl;
   }
}

int main(int argc, char *argv[])
{
   std::cout << "Start.." << std::endl;

   // ★精度 ここを調整する
   int digits = 1000000;

   // 精度
   mpf_set_default_prec(digits * log2(10));
   int recursion_count = log2(digits);

   Parameter param;
   mpf_class pi;

   for (int i = 1; i <= recursion_count; i++)
   {
      pi = gauss_Legendre_algorithm(param);
   }

   std::cout << "有効桁 = " << digits << std::endl;
   std::cout << "計算回数 = " << recursion_count << std::endl;

   mp_exp_t exp;
   displayWrappedText(pi.get_str(exp), 100);

   std::cout << "End.." << std::endl;
   return 0;
};

結果　（長いので省略）

有効桁 = 1000000
計算回数 = 19
3141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067
9821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819
6442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127
3724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609
4330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491
2983367336244065664308602139494639522473719070217986094370277053921717629317675238467481846766940513
2000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279689258923
5420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185
.
.
.

円周率表を使って最初の１００桁まで調べてみたが、合っているので良しとする。
https://www.tstcl.jp/ja/randd/pi.php

さらに精度を上げていく

結構早かった、１００万桁まで（出力を除くと）１秒位だ、
調子にのって１０００万桁まで求めてみると計算時間はだいたい　２５秒
さらに調子にのって１億桁まで求めてみると計算時間はだいたい　６分

処理時間は指数関数的に伸びていくような感じではある。

M2 MacBookAir での結果です。マルチスレッド化やGPUを利用すれば更に高速化できるが、その場合はアルゴリズムを変更する必要があると思われる。Cudaを使って高速化してみたいところではある。