円周率とベルヌーイ数

はじめに

2022年 Simon Plouffe氏の論文「A formula for the n’th digit of 𝜋 and 𝜋^n」より

円周率の特定桁の数字がわかる、というものです。
ベルヌーイ数が登場する、円周率の公式です。

計算式　その１

下記のMartin Bauer氏の投稿によると、「^(-n)」ではなくて、「^(-2n)」じゃないかという指摘があります。
確かに、前者ではうまく計算できませんでした。後者ではうまくいく。

ベルヌーイ数

式中にBが登場し、コレはベルヌーイ数を表すそうです。
小数版と整数版のプログラム。
この数列の偶数番目を使う円周率の式、が冒頭の研究で出てくる式です。

void setup() {
  int N = 12;
  double[] B = calcBernoulli(N);

  for (int n = 0; n <= N; n++) {
    println("B_" + n + " = " + B[n]);
  }
  exit(); // 終了（ウィンドウ不要なら）
}

// ベルヌーイ数の計算（double近似）
double[] calcBernoulli(int nmax) {
  double[] B = new double[nmax + 1];
  B[0] = 1.0;
  for (int n = 1; n <= nmax; n++) {
    double sum = 0.0;
    for (int k = 0; k < n; k++) {
      sum += binomial(n + 1, k) * B[k];
    }
    B[n] = -sum / (n + 1);
  }
  return B;
}

// 二項係数 nCk
long binomial(int n, int k) {
  if (k < 0 || k > n) return 0;
  if (k == 0 || k == n) return 1;
  long res = 1;
  for (int i = 1; i <= k; i++) {
    res = res * (n - i + 1) / i;
  }
  return res;
}

B_0 = 1.0
B_1 = -0.5
B_2 = 0.16666666666666666
B_3 = -0.0
B_4 = -0.033333333333333305
B_5 = -7.401486830834377E-17
B_6 = 0.023809523809523885
B_7 = 6.938893903907228E-17
B_8 = -0.03333333333333364
B_9 = 1.7763568394002506E-16
B_10 = 0.07575757575757645
B_11 = -1.4802973661668755E-15
B_12 = -0.25311355311355305

void setup() {
  int N = 12;
  Fraction[] B = calcBernoulliFraction(N);

  for (int n = 0; n <= N; n++) {
    println("B_" + n + " = " + B[n]);
  }
  exit();
}

// 分数型クラス
class Fraction {
  long num, den;
  Fraction(long n, long d) {
    if (d == 0) throw new IllegalArgumentException("Zero denominator");
    // 約分と符号処理
    long g = gcd(Math.abs(n), Math.abs(d));
    n /= g;
    d /= g;
    if (d < 0) { n = -n; d = -d; }
    num = n;
    den = d;
  }
  Fraction add(Fraction f) {
    return new Fraction(num * f.den + f.num * den, den * f.den);
  }
  Fraction sub(Fraction f) {
    return new Fraction(num * f.den - f.num * den, den * f.den);
  }
  Fraction mul(Fraction f) {
    return new Fraction(num * f.num, den * f.den);
  }
  Fraction div(Fraction f) {
    return new Fraction(num * f.den, den * f.num);
  }
  public String toString() {
    if (num == 0) return "0";
    if (den == 1) return "" + num;
    return num + "/" + den;
  }
  long gcd(long a, long b) {
    return b == 0 ? a : gcd(b, a % b);
  }
}

// nまでのベルヌーイ数（分数型）を返す
Fraction[] calcBernoulliFraction(int nmax) {
  Fraction[] B = new Fraction[nmax + 1];
  B[0] = new Fraction(1, 1);
  for (int n = 1; n <= nmax; n++) {
    Fraction sum = new Fraction(0, 1);
    for (int k = 0; k < n; k++) {
      Fraction coef = new Fraction(binomial(n + 1, k), 1);
      sum = sum.add(coef.mul(B[k]));
    }
    B[n] = sum.div(new Fraction(-(n + 1), 1));
  }
  return B;
}

// 二項係数 nCk
long binomial(int n, int k) {
  if (k < 0 || k > n) return 0;
  if (k == 0 || k == n) return 1;
  long res = 1;
  for (int i = 1; i <= k; i++) {
    res = res * (n - i + 1) / i;
  }
  return res;
}

B_0 = 1
B_1 = -1/2
B_2 = 1/6
B_3 = 0
B_4 = -1/30
B_5 = 0
B_6 = 1/42
B_7 = 0
B_8 = -1/30
B_9 = 0
B_10 = 5/66
B_11 = 0
B_12 = -691/2730

計算式　その１　結果

特定の桁ではなくて、円周率そのものが計算結果に登場しました。
桁溢れが心配ですが、nが１増えると、２桁精度が高まるので、面白いと思いました。

計算式　その２

ベルヌーイ数の奇数が０で、偶数だけ使うことから、「その１」の式ができているようで、
「n=偶数」とするならば、以下のよりシンプルな表記ができるもよう。

式としては、等価ですが、こちらの式は、ベルヌーイ数がマイナスのときに、エクセル計算できないので、プラスにしました。（黄色セル）