Integer Fraction Approximation of Pi

One day when I was in high school, I forgot my calculator for an Algebra test. I didn't really need it, except for one question where I had to use π. Without my calculator I used the fraction 22/7, which is a good approximation, but my answer was slightly off, so points were deducted (unfairly, in my opinion).

It's been [redacted] years since high school, but for some reason I started thinking about that again. 22/7 isn't a very precise approximation of π, in fact it differs by 0.001264489, but maybe larger integer numbers could give a more precise approximation. I thought I'd play around with some code to see if I could do better than 22/7.

Full disclosure: the code here works, but I do not understand what turns out to be a periodic relationship between locally optimal numerators/denominators. In the end of this post, I veer into a bit of phenomenology, because I think it's cool, but I don't really understand it.

Define the Problem:

To describe the question more precisely, I wanted to know what two integers, when divided, would give the most precise approximation of π that I could use on an imaginary math test where I forgot my calculator.

Preliminary Considerations:

1. The approximation would likely become more accurate with larger numerator and denominator values, so there is not likely a "best" approximation.

2. I would just need to  I needed a truth value of π. The first 100 digits after the decimal are:

3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679

Perfect. Unfortunately, Python's `float` class only supports 15 decimal places, so this value has to be stored as a decimal object.

Here's my code for calculating the best fractional approximation of π:

from math import floor

import decimal

dec_pi = decimal.Decimal('3.141592653589793238462643383279502884'
                         '19716939937510582097494459230781640628'
                         '62089986280348253421170679')

def approx_pi(max_denom=100):

    final_diff = None

    fin_numer = None

    fin_denom = max_denom

    for denom in range(1, max_denom + 1):
        numer = floor(denom * dec_pi)

        for num in (numer, numer + 1):
            abs_diff = abs(dec_pi - (decimal.Decimal(num) / decimal.Decimal(denom)))

            if abs_diff < final_diff:
                final_diff = abs_diff
                fin_numer = num
                fin_denom = denom

    return f'{fin_num}/{fin_denom}'

Results:

Max Denominator	Best Numerator	Best Denominator	Decimal Value
30	22	7	3.14285714286
100	311	99	3.14141414141
300	355	113	3.14159292035
1000	355	113	3.14159292035
3000	355	113	3.14159292035
10000	355	113	3.14159292035
30000	94053	29938	3.14159262476
100000	312689	99532	3.14159265362
300000	833719	265381	3.14159265358
1000000	3126535	995207	3.14159265359

The fraction 355/113 is called Milu, discovered by Chinese astronomer Zu Chongzhi in the 5th century AD. This ratio warrants its own name because a better approximation doesn't appear until the denominator is almost 30,000. I wonder if Zu Chongzhi knew how much larger numbers would have to be to get a better approximation than Milu. Probably. In case you were wondering, Zu had paper but he didn't have Python. The guy also calculated the length of a year to within around 20 seconds of our modern understanding (365.24219878 days). Such a bad ass.

Anyway, to calculate the best fractional approximation of π, I'm analyzing only two numerators for each denominator - that which gives a value just less than π, and that which gives a value just greater than π.  The lower numerator is given by:

numer = floor(denom * dec_pi)

In Big O notation, my code should have a run time of ~O(2n). In other words, the run time for my code scales linearly with the `max_denom` input. This seems to hold, in that it takes ~1.5 seconds per million iterations through the numerator `for` loop. That is to say the code is optimal within reason, and seems to be doing the number of calculations I expect.

Before I made the numerator optimization I described above, I was testing every numerator from 1 to 4x the denominator just to get the code to work, knowing I'd have to optimize later. I happened to mention this to my brother Richard, and he suggested that I use binary search to find the best numerator for any given denominator, or the best denominator for any numerator. I don't think a binary search implementation would have been faster than the code above, but it led me down a weird rabbit hole:

The Part I Don't Understand:

I knew that binary search required an ordered array and just started trying stuff in an effort to look for a strictly increasing relationship of something to something else. Reader, I also appreciate how crazy that sounds, but what I lacked in direction, I more than made up for in enthusiasm.

At some point, I wanted to know if the the best approximation of π strictly improved with larger denominators.  In other words, assuming we have the best integer numerator for any denominator, does the approximation strictly improve as the denominator gets larger? You would think so, but no it doesn't. Here's what I got for denominators up to 1000, wherein the y-axis is the numerical difference between the quotient of the two integers and the true value of π:

I was not expecting such a periodic relationship between denominators and approximation to π. On the bottom of the graph, there are eight data points that stick out among the rest (green box in the plot below). These occur at denominators 113, 226, 339, 452, 565, 678, 791, and 904, which are all integer multiples of 113. But none are better than 355/113!

Furthermore, if you trace what look like overlapping curves in the graph above, points comprising each curve are separated by 7. For example, the seemingly related data points that increase from 233 (see red oval in the graph below) are 233, 240, 247, 254, 261, 268 and 275. Both 113 and 7 are "locally optimal" denominators, in that 22/7 is a closer approximation to pi than can be produced by using 3,4, 5, 6, 8, 9, or 10 as a denominator, and 355/113 is a closer approximation to π than can be produced by using nearby integers as the denominator.

I did the same for numerators. Looking at the relationship between numerators and the best approximation of π, I got a similar graph. Based on plots of denominators vs. deviation of their best approximation to π, I would expect that graphing optimal deviation from π vs. numerators would give me local minima with a period of 355, and that's exactly the case.

Pseudo-curves (for want of a better word) in the region seem to also be separated by 7, but the periodicity of local minima (the eight points at the bottom of the graph) indeed have a periodicity of 355.

For the moment, that's all I have. I don't understand the periodic relationship between numerators or denominators and their best approximations to π, but there it is. π is just a weird number.