1.4.1 円周率の計算

卒研で何人かの学生が円周率の数値計算を卒業研究レポートのテーマにした。

1976年に E. Salamin と R. P. Brent により発表された画期的な 円周率の計算アルゴリズム (算術幾何平均反復に基づく) は、 Gauss, Legendre によって知られていた楕円積分論の事実と、 1971年に発表された Schönhage–Strassen の高速乗算法を 組み合わせたものである。

これについては、 鎌田・吉本[4]や、 桂田[5]を見よ。

一松先生が数研通信に https://www.chart.co.jp/subject/sugaku/suken_tsushin/77/77-8.pdfというのを書いているが、あまり分かりやすくない。

2017年9月15日、「現象と数学」で少し話した。

「算術幾何平均で遊ぶ」
http://nalab.mind.meiji.ac.jp/~mk/20170915.pdf

/*
 *   g++ -I/opt/local/include test_pi_1mega.cpp -L/opt/local/lib -lmpfr
 *
 *   Mathematica があれば
 *     NumberForm[N[Pi, 1000000 + 1], DigitBlock -> 10]
 *   として、小数第100万位まで表示できるので、照らし合わせると良い。
 * 1000万桁100秒: 96.984u 5.236s 1:42.75 99.4%	0+0k 103+4io 103pf+0w
 * 1億桁        1634.801u 85.617s 29:25.04 97.4%	0+0k 207+11io 45pf+0w
 */

#include <iostream>
#include <string>
#include "mpreal.h"

using namespace std;

int main(void)
{
  int n, maxn = 50;
  using mpfr::mpreal;
  // const int digits = 1000 * 1000 * 180 + 50;
  const int digits = 1000 * 1000 + 50;
  mpreal::set_default_prec(mpfr::digits2bits(digits));
  mpreal two, a, b, c, ap1, c2, mypi, p, s, old, error, d50;
  string str;
  int i;
  a = 1; two = 2; b = 1 / sqrt(two); c = b;
  s = 0;
  p = 1;
  for (n = 0; n <= maxn; n++) {
    old = mypi;
    c2 = c * c;
    s += p * c2;
    p *= 2;
    ap1 = (a + b) / 2;
    b = sqrt(a * b);
    c = c2 / (4 * ap1);
    a = ap1;
    mypi = 2 * a * a / (1 - s);
    error = abs(mypi - old);
    cout.precision(10);
    cout << "n = " << n << ", 誤差=" << error << endl;
    if (error == 0)
      break;
  }
  cout.precision(digits);
  str = mypi.toString();
  cout << str.substr(0, 2) << endl;
  for (n = 0; n < digits; n += 50) {
    cout << "  " << str.substr(n + 2, 50) << endl;
  }
  return 0;
}

n = 0, 誤差=2.914213562
n = 1, 誤差=0.2263656881
n = 2, 誤差=0.001013395691
n = 3, 誤差=7.376250956e-09
n = 4, 誤差=1.831306085e-19
n = 5, 誤差=5.472109146e-41
n = 6, 誤差=2.40612213e-84
n = 7, 誤差=2.308580715e-171
n = 8, 誤差=1.058626539e-345
n = 9, 誤差=1.110954934e-694
n = 10, 誤差=6.111788737e-1393
n = 11, 誤差=9.244416653e-2790
n = 12, 誤差=1.057232713e-5583
n = 13, 誤差=6.913088685e-11172
n = 14, 誤差=1.477814384e-22348
n = 15, 誤差=3.376546688e-44702
n = 16, 誤差=8.813373486e-89410
n = 17, 誤差=3.002256862e-178825
n = 18, 誤差=1.741917617e-357656
n = 19, 誤差=2.931948617e-715319
n = 20, 誤差=4.153205794e-1430645
n = 21, 誤差=4.166846002e-2861297
n = 22, 誤差=2.097130002e-5722601
n = 23, 誤差=2.656018619e-11445210
n = 24, 誤差=2.130161926e-22890428
n = 25, 誤差=6.85086825e-45780865
n = 26, 誤差=3.543085308e-91561738
n = 27, 誤差=0
3.
  14159265358979323846264338327950288419716939937510
  58209749445923078164062862089986280348253421170679
  82148086513282306647093844609550582231725359408128
(以下略)
3113.464u 197.230s 56:28.35 97.7%0+0k 74+35io 48pf+0w