Delphiで、もっと大きな桁の円周率を計算する

本ブログは、こちらの続きで前回 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;

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;

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;

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

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

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

delphi big number pi computation time 9809679 2

上図のグラフのように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

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;

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桁まで検証しています。

Reduce development time and get to market faster with RAD Studio, Delphi, or C++Builder.
Design. Code. Compile. Deploy.
Start Free Trial Upgrade Today

Free Delphi Community Edition Free C++Builder Community Edition

Delphiで、もっと大きな桁の円周率を計算する

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

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

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

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

サンプルコードの入手

オイラー数

コメントコメントをキャンセル

検索

新着記事

RAD Studio 12.1 Athens Patch 1 のリリース

新しいEmbarcadero Quality Portal 運用開始のお知らせ

RAD Studio 12.1 Athensの提供開始

Popular Posts

DelphiおよびC++Builder 11.3 Community Editionがリリースされました!

RAD Studio 11.2 Alexandria リリースのお知らせ

RAD Studio 11.1 Alexandriaリリースのお知らせ

Windows11やApple M1に対応した新バージョン『RAD Studio 11 Alexandria』提供開始

RAD Studio 11新機能ウェビナーを9月14日に開催

カテゴリー

人気の記事

Google Book APIを使ってISBNコードから書籍情報を取得する

基礎から学べるFireDACデータアクセス再入門（第３回）

最新のC++Buidlerで旧バージョンのプロジェクトを利用する方法

RAD Studio 10.4.2で提供される2つの新しいVCLコントロール

RAD Studio 10.4の新機能: 新しいVCL TEdgeBrowserコンポーネント

不明なフィード

不明なフィード

不明なフィード

Categories

エンバカデロをフォロー

jp custom

Delphiで、もっと大きな桁の円周率を計算する

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

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

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

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

サンプルコードの入手

オイラー数

コメントコメントをキャンセル

Join Our Global Developer Community

検索

新着記事

Popular Posts

カテゴリー

人気の記事

Categories

エンバカデロをフォロー

jp custom