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;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;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).
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.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274274663919320030599218174135966290435729003342952605956307381323286279434907632338298807531952510190115738341879307021540891499348841675092447614606680822648001684774118537423454424371075390777449920695517027618386062613313845830007520449338265602976067371132007093287091274437470472306969772093101416928368190255151086574637721112523897844250569536967707854499699679468644549059879316368892300987931After 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;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.