本ブログは、こちらの続きで前回 Rudy VelthuisのBigNumbersライブラリを使用した円周率の計算を行いましたが、残念ながら20桁程度しか正確に取得できませんでした。今回はもっと大きな桁の円周率の計算が行えるように、まずBigNumbersライブラリのアップデートを行うことにしました。前回でも少し触れましたがBigNumbersライブラリの作者であるRudy Velthuisは数年前に亡くなっており、Rudyのライブラリは4年間メンテナンスされていませんでした。その間、いくつかのフォーク(派生)やプルリクエストはありましたが、Rudyのライブラリの新しい正規なホームではありません。そこで筆者は、すべてのプルリクエストとフォークをGitHubのTurbo Pack Organizationの新しいフォークにマージしました。そしてRudyのライブラリに関する記事を集めてWikiに追加し、さらに一般的なクリーンアップとアップデートを行いました。

その苦労の甲斐があり、Turbo PackのメンテナであるRoman Kassebaumと協力して、筆者のマシン以外にもビルドしてインストールできるようにしました。

RudyのBig Numbersライブラリの新しい正規のホームは、こちらです。

RudyのBig Numbersライブラリが正しく動作していることを確認した後、円周率計算の実装を書き直しました。最初はChatGPTとBingに頼ってヒントやアルゴリズムを模索しましたが、期待した精度の結果が得られなかったため、他のコンピュータ言語で紹介されている例題を見つけて、それをDelphiコードへコンバートすることを試みました。その結果、10万桁の円周率を求めるコードで、かなりの速度で計算できるようになりました。以下でもっと詳しく紹介いたします。

チュドノフスキー公式 × BigDecimals

チュドノフスキーの公式は、1988年にチュドノフスキー兄弟によって発表され、2022年3月21日に世界記録となった100兆桁の円周率の計算に使用されました。

計算する桁数を指定すると、推定値から始まり、小数点以下の桁数が変わらなくなる近似値まで求め続けます。内部的にはさらに6桁の精度が使用されますが、数値が決まると指定された精度に丸められます。DelphiでビルドしたWindowsアプリでは、32ビット版では途中でメモリ不足になりましたが、64ビット版では10万桁まで達することができました。

  
function Chudnovsky(Places: Integer):BigDecimal;
begin
// Use +6 internally for calculations
  var internalPrecision := MaxInt;
  if Places <= MaxInt - 6 then
    internalPrecision := Places + 6;
 
  var lastSum: BigDecimal;
  var t := BigDecimal.Create(3);
  var sum := BigDecimal.Create(3);
  sum.DefaultPrecision := internalPrecision;
  sum.DefaultRoundingMode := rmFloor;
  var n := BigInteger.One;
  var d := BigInteger.Zero;
  var na: UInt64 := 0;
  var da: UInt64 := 24;
  while sum  lastSum do
  begin
    lastSum := sum;
    n := n + na;
    na := na + 8;
    d := d + da;
    da := da + 32;
    t := ((t * n)/d);
    sum := (sum + t).RoundToPrecision(internalPrecision);
  end;
  Result := sum.RoundToPrecision(Places);
end;

RudyのBigDecimalsの実装では、あまり効率的でないことが指摘されているので、BigIntegersに切り替えることで、もっと高速に計算できるかどうか興味がありましたので、実際に試してみました。以下で紹介いたします。

ベイリー＝ボールウェイン＝プラウフ公式 × BigIntegers

1995年にサイモン・プラウフによって発見されたBBP公式は、整数を用いて円周率の連続した桁を生成します。そのため円周率全体を計算するのではなく、それぞれの桁を計算します。以下のコード実装では、これらの桁を配列に追加しています。

function BBPpi(Places: UInt64): TDigits;
// Bailey-Borwein-Plouffe
begin
  SetLength(Result, Places);
 
  var idx: Uint64 := 0;
  var q := BigInteger.One;
  var r := BigInteger.Zero;
  var t := BigInteger.One;
  var k := BigInteger.One;
  var n := BigInteger.Create(3);
  var l := BigInteger.Create(3);
 
  while true do
  begin
    if 4*q+r-t = places then break;
      var newR := 10 * (r - n * t);
      n := ((10 * (3 * q + r)) div t) - 10 * n;
      q := q * 10;
      r := newR;
    end
    else
    begin
      var newR := (2 * q + r) * l;
      var newN := (q * (7 * k)+2+(r * l)) div (t * l);
      q := q * k;
      t := t * l;
      l := l + 2;
      k := k + 1;
      n := newN;
      r := newR;
    end;
  end;
end;

配列に格納するのではなく、各桁をそのまま出力するように修正することは簡単にできますが、その場合、深い桁に行くほど生成速度が遅くなります。各桁ごとにコールバックを追加することも考えましたが、それよりも毎回先頭から始める必要がないように、範囲内の桁を生成するように変更するほうが良いと考えました。しかしそれを試すのは、また別の機会にしたいと思います。

各プラットホームの演算時間は?

パフォーマンスはコンピュータのハードウェアによって異なりますが、どのようにスケーリングされ、どのアルゴリズムが最速なのか確認したかったので、同じハードウェア構成を使用し、Win32、Win64、および Linux64 の各プラットホームでアプリケーションを実行し、演算時間を計測しました（Linuxは、UbuntuでWSL2を使用）。

上図のグラフのようにWin64が最速で、BigIntegerを使用したBBPの実装が一番でした。

モバイルデバイスでの円周率計算

デスクトップPCではなくモバイルデバイスで、円周率の数千桁を計算する需要もあるかもしれません。

そこでAndroid上で動作させるためのFMX版のアプリを作成しました。実際にAndroidでは予想以上に高速でした。わざわざベンチマークを取ったり、何桁の計算まで可能か確認していませんが、理論的には、おそらく利用可能なメモリ量に依存します。

サンプルコードの入手

筆者はGitHubのSamplesフォルダに今回紹介したすべてのコードをアップロードしています。また、今回紹介した2つのアルゴリズムが期待通りに動作することを検証するためのユニットテストも用意しています。様々な桁数のハッシュを使用してテストし、100K桁の円周率も保存しています。もしGitHubからインストールしたくない場合は、近日中にGetItから入手できるようにインストーラを用意しますので、どうかお待ちください。

オイラー数

Rudyのライブラリを理解したら、円周率以外の別の無理数も試してみたくなりました。円周率は比較的簡単に生成できましたが、オイラー数は少し作業が必要でした。

Precision: 200
√2, Square root of 2, Pythagoras constant:
1.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727
√3, Square root of 3, Theodorus' constant:
1.7320508075688772935274463415058723669428052538103806280558069794519330169088000370811461867572485756
√5, Square root of 5:
2.2360679774997896964091736687312762354406183596115257242708972454105209256378048994144144083787822749
φ, Phi, Golden ratio, (1 + √5)/2:
1.61803398874989484820458683436563811772030917980576286213544862270526046281890244970720720418939113745
Common logarithm of 2:
0.301029995663981
Natural logarithm of 2:
0.693147180559945
------
Euler's Number to 500 digits
2.71828182845904523536028747135266249775724709369995957496696762772407663035354759457138217852516642

いくつかの異なるアルゴリズムで何度も繰り返した結果、オイラー数の実装コードは、以下の通りです。

function Euler(const Digits: Integer): BigDecimal;
begin
  Result := 1;
  var InternalPrecision := MaxInt;
  if Digits <= MaxInt - 2 then
    InternalPrecision := Digits + 2;
 
  Result.DefaultPrecision :=InternalPrecision;
  var lastFactorial := BigInteger.One;
  var lastResult := BigDecimal.Zero;
  var iteration: UInt64 := 1;
  while true do
  begin
    lastResult := Result;
    lastFactorial := lastFactorial * iteration;
    Result := (Result + BigDecimal.One / lastFactorial).RoundToPrecision(InternalPrecision);
    if lastResult = Result then break;
    inc(Iteration);
  end;
  //Writeln(Format('%d digits took %d iterations',[Digits, iteration]));
  Result.DefaultRoundingMode := rmFloor;
  Result := Result.RoundToPrecision(Digits);
end;

今回試した手順では、アルゴリズム的には上記で紹介したチュドノフスキーの実装に似ています。 この式は階乗を使いますが、順番に行うので、さらに高速で、10K桁まで検証しています。