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.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 |
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.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 |
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.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 |
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

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 |
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.

Design. Code. Compile. Deploy.

Start Free Trial Upgrade Today

Free Delphi Community Edition Free C++Builder Community Edition

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

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.

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!

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

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