Good rational approximations of pi

Question

Write a program that prints out all the good rational approximations of pi with denominator < 1000000, in increasing denominator order. a/b is a "good rational approximation" of pi if it is closer to pi than any other rational with denominator no bigger than b.

The output should have 167 lines total, and start and end like this:

3/1
13/4
16/5
19/6
22/7
179/57
...
833719/265381
1146408/364913
3126535/995207

Shortest program wins.

Answer 1

Golfscript, 71 70 69 chars

2\!:^2^..292^15.2/3]{(.)2/.9>+{\+.((}*;.}do;;]-1%{^0@{2$*+\}/"/"\n}/;

(Assumes that you don't pass it anything on stdin)

I don't want to hear any more whinging by people who don't have built-in constants for pi. I don't even have floating point numbers!

See http://en.wikipedia.org/wiki/Continued_fraction#Best_rational_approximations for the background.

# No input, so the stack contains ""
2\!:^2^..292^15.2/3]
# ^ is used to store 1 because that saves a char by allowing the elimination of whitespace
# Otherwise straightforward: stack now contains [2 1 2 1 1 1 292 1 15 7 3]
# Pi as a continued fraction is 3+1/(7+1/(15+1/(...)))
# If you reverse the array now on the stack you get the first 10 continuants followed by 2
# (rather than 3)
# That's a little hack to avoid passing the denominator 1000000

{
    # Stack holds: ... [c_n c_{n-1} ... c_0]
    (.)2/.9>+
    # Stack holds ... [c_{n-1} ... c_0] c_n (1+c_n)/2+((1+c_n)/2 > 9 ? 1 : 0)
    # (1+c_n)/2 > 9 is an ad-hoc approximation of the "half rule"
    # which works in this case but not in general
    # Let k = (1+c_n)/2+((1+c_n)/2 > 9 ? 1 : 0)
    # We execute the next block k times
    {
        # ... [c_{n-1} ... c_0] z
        \+.((
        # ... [z c_{n-1} ... c_0] [c_{n-1} ... c_0] z-1
    }*
    # So we now have ... [c_n c_{n-1} ... c_0] [(c_n)-1 c_{n-1} ... c_0] ...
    #                    [(c_n)-k+1 c_{n-1} ... c_0] [c_{n-1} ... c_0] c_n-k
    ;
    # Go round the loop until the array runs out
    .
}do

# Stack now contains all the solutions as CFs in reverse order, plus two surplus:
# [2 1 2 1 1 1 292 1 15 7 3] [1 2 1 1 1 292 1 15 7 3] ... [6 3] [5 3] [4 3] [3] [2] []
# Ditch the two surplus ones, bundle everything up in an array, and reverse it
;;]-1%

# For each CF...
{
    # Stack holds ... [(c_n)-j c_{n-1} ... c_0]
    # We now need to convert the CF into a rational in canonical form
    # We unwind from the inside out starting with (c_n)-j + 1/infinity,
    # representing infinity as 1/0
    ^0@
    # ... 1 0 [c_n-j c_{n-1} ... c_0]
    # Loop over the terms of the CF
    {
        # ... numerator denominator term-of-CF
        2$*+\
        # ... (term-of-CF * numerator + denominator) numerator
    }/

    # Presentation
    "/"\n
    # ... numerator "/" denominator newline
}/

# Pop that final newline to avoid a trailing blank line which isn't in the spec
;

Answer 2

Mathematica, 67 63

This isn't going to be fast, but I believe it is technically correct.

Round[π,1/Range@1*^6]//.x_:>First/@Split[x,#2≥#&@@Abs[π-{##}]&]

Round[π, x] gives the closest fraction to π in steps of x. This is "listable" so Round[π,1/Range@1*^6] does this for all fractions down to 1/10^6 in order. The resulting list with many "bad" rational approximations is then repeatedly (//.) processed by removing any elements which are farther from π than the preceding one.

Answer 3

Perl, 77 chars

$e=$p=atan2 0,-1;($f=abs$p-($==$p*$_+.5)/$_)<$e&&($e=$f,say"$=/$_")for 1..1e6

A minor challenge is that Perl doesn't have a built-in π constant available, so I first had to calculate it as atan2(0,-1). I'm sure this will be beaten by languages more suited for the job, but it's not bad for a language mainly designed for text processing.

Answer 4

Python, 96 93 89 characters

a=b=d=1.
while b<=1e6:
 e=3.14159265359-a/b;x=abs(e)
 if x<d:print a,b;d=x
 a+=e>0;b+=e<0

Python, 95 93 characters, different algorithm

p=3.14159265359;d=1
for a in range(3,p*1e6):
 b=round(a/p);e=abs(p-a/b)
 if e<d:print a,b;d=e

note: It was less characters to write p=3.14159265359; than from math import*. Darn those wordy imports!

Answer 5

JS (95 characters)

for(i=k=1,m=Math;i<1e6;i++)if((j=m.abs((x=m.round(m.PI*i))/i-m.PI))<k)k=j,console.log(x+'/'+i)

It does print 167 lines.

Answer 6

Ruby 1.9, 84 characters

m=1;(1..1e6).map{|d|n=(d*q=Math::PI).round;k=(n-q*d).abs/d;k<m&&(m=k;puts [n,d]*?/)}

Answer 7

C99, 113 characters

main(d,n){double e=9,p=2*asin(1),c,a=1;for(;n=d*p+.5,c=fabsl(p-a*n/d),d<1e6;++d)c<e&&printf("%d/%d\n",n,d,e=c);}

Need to compile with -lm, and probably full of undefined behaviour, but it works for me.

Answer 8

Scala - 180 chars

import math._
def p(z:Int,n:Int,s:Double):Unit=
if(n==1e6)0 else{val q=1.0*z/n
val x=if(abs(Pi-q)<s){println(z+"/"+n)
abs(Pi-q)}else s
if(Pi-q<0)p(z,n+1,x)else p(z+1,n,x)}
p(3,1,1)

// ungolfed: 457

val pi=math.Pi
@annotation.tailrec
def toPi (zaehler: Int = 3, nenner: Int = 1, sofar: Double=1): Unit = {
  if (nenner == 1000000) () 
  else {
    val quotient = 1.0*zaehler/nenner
    val diff = (pi - quotient)
    val adiff= math.abs (diff)
    val next = if (adiff < sofar) {
      println (zaehler + "/" + nenner) 
      adiff 
    }
    else sofar
    if (diff < 0) toPi (zaehler, nenner + 1, next) 
    else toPi (zaehler + 1, nenner, next) 
  }  
}

The tailrec annotation is just a check, to verify, that it is tail-recursive, which is often a performance improvement.

Answer 9

Mathematica 18 17 chars

I chose to use, as a measure of "best", the number of terms in a continued fraction representation of π. By this criterion, the best rational approximations of π are its convergents.

There are 10 convergents of π with a denominator less than one million. This is fewer than the requested 167 terms, but I am including it here because it may be of interest to others.

Convergents[π, 10] 

(* out *)
{3, 22/7, 333/106, 355/113, 103993/33102, 104348/33215, 208341/66317,
312689/99532, 833719/265381, 1146408/364913}

If you really want to see the denominator for the first convergent, it will cost an additional 11 characters:

Convergents[π, 10] /. {3 -> "3/1"}
(* out *)
{"3/1", 22/7, 333/106, 355/113, 103993/33102, 104348/33215,
208341/66317, 312689/99532, 833719/265381, 1146408/364913}

---

For those who are interested, the following shows the relations among the convergents, partial quotients, and continued fraction expression of convergents of π:

Table[ContinuedFraction[π, k], {k, 10}]
w[frac_] := Row[{Fold[(#1^-1 + #2) &, Last[#], Rest[Reverse[#]]] &[Text@Style[#, Blue, Bold, 14] & /@ ToString /@ ContinuedFraction[frac]]}];
w /@ FromContinuedFraction /@ ContinuedFraction /@ Convergents[π, 10]

Please excuse the inconsistent formatting of the continued fractions.

Answer 10

J, ⁶⁹ 65

New

]`,@.(<&j{.)/({~(i.<./)@j=.|@-l)@(%~(i:3x)+<.@*l=.1p1&)"0>:_i.1e3

Still a brute force approach but much faster and a tad shorter.

Old

A simple "brute force":

(#~({:<<./@}:)\@j)({~(i.<./)@j=.|@-l)@(%~(i:6x)+<.@*l=.1p1&)"0>:i.1e3

make a list of a/bs and then discard those that are farther from π for some b'<b.

Note: Change 1e3 to 1e6 for the full list. Go do something else and return later.

Answer 11

C# 140 129 chars

double n=3,d=1,e=d;while(n<4e5){double w=n/d-Math.PI,a=Math.Abs(w);if(a<e){e=a;Console.WriteLine(n+"/"+d);}if(w>0)d++;else n++;}

Uncompressed code

var numerator = 3d;
var denominator = 1d;
var delta = 4d;
while (numerator < 4e5) 
{
    var newDelta = (numerator / denominator) - Math.PI;
    var absNewDelta = Math.Abs(newDelta);
    if (absNewDelta < delta)
    {
        delta = absNewDelta;
        Console.WriteLine(string.Format("{0}/{1}", numerator, denominator));
    }

    if (newDelta > 0)
    {
        denominator++;
    }
    else
    {
        numerator++;
    }
}