Calculating pi with C#
published: Sat, 15-Mar-2008 | updated: Wed, 2-Jan-2019
A couple of evenings ago, my wife was reading an article about a federal prosecutor, when she suddenly exclaimed "How many digits of π do you know? Because this lawyer knows 500 digits." I had to admit that I'd only memorized π to eight decimal places, 3.14159265, but I immediately jumped in with "I know how to write a program that can print out thousands of digits of π." So she threw down the gauntlet.
OK, so my boast was a bit of a stretch: although I'd never written such a program before, I did remember the basics, something about Taylor series.
Back in 1706, John Machin, then a professor of astronomy at Gresham College in London, devised a formula for π that would be used for hundreds of years afterwards to calculate ever more digits of π (he himself used it to calculate 100 decimals of π):
π/4 = 4·arctan(1⁄5) - arctan(1⁄239)
Doesn't look like much improvement — you now have the problem of calculating two arctan values instead — but there's a celebrated rapidly converging infinite Taylor series for doing just that:
arctan(x) = x - x3⁄3 + x5⁄5 - x7⁄7 + x9⁄9 - ... for |x| ≤ 1
It converges faster the closer x is to 0, so in Machin's formula for π, the arctan(1/239) term converges much faster than the arctan(1/5) term.
(The Taylor series were invented by Brook Taylor, who was a contemporary of Machin's and a mathematician at Cambridge. The Taylor series were, at that time, a brand new analytic tool.)
OK, now that we have enough mathematics, we have to think about devising a program to calculate π. The first thing to note is that we can't use any of the primitive types in C#: they just don't have enough precision (the double type, for example, only has 15 significant digits, and we want to calculate, say, 10000). We need a large precision number library, but there's nothing like that in the .NET Framework, and so the first thing we'll have to do is write a BigNumber class.
If you look at the series for arctan, you'll see that you can use a shortcut. Each term in the series consists of a power-of-x, divided by some constant. The powers-of-x increase by taking the previous power-of-x and multiplying it by x2. Since we shall be using integer reciprocals, we can see that to get the next power of x, we merely divide by the square of that integer (in our case, 5 or 239). In other words, our BigNumber class doesn't have to have a "divide this BigNumber by another BigNumber" method — at which point I would have discretely dropped the entire project, thank you very much — but can make do with a "divide by an integer" method instead. This, to be utterly frank, is a much easier proposition. And we can use it for the divide-by-a-constant parts of each term as well.
The BigNumber class will need addition and subtraction too, and also a multiply by an integer operation.
What we'll do is make our BigNumber class use a fixed point representation. The "digits" in our representation are full unsigned 32-bit integers (we'll be doing calculations base 4,294,967,296, if you like). There will be one "digit" before the radix point, and some large number of digits after the radix point. We shall in due course have to write a conversion routine to get our answer base 4,294,967,296 into our usual base 10.
There is one more issue we have to solve first. We'd like to specify the maximum number of decimal digits to our BigNumber class in order to size it, rather than the number of digits base 4,294,967,296. Since that's 232, which is just a bit more than 109.6 we can divide the maximum number of decimal digits by 9.6 (or multiply by .104) to get the number of digits base 232. Add one for the digits before the decimal point, plus one for rounding errors in the conversion to base 10, and we should be good to go. So, for example, to get 10000 digits of π after the decimal point, we would need to allocate an array of 1042 UInt32 integers.
First up then is the constructor.
public class BigNumber { private UInt32[] number; private int size; private int maxDigits; public BigNumber(int maxDigits) { this.maxDigits = maxDigits; this.size = (int) Math.Ceiling((float)maxDigits * 0.104) + 2; number = new UInt32[size]; } public BigNumber(int maxDigits, UInt32 intPart) : this(maxDigits) { number[0] = intPart; for (int i = 1; i < size; i++) { number[i] = 0; } }
As you can see there are two: the first allocates the requires UInt32 array and just initializes the number to zero (by default), and the second, initializes the number to an integer.
Next we'll write the Add() method for adding two BigNumbers. Just in case, we'll make sure that the two numbers have the same size, and they're not the same object (we'll follow this check with other methods).
private void VerifySameSize(BigNumber value) {
if (Object.ReferenceEquals(this, value))
throw new Exception("BigNumbers cannot operate on themselves");
if (value.size != this.size)
throw new Exception("BigNumbers must have the same size");
}
public void Add(BigNumber value) {
VerifySameSize(value);
int index = size - 1;
while (index >= 0 && value.number[index] == 0)
index--;
UInt32 carry = 0;
while (index >= 0) {
UInt64 result = (UInt64)number[index] +
value.number[index] + carry;
number[index] = (UInt32)result;
if (result >= 0x100000000U)
carry = 1;
else
carry = 0;
index--;
}
}
In essence the code merely does the same kind of addition we learned at school: add up the digits from right to left, and remember to take care of the carry as we do so. It's complicated a little bit in that the addition of the digits has to be done in 64-bit arithmetic, but other than that, nothing too spectacular.
Next is subtraction. Here we assume that the value being subtracted is less than the one being subtracted from. In other words, I'm assuming that there are no negative numbers; yet again, a simplifying tactic since this is not going to be a universal BigNumber class useful for a whole bunch of calculations, just one that's good for calculating π (and maybe other simple Taylor series, like calculating e).
public void Subtract(BigNumber value) {
VerifySameSize(value);
int index = size - 1;
while (index >= 0 && value.number[index] == 0)
index--;
UInt32 borrow = 0;
while (index >= 0) {
UInt64 result = 0x100000000U + (UInt64)number[index] -
value.number[index] - borrow;
number[index] = (UInt32)result;
if (result >= 0x100000000U)
borrow = 0;
else
borrow = 1;
index--;
}
}
Again, this is standard school-level subtraction. To ensure no negative numbers are calculated as I go through, I assume that a borrow is going to be needed with each digit position. I do the subtraction of each digit, and then check to see if the borrow actually was used or not.
We now visit multiplication by an integer.
public void Multiply(UInt32 value) {
int index = size - 1;
while (index >= 0 && number[index] == 0)
index--;
UInt32 carry = 0;
while (index >= 0) {
UInt64 result = (UInt64)number[index] * value + carry;
number[index] = (UInt32)result;
carry = (UInt32)(result >> 32);
index--;
}
}
Like before, nothing too new here: it's just like multiplying an ordinary decimal number by 7, say. We work from the right to the left using 64-bit arithmetic, making sure we calculate the carry with every digit multiplication. We don't worry about overflow, although it's simple to spot, since, again, we're limiting this to our very narrow scenario.
Division by an integer, next.
public void Divide(UInt32 value) {
int index = 0;
while (index < size && number[index] == 0)
index++;
UInt32 carry = 0;
while (index < size) {
UInt64 result = number[index] + ((UInt64)carry << 32);
number[index] = (UInt32) (result / (UInt64)value);
carry = (UInt32)(result % (UInt64)value);
index++;
}
}
Here we work from the left to the right, dividing each 32-bit digit by the divisor, and working out the carry using the mod operator. In essence, and I repeat myself, much like we learned long division at school.
Now we have the basic operations written on a BigNumber, we can now write a routine to print such a number in decimal.
public void Assign(BigNumber value) {
VerifySameSize(value);
for (int i = 0; i < size; i++) {
number[i] = value.number[i];
}
}
public string Print() {
BigNumber temp = new BigNumber(maxDigits);
temp.Assign(this);
StringBuilder sb = new StringBuilder();
sb.Append(temp.number[0]);
sb.Append(System.Globalization.CultureInfo.CurrentCulture.NumberFormat.CurrencyDecimalSeparator);
int digitCount = 0;
while (digitCount < maxDigits) {
temp.number[0] = 0;
temp.Multiply(100000);
sb.AppendFormat("{0:D5}", temp.number[0]);
digitCount += 5;
}
return sb.ToString();
}
This routine works as follows. Convert the integer part to decimal (using the standard ToString() method) and add the decimal point. Now go into a loop until we've extracted the maximum number of digits the BigNumber holds. Set the integer part to zero, multiply the number that's left by 100000 to get the next five decimal digits, and convert the integer part to those five decimal digits (we must make sure we keep the leading zeros!). Whip round the loop again and again until we have calculated the required number of decimal digits.
We're all set now for the pièce de résistance: calculating the arctan of a number, where that number is assumed to be the reciprocal of an integer.
public bool IsZero() {
foreach (UInt32 item in number) {
if (item != 0)
return false;
}
return true;
}
public void ArcTan(UInt32 multiplicand, UInt32 reciprocal) {
BigNumber X = new BigNumber(maxDigits, multiplicand);
X.Divide(reciprocal);
reciprocal *= reciprocal;
this.Assign(X);
BigNumber term = new BigNumber(maxDigits);
UInt32 divisor = 1;
bool subtractTerm = true;
while (true) {
X.Divide(reciprocal);
term.Assign(X);
divisor += 2;
term.Divide(divisor);
if (term.IsZero())
break;
if (subtractTerm)
this.Subtract(term);
else
this.Add(term);
subtractTerm = !subtractTerm;
}
}
We need two temporary BigNumbers here, one to keep track of the
powers-of-x parts, X, and the other to calculate each
term of the series, term. We start off by calculating the
first term of the series and putting it into X. We can
then square the integer reciprocal. And then we're off into the loop.
We divide the powers-of-x part by the square of the reciprocal, and
assign it to the term temporary. We then divide
term by the new divisor (this increments by 2 each time
through the loop to count off the odd numbers) and then we either add
or subtract this new term of the series into our running total,
alternating each time through the loop. The loop ends once the new
value for term has underflowed and is zero.
Finally, fanfare, we can write the actual calculation of π.
class Program { static void Main(string[] args) { BigNumber x = new BigNumber(1000); BigNumber y = new BigNumber(1000); x.ArcTan(16, 5); y.ArcTan(4, 239); x.Subtract(y); Console.WriteLine(x.PrintAsTable()); Console.ReadLine(); } }
(Here PrintAsTable() produces a very stylized output of the value as a grid of 10 digit blocks.) The output for the first 1000 decimal digits of π is:
3. 1415926535 8979323846 2643383279 5028841971 6939937510 5820974944 5923078164 0628620899 8628034825 3421170679 8214808651 3282306647 0938446095 5058223172 5359408128 4811174502 8410270193 8521105559 6446229489 5493038196 4428810975 6659334461 2847564823 3786783165 2712019091 4564856692 3460348610 4543266482 1339360726 0249141273 7245870066 0631558817 4881520920 9628292540 9171536436 7892590360 0113305305 4882046652 1384146951 9415116094 3305727036 5759591953 0921861173 8193261179 3105118548 0744623799 6274956735 1885752724 8912279381 8301194912 9833673362 4406566430 8602139494 6395224737 1907021798 6094370277 0539217176 2931767523 8467481846 7669405132 0005681271 4526356082 7785771342 7577896091 7363717872 1468440901 2249534301 4654958537 1050792279 6892589235 4201995611 2129021960 8640344181 5981362977 4771309960 5187072113 4999999837 2978049951 0597317328 1609631859 5024459455 3469083026 4252230825 3344685035 2619311881 7101000313 7838752886 5875332083 8142061717 7669147303 5982534904 2875546873 1159562863 8823537875 9375195778 1857780532 1712268066 1300192787 6611195909 2164201989
(I am indebted to J.W. Stumpel and his page on calculating π, especially his simple C program, for this article.)
