Find closest approximation of Pi for given range of denominators

Question

I am looking for an algorithm (preferably in Go or C) to find the closest common fraction n/d to Pi out of an inclusive range of possible denominators (dmin,dmax with 1 <= dmin <= dmax <= 1e15). If there are multiple common fractions with the same distance to Pi, i want to find the one with the smallest denominator.

Note: A bruteforce approach is not efficient enough, therefore i am looking for a smarter / more efficient solution.

Example: For dmin=1 and dmax=10 the closest common fraction is 22/7 with a distance to Pi of approx 0.001

First thought: Looking at the farey sequence, we could find the closest approximation for all denominators up to dmax. Unfortunately that result does not fulfill the constraint of dmin.

Answer 1

I don't have time for a full answer, but here is a partial answer. This technique uses the concepts of continued fractions--there is much about them online. I'll ignore your value dmin, which is not used below.

Get the continued fraction expansion of pi to as many places as you need. For your bound of dmax <= 1e15 you need only the first 28 numbers, which are

[3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 2, 1, 84, 2, 1, 1, 15, 3, 13]

Use a short loop to find the convergents for pi that have denominators just below and just above dmax. In Python that would be

pi_cont_frac = [3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 
                3, 1, 14, 2, 1, 1, 2, 2, 2, 2,
                1, 84, 2, 1, 1, 15, 3, 13]
denomlo, denomhi = 1, 0
numlo, numhi = 0, 1
for q in pi_cont_frac:
    denomlo, denomhi = denomhi, q * denomhi + denomlo
    numlo, numhi = numhi, q * numhi + numlo
    if denomhi > dmax:
        break

Some software, such as Microsoft Excel, would use the fraction numlo/denomlo, but there may be better approximation than that. Now find the value of natural number r that makes denomhi - r * denomlo just below (or equal to) dmax.

Then either numlo/denomlo or (denomhi - r * denomlo)/(denomhi - r * denomlo) is your desired closest fraction to pi. Just check which one is closer.

This algorithm is of order log(dmax), and due to the properties of pi it is usually much lower. For dmax <= 1e15 it takes 28 loops but a few more clean-up statements.

You could make a faster algorithm by pre-calculating and storing the convergents (values of numhi and denomhi) and doing a search for the value of denomhi just above dmax. This also only takes 28 numbers, but you would need this for both the numerators and the denominators. A binary search would take at most 5 steps to find it--practically instantaneous. Yet another possibility using more storage and less calculation would be to store all intermediate fractions. That storage would go into the hundreds, at least three hundred. If you don't like that stored list for the continued fraction expansion of pi, you could use the value of pi to calculate that on the fly, but using double precision (in C) would get you only to the 28 numbers I showed you.

For more research, look up continued fractions and intermediate fractions.

Answer 2

The more general question is how to represent a fraction p/q as n/d where d is in a desired range drange = [dmin, dmax], q not in drange, gcd(p,q) = 1. I imagine the following refined conditions for the problem:

1. gcd(n, d) = 1 (so 44/14 is not acceptable for a solution to approximating pi as a fraction with a denominator between 10 and 15; the (brute force) answer is 41/13)
2. dmin > q this doesn't apply for the case of pi since the denominator is arbitrarily long but it could apply to the case of finding the nearest fraction to 30/43 with a denominator between 160 and 170; brute force answer is 118/169.
3. q > dmax and Daulton's method succeeds because either it a) serendipitously gives d in the desired range or b) succeeds unconditionally if 2*(dmin - 1) < dmax since answer a/b = m*a/m*b = n/d which is in range.

So the hard cases are

1. when q > dmax but Daulton method either doesn't give an answer that satisfies the conditions because d is less than dmin
2. q < dmax where Daulton's method does not apply

In either case we can see why this will be hard to resolve by thinking about the Stern-Brocot table of fractions. Daulton's method will take us to some denominator smaller than desired, but the answer lies deeper down in the table and we have no way to know which path will take us to the closest representation of p/q that we are seeking.

We can do better than checking all possibilities, however.

Consider that the fraction is 0 < x < 1. Define a method to tell the multiple of n closest to a fraction:

def nfrac(x, n, reduce=1):
    """Return the multiple of 1/n that is closest to x
    which is not necessarily the fraction closest to x
    with denom less than n;
    """
    from math import gcd
    scaled = int(round(x * n))
    m, r = divmod(scaled, n)
    if reduce:
        g = gcd(r, n)
        r, n = r//g, n//g
    return m*n+r, n

If we calculate the nearest fractions with denominators of dmin and dmax we will see the range of numerators that we must consider. These will be fewer than the number of denominators and thus will require less work. This is put together using SymPy

def inrange(x, a, b):
    w, r = divmod(a, x)
    if not r:
        return True
    if a < (w + 1)*x <= b:
        return True
    return False

def dlimit(f, n, m):
    if f > 1:
        w = int(f)
        a, b, c = dlimit(f - w, n, m)
        return a + b*w, b, c
    p, q = f.numerator, f.denominator
    if q > m:
        approx = f.limit_denominator(m)
        p = approx.numerator
        q = approx.denominator
    if inrange(q, n, m):
        w, r = divmod(n, q)
        if r:
            w += 1
        return w*p, w*q, 'easy'
    N,_ = nfrac(f, n, 0)
    M,_ = nfrac(f, m, 0)
    af = 1/f
    p = [nfrac(af, i, 0) for i in range(N, M + 1)]
    if p[0][0] < n:
        p = p[1:]
    if p[-1][0] > m:
        p = p[:-1]
    assert all(n <= i[0] <= m for i in p)
    return tuple(reversed(sorted([i for i in p if i[0]],
        key=lambda x: abs(f - x[1]/x[0]))[0])) +( (N, M),)

>>> import math
>>> from fractions import Fraction
>>> dlimit(Fraction(math.pi),160000,161700)
(507895, 161668, (22655, 22896))

So instead of considering the nearest fraction to 1700 different denominators, the nearest fraction to numerators in range [22655, 22896] -- 241 of them -- are considered instead. This result was checked against the brute-force method and found to be the same. It is not the same as what would be obtained from the method suggested by Daulton:

>>> Fraction(math.pi).limit_denominator(161700)
Fraction(312689, 99532)