Take the 2-minute tour ×
Programming Puzzles & Code Golf Stack Exchange is a question and answer site for programming puzzle enthusiasts and code golfers. It's 100% free, no registration required.

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.

share|improve this question

11 Answers 11

up vote 21 down vote accepted

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
;
share|improve this answer
1  
Well, technically GolfScript has both floating point numbers and constant for PI. It's called "#{Math.PI}". –  xfix May 28 '12 at 17:52
2  
@GlitchMr, in what way is a string a floating point number? –  Peter Taylor Aug 17 '12 at 12:17
    
I'd really like to see this unrolled with comments. –  primo Nov 29 '12 at 11:23
    
Amazing. The very first line 2\!:^2^..292^15.2/3] already blew my mind. –  primo Dec 10 '12 at 15:27
    
@PeterTaylor Tied. Can we do better? –  Eelvex Feb 20 at 6:09

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.

share|improve this answer
    
Pretty cool, but I can't test it because I don't have Mathematica. –  Keith Randall Dec 20 '11 at 17:35
    
@Keith, here's the logic. Round[Pi, x] gives the closest fraction to Pi in steps of x. This is "listable" so Round[Pi,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 pi than the preceding one. –  Mr.Wizard Dec 20 '11 at 18:41
    
Mathematica beating GolfScript. Neat. –  SpellingD Nov 29 '12 at 17:47
    
In 61: Select[Round[f=Pi,1/Range@1*^6],If[#<f,f=#;True]&@Abs[#-Pi]&] ... but useless given the dominant bias –  belisarius Feb 21 at 3:28
    
Yarr, Matie. Thar be magic in this code. –  Michael Stern Feb 21 at 3:31

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.

share|improve this answer
1  
You could change 999999 to 1e6 and save 3 chars. –  M42 Dec 10 '11 at 10:28
    
@M42: Thanks! Down to 82 chars now. –  Ilmari Karonen Dec 10 '11 at 16:05
    
Really nice, $= to get integer. Sorry, I can't upvote twice. –  M42 Dec 11 '11 at 13:36
    
I can't get this to run: String found where operator expected at prog.pl line 1, near "say"$=/$_"" –  Keith Randall Dec 20 '11 at 17:25
    
@KeithRandall: You need the -M5.01 switch (and Perl 5.10.0 or later) for the say command. Sorry for not mentioning that. –  Ilmari Karonen Dec 20 '11 at 22:34

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.

share|improve this answer

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]*?/)}
share|improve this answer
    
@Peter Taylor You're right. You have to use Ruby 1.9. –  Howard Dec 10 '11 at 11:46

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!

share|improve this answer
1  
Some shortenings: 1.0 -> 1., 10**6 -> 1e6 –  Keith Randall Dec 10 '11 at 6:34
    
I have updated with your improvements. Thanks much. –  Steven Rumbalski Dec 10 '11 at 6:40
    
@KeithRandall, but the second of those makes the output violate the spec. –  Peter Taylor Dec 10 '11 at 7:32
    
In second approach there is no need for variable p. That is 4 chars. –  Ante Dec 10 '11 at 14:23
    
@PeterTaylor: I don't understand. How does it violate the spec? –  Steven Rumbalski Dec 11 '11 at 0:48

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.

share|improve this answer

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.

share|improve this answer
    
I can't get this to work: pi.scala:1 error: not found: value math –  Keith Randall Dec 20 '11 at 17:34
    
Are you using Scala 2.8? –  user unknown Dec 21 '11 at 18:49
    
My scala says "unknown version", weird. On ideone.com, they use 2.8.0 and I still get errors. –  Keith Randall Dec 21 '11 at 20:19
    
Try it at simplyscala.com - works for me. For scala-2.8, replacing math with Math might be sufficient. I mentioned simplyscala on this metathread, if you happen to search for it again: meta.codegolf.stackexchange.com/a/401/373 –  user unknown Dec 21 '11 at 20:48
    
OK, that works. –  Keith Randall Dec 22 '11 at 0:52

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]

continued fractions

Please excuse the inconsistent formatting of the continued fractions.

share|improve this answer
    
That's about half-way to a solution, but it's the easiest half. My GolfScript solution hard-codes a suitable representation of the the continued fraction in only 2 characters more. –  Peter Taylor Aug 17 '12 at 12:23
    
But you didn't use continued fractions for your solution to this question, did you? –  David Carraher Aug 17 '12 at 17:21
    
Yes. It was the obvious way to do it. –  Peter Taylor Aug 17 '12 at 19:42
    
In addition to being concise, this is much faster than most or all of the other solutions that have been posted. –  Michael Stern Feb 21 at 3:36

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++;
    }
}
share|improve this answer
2  
var is not always your friend. By removing it in favour of double you gain the ability to merge declarations, lose the requirement to use double literals, and can save 16 chars. OTOH the question asks for a program, so you'll lose a few to adding a class declaration and a Main method. –  Peter Taylor Dec 11 '12 at 9:36

J, 69 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.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.