It’s been a few years since I’ve had π day here on the C blog. Because I schedule posts for Saturdays, only thrice has a post fallen on π day: In 2015, 2020, and now in 2026.
Today is known as π day because it’s 3/14, which are the first three digits of π: 3.14159… So why not code yet another program that uses some clever mathematical trick to generate the value of π?
The secret, of course, is to locate some famous mathematician’s efforts to use a strange calculation that generates the value of π. In 2020, I used the Wallis Product. This time, I’m thieving the Leibniz Formula. Here is the mumbo jumbo artfully depicted in Figure 1:
Figure 1. The Leibniz formula for π.
Do you see a pattern? The denominators in the fractions show odd numbers: 3, 5, 7, 9, and up. The signs alternates between plus and minus. The trailing three dots denotes an infinite series.
My feeble math brain picked up on a few things right away.
First, the value 1 can be re-written as 1/1, shown in Figure 2.
Figure 1. Re-writing the value 1 as 1/1 in the Leibniz formula for π.
Second, because the result is written as π/4, it’s possible to multiple both sides of the equation by four. The result is shown in Figure 3:
Figure 3. After multiplying both sides by 4, the equation calculates the full value of π and is easier to code in C.
To work through this series you must set a denominator value starting with one and then increasing by two each operation, resulting in the series of odd numbers at the bottom of the fraction. The top value (the numerator) is always going to be four.
Further, the sign must be flipped between each operation, first subtracting the next value and then adding it. I covered this operation in an Exercise blog post from 2024.
My attempt to use the Leibniz formula to calculate the value of π, is shown in this code:
2026_03_14-Lesson.c
#include <locale.h>
#include <wchar.h>
#define PI_CHAR 0x03c0
#define ACCURACY 500000
int main()
{
float pi = 0.0;
int loop,denominator,sign;
/* use Leibniz formula to calculate pi */
denominator = 1;
sign = 0;
for( loop=0; loop<ACCURACY; loop++ )
{
pi += (sign%2?-1:1) * 4.0/denominator;
denominator+=2;
sign++;
}
/* output wide character 'pi' */
setlocale(LC_ALL,"");
wprintf(L"%lc = %.5f\n",PI_CHAR,pi);
/* clean-up */
return 0;
}
This program outputs wide characters (Unicode), so the stdio.h header file isn’t needed. Instead, the headers required for wide characters are used.
The two defined constants are the π Unicode character, PI_CHAR, and ACCURACY, which sets the number of times the main loop repeats. The greater the value of ACCURACY, the more accurate the result, though 500,000 is good enough for this program. (And it’s good enough for plotting the course of Jupiter’s orbit to within a few meters.)
The code uses a for loop to repeat ACCURACY times. Each time, the value of variable pi is increased by the value of sign (which determines whether to add or subject the next value) multiplied by 4.0 divided by the value of denomniator:
pi += (sign%2?-1:1) * 4.0/denominator;
The expression 4.0/denomniator need not be typecast to a float, at least not when using the clang compiler.
The value of denominator is increased by two: denominator+=2;, and sign is incremented: sign++;
The loop repeats, honing the value of pi with each iteration. When it’s done, the value of variable pi is output:
π = 3.14159
Though I may never understand the math, I do enjoy coding these π exercises — so much so, that I do another one next week, even though it’s not π day.
I implemented this a while ago and with 400000 iterations I got 3.14159(01)535897932(59).
The digits between () are wrong.
I handled the alternating +- a bit differently. It has to start positive as it’s flipped before each use.
long double multiplier = 1.0L;
long double denominator = 1.0L;
.
.
.
pi += ( (4.0L / (denominator += 2.0L)) * (multiplier *= -1.0L) );
When is the next Saturday?
Well, it’s Saturday again next week. But the next Saturday March 14 happens in 2037. I don’t know if I’ll still be doing this then….
I’ve just checked my diary and you’re right, there’s another Saturday this time next week. What an amazing coincidence!
Maybe you could write a program to find every pi Saturday into the far distant future. You might also like to look at the work of the Chudnovsky brothers.
On it!
Having been a fan of Srinivasa Ramanujan since my youth, I thought I should
give Ramanujanʼs 1/ series a try.
As can be seen on the linked Wikipedia page, Ramanujanʼs 1914 formula looks quite involved, with (4·k)!÷((k!)⁴·396⁴ᵏ) being the most computationally intensive part.
Compared to a naïve implementation, this part of the formula can fortunately be sped up significantly by relying on the following recurrence relation:
s₀ = 1
sₖ₊₁ = sₖ·(4k+1)(4k+2)(4k+3)(4k+4) ÷ ((k+1)⁴·396⁴)
With that, Ramanujanʼs formula can be simplified to (in LaTeX math notation):
\frac{1}{\pi} = \frac{2\sqrt{2}}{9801} \sum_{k=0}^{\infty} (1103 + 26390k)s_k
Implementing this variant of the equation by relying on Cʼs
doublefloating-point type might then look as follows:static double compute_pi_ramanujan_double (size_t decimal_digits)
{ double sum = 0.0, s_k = 1.0;
double epsilon = pow (10.0, -(double)decimal_digits);
size_t k;
for (k = 0; ; ++k)
{
double term = s_k * (1103ul + 26390ul*k);
if (term <= epsilon)
break;
sum += term;
{ double signature_index = 4.0 * k,
index_plus_1_sq = (k + 1.0) * (k + 1.0);
double fraction =
((signature_index + 1.0) * (signature_index + 2.0)
* (signature_index + 3.0) * (signature_index + 4.0))
/ (index_plus_1_sq * index_plus_1_sq
* 156816ul * 156816ul);
s_k *= fraction;
}
}
return (1.0 / (sum * (2.0*sqrt(2.0) / 9801ul)));
}
Since Ramanujanʼs formula increases the number of digits by ~7.98 digits in each iteration, Cʼs
doubledata type is only suitable for k up to 2.To allow for the calculation of any number of digits of , the following program relies on GNU GMP (under Linux) / MPIR (under Windows) for efficient support of arbitrary-precision arithmetic.
With that, calculating the first 100,000 digits takes ~10.992 seconds on my machine. Calculating the first 1,000,000 digits already takes ~1455.476 seconds, however.
Without further optimizations, such as a divide-and-conquer algorithm like Binary splitting, I therefore wouldnʼt recommend to go any further than this with the linked implementation.
Note: the given program works under Linux and Windows. For Windows, a pre-compiled version of the MPIR library can be found at this link, which can then be used with VS2008, VS2022 and VS2026.
Wow
Thanks for the “wow”… but to be honest, I donʼt really understand the underlying mathematics. Thatʼs well above my pay grade. But regardless of how all the constants in the given formula can be justified, it seems to work out in the end. Thatʼs good enough for me ☺
(If I can find time to do so, I may even try to implement the mentioned divide-and-conquer strategy. I donʼt know yet… weʼll see.)
I admire your candor. I don’t understand the math thanks to a poor education. Still, I’m delighted when I can code something using this type of math and it works!