Item — HN Reader

I wrote a quick script to walk through Q, spitting out fractions with high "quality". Here are some early values:

  22/7: 3.42928833728
  355/113: 3.20195874233
  710/226: 2.79251047296

A quick search through 1/1 - 100,000/100,000 shows no other interesting values. My algorithm is extremely naive though (brute force), a faster search could potentially be used to find other high quality values much more quickly.

Update - A slightly optimized algorithm, and nothing special found with denominators up to the hundreds of millions. I know this means nothing wrt actual mathematics though...

4 comments

http://en.wikipedia.org/wiki/Continued_fraction#Best_rationa... will teach you how to find the best approximations in no time.

If you understand that, read http://blog.wolfram.com/2011/06/30/all-rational-approximatio... :-)

By your output you're performing some redundant checks. 355/113 == 710/226. Not sure how you're performing the search, but removing those may speed things up.

I think the cost of removing them (checking every value to see if it's a multiple of another) is greater than the benefit of removing them. My search terminated though (with no interesting finds) - my fractional values were getting too small and I got a ZeroDivisionError.

You can easily generate all fractions in lowest terms within a disk by using a recursive Farey sequence generator.

Reading about Farey sequences led me to http://en.wikipedia.org/wiki/Stern%E2%80%93Brocot_tree. This could be used to perform a search of the space of rational numbers with denominator < threshold with each number examined only once.

Edit: Oops, suggested an infinitely long search

The Wilf-Calkin tree also does that:

http://en.wikipedia.org/wiki/Calkin%E2%80%93Wilf_tree

I refer to it as Wilf-Calkin because that's how Neil Calkin refers to it, despite the alphabetical ordering on the Wikipedia page.

Even better, there's an explicit formula for generating rationals, without repeats:

http://en.wikipedia.org/wiki/Calkin%E2%80%93Wilf_tree#Breadt...

Neat, that was in a tab I hadn't reached yet. Though for the purposes of this sort of search, the Stern-Brocot tree has the nice property of being a binary search tree. Allowing you to constrain the breadth as well as the depth of the search.

I had a quick go too:

    num     den     difference to pi
    1	1	2.1415926535898
    3	1	0.14159265358979
    22	7	0.0012644892673497
    333	106	8.3219627529107e-05
    355	113	2.6676418940497e-07
    103993	33102	5.7789062424263e-10
    104348	33215	3.3162805834763e-10
    208341	66317	1.2235634727631e-10
    312689	99532	2.914335439641e-11
    833719	265381	8.7152507433075e-12
    1146408	364913	1.6107115641262e-12
    4272943	1360120	4.0412118096356e-13
    5419351	1725033	2.2204460492503e-14
    80143857	25510582	4.4408920985006e-16

Lua script to generate (dirty):

    local max = 100000000
    local c = {}
    local min , j = math.huge
    for i=1,max do
        local a = math.pi-math.abs(math.cos(i))
        if a < min then
            min = a
            j = i
            table.insert ( c , j )
        end
    end

    local min , j = math.huge
    for k=1,#c do
        for i=1,max do
            local a = math.abs(math.pi-c[k]/i)
            if a < min then
                min = a
                j = i
                print(c[k],j,min)
            end
        end
    end

It's not actually the difference to pi that matters, it's the "quality" of the approximation, as described in the article. You can always get better approximations by using larger denominators, it's just that some approximations are, given the size of the denominator, unreasonably good.

355/113 is one of those.

True, I didn't do the quality measure they talk of. I was just fascinated by:

"For example, use any scientific calculator to compute cos(355) in radians. The oddball result is due to the freakish closeness of 355/113 to pi."

And wanted to see what else would fit this criteria...

The numbers that you're producing can be better produced by the following Python code:

    from fractions import Fraction as frac
    def continued_fraction(x, n=10):
        """Represent x as a continued fraction out to n places."""
        last = int(x)
        out = []
        for i in range(n):
            x = 1/(x - last)
            last = int(x)
            out.append(last)
        return out

    def bras(x, n=10):
        """Compute the best rational approximations to x."""
        base = int(x)
        nums = continued_fraction(x, n)
        S = lambda depth, i = 0: \
            frac(1, nums[i]) if depth == 0 else 1 / (nums[i] + S(depth - 1, i + 1))
        return ', '.join(str(x) for x in (base,) + tuple(base + S(k) for k in range(n)))

    from decimal import Decimal
    print(bras(Decimal('3.1415926535897932384626433832795028841971693993')))
    # 3, 22/7, 333/106, 355/113, 103993/33102, 104348/33215, 208341/66317,
    # 312689/99532, 833719/265381, 1146408/364913, 4272943/1360120

These are the best rational approximations to pi, which are determined by the continued fraction of pi:

    pi ~=  3 + 1/(7 + 1/(15 + ...))
    pi = 3 + [7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 2, ...]

Large numbers are opportunities to "grab a lot of extra digits" and thus the large qualities come when we truncate after 7, 15, 292, and 14.

You may wonder what the "most irrational" number is, in the sense that all of its best-rational-approximations are of low quality. That distinction belongs to:

    phi = (1 + Decimal(5).sqrt()) / 2 # golden ratio
    continued_fraction(phi) # [1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
    bras(phi) # 1, 2, 3/2, 5/3, 8/5, 13/8, 21/13, 34/21, 55/34, 89/55, 144/89

Those last numbers you may recognize from your intro to programming as the Fibonaccis.

Here's some simpler code that computes the same result -- I'm too tired and dumb to work out if it will for any number:

    def rationalizations(x):
        """Generate good rational approximations of x in order of
        increasing denominator."""
        assert 0 <= x
        ix = int(x)
        yield ix, 1
        if x != ix:
            for numer, denom in rationalizations(1.0/(x-ix)):
                yield denom + ix * numer, numer

    import itertools, math
    def show(x, n): return list(itertools.islice(rationalizations(x), n))

    print show(math.pi, 11)
    # [(3, 1), (22, 7), (333, 106), (355, 113), (103993, 33102),
    # (104348, 33215), (208341, 66317), (312689, 99532),
    # (833719, 265381), (1146408, 364913), (4272943, 1360120)]

From https://github.com/darius/sketchbook/blob/master/misc/ration...

It's true, my code is longer because it computes an explicit continued fraction representation (and also because it's not written as two generator scripts).

I think that code must be doing something wrong. 355/113 = 710/226 but you have different values for both