Getting Big with Pi in Delphi

When I was very young and first learned about Pi they told me 3.14 was a good approximation, but it was an irrational number that went on forever, and people with computers were able to calculate many digits. I had a computer (Commodore Vic 20) and wanted to see how many digits of Pi I could calculate. This is when I discovered floating point number precision limitations, which foiled my first attempt. The idea of calculating more digits of Pi has always itched in the back of my mind, but CPU based floating point precision was never enough.

In my 2023 Pi day post I experimented with Rudy Velthuis’ arbitrary precision big numbers library and still only got 20 digits, but I wanted more. As you may know, Rudy Velthuis passed away a few years ago, so his library has remained unmaintained for four years. There were a few forks, and a pull request, but no new canonical home of his library. First thing I did was merge all the pull requests and forks into a new fork in the Turbo Pack organization on GitHub, then I collected up Rudy’s articles on his libraries, and added them to the wiki, and did some additional general clean up and updating. I worked with Roman Kassebaum, the maintainer of Turbo Pack, to make sure it built and installed on more than just my machine.

[ The new canonical home of
Rudy’s Big Numbers library
github.com/TurboPack/RudysBigNumbers ]

Once I was sure Rudy’s Big Numbers libraries were working correctly, I went about rewriting my Pi calculation implementations. I used ChatGPT and Bing at first, but they never got as much accuracy as I wanted. Instead I ended up searching for more details and examples in other languages, and attempted to translate them to Delphi. I ended up with two results, both calculate a verified 100k digits of Pi, and with a respectable speed.

Table of Contents

Chudnovsky and BigDecimals

The Chudnovsky algorithm was published by the Chudnovsky brothers in 1988. It was used in the world record calculation of 100 trillion digits of Pi on March 21, 2022.

You specify the number of places you want it to calculate too, and then it starts with an estimate, and continues improving it until there are no changes out to that many decimal places. Internally is uses an additional six digits of precision, but then rounds it down to the specified precision once it settles on a number. The 32-bit version ran out of memory, but with 64-bit Windows version I reached 100K digits.

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’s implementation of BigDecimals points out it isn’t very efficient, so I was curious if I could calculate faster by switching to BigIntegers.

Bailey-Borwein-Plouffe and BigIntegers

Discovered in 1995, the BBP algorithm works with integers to generate successive digits of pi. So it doesn’t calculate Pi as a whole, but just each individual digit. My implementation adds those digits to an array.

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 < n*t then
    begin
      result[idx] := n.AsInt64; // It is just a byte
      inc(idx);
      //write(n.ToString[1]);
      if idx >= 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 < n*t then

begin

result[idx] := n.AsInt64; // It is just a byte

inc(idx);

//write(n.ToString[1]);

if idx >= 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;

You can easily modify it to just output each digit as it goes, instead of storing it in the array. If you do then you will see the rate of generation slows down the deeper you go. I thought about adding a call back for each digit, but more likely I will modify it to generate digits in a range, so you don’t need to start at the top each time. But I’ll leave that for another day.

Calculation Time

Performance will vary based on your computer hardware, but I wanted to see how it scaled and which algorithm was faster. I ran each on Win32, Win64, and Linux64, using the same hardware (Linux was under Ubuntu with WSL2).

delphi big number pi computation time 9809679

Win64 was the fastest, with BBP implementation using BigIntegers being the leader.

Mobile Pi Calculation

Sometimes you need to calculate a few thousand digits of Pi while you are away from your computer. So to make sure this all works under Android, I created a FMX version. It was actually faster on Android than I expected, but I didn’t bother benchmarking it or seeing how many digits it could generate. In theory it is only limited by your available memory.

Getting the Code

I added all of my code to the samples folder on GitHub. There are also unit tests that verify both algorithms work as expected. It uses a series of hashes to test it to various numbers of digits, as well as a a saved 100K digits of Pi. If you don’t want to install from GitHub then stay tuned to GetIt where we will have an installer soon.

More than Pi

Once I figured out Rudy’s libraries I wanted to try other irrational numbers besides Pi. They were all pretty easy to generate, but Euler’s number required a little work.

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.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274274663919320030599218174135966290435729003342952605956307381323286279434907632338298807531952510190115738341879307021540891499348841675092447614606680822648001684774118537423454424371075390777449920695517027618386062613313845830007520449338265602976067371132007093287091274437470472306969772093101416928368190255151086574637721112523897844250569536967707854499699679468644549059879316368892300987931

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.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274274663919320030599218174135966290435729003342952605956307381323286279434907632338298807531952510190115738341879307021540891499348841675092447614606680822648001684774118537423454424371075390777449920695517027618386062613313845830007520449338265602976067371132007093287091274437470472306969772093101416928368190255151086574637721112523897844250569536967707854499699679468644549059879316368892300987931

After a few different algorithms and multiple iterations I ended up with the following for Euler’s number

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;

The process is similar to my Chudnovsky implementation in that it continues to improve it’s estimate until there are no changes within the specified number of digits. The formula does use factorials, but does them in sequence, which is much faster. I verified it to 10K digits as well.

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

5 Comments


David Green

April 12, 2023 at 8:36 pm

Hi Jim
I enjoyed your blog about pi. One of the first, if not the first, efforts at multiple length calculation was in 1951, on the Pilot ACE in England. They calculated e to 306 decimal figures, not a bad feat given the computer then had less than 300 words of storage for code and data. If you are interested, I have resurrected the lost code, on http://uraone.com/e306/e306.html
Regards
David Green
- Reply
  
  Jim McKeeth
  
  April 13, 2023 at 1:54 pm
  
  Thanks, glad there is someone else who enjoyed it. This is something I’ve been interested in for a long time, and all of my friends and family had to put up with me going on and on about it.

Charles Stack

April 13, 2023 at 5:32 pm

I’ve been using Rudy’s BigInteger code for some cryptographic and combinatoric code I’m working on.

I stayed away from BigDecimals because of the performance and, well, I didn’t really need it for what I was doing.

I’ll have to take a look at the fork you have compiled. Perhaps, in time, I’ll check in some of the code I’ve been working on.

Nice job, Jim!

breitsch breitsch

April 14, 2023 at 3:05 pm

Fascinating is the combination of Bigint with Bigfloat in rounding. Did also code and test a script version in maXbox as a combination of bigint and bigfloat with another precompiled runtime BigX library from Delphi for Science and Delphi for Fun:
https://github.com/breitsch2/RudysBigNumbers/blob/main/Samples/Pi/1207_cryptosystem2test_bigpi1.txt
The hash of Pi1000 is the same so it works, speed is another story as the script run as interpreter.
By the way Pi is an irrational number, and its sequence is as we know infinite and non-repeating.
Feynman thought that if he could memorize Pi up to the 762nd decimal place, he could trick people into thinking it was rational and say “999999…
Max

breitsch breitsch

April 17, 2023 at 4:37 am

Hi Jim
Just a note concerning the file: MathConstants/Euler10k.txt
its more than 10k , its 11141
abigpi2:= BIGE2;
SetLength(abigpi2, Length(abigpi2) – 1140);
Greetings, Max

Getting Big with Pi in Delphi

Chudnovsky and BigDecimals

Bailey-Borwein-Plouffe and BigIntegers

Calculation Time

Mobile Pi Calculation

Getting the Code

More than Pi

About author

Jim McKeeth

5 Comments

Leave a ReplyCancel reply

Search

Something Fresh

How To Live Track The ISS International Space Station

Read The Webinar: Why You Should Leave Delphi 7 Once And For All

Upgrading C++Builder 12.2: Tip #7, bcc64x tips!

Popular Posts

RAD Studio 12.2 Athens Inline Patch 1 Available

Announcing the Availability of RAD Studio 12.2 Athens

Delphi 12 And C++Builder 12 Community Editions Released!

InterBase ODBC Driver on GitHub

Embarcadero Partners with Raize Software for KSVC Maintenance

Categories

Unknown Feed

Unknown Feed

Categories

Useful Links

Follow us

Getting Big with Pi in Delphi

Chudnovsky and BigDecimals

Bailey-Borwein-Plouffe and BigIntegers

Calculation Time

Mobile Pi Calculation

Getting the Code

More than Pi

About author

5 Comments

Leave a ReplyCancel reply

Join Our Global Developer Community

Search

Something Fresh

Popular Posts

Categories

Categories

Useful Links

Follow us