Created
December 20, 2012 04:55
-
-
Save jdpage/4343048 to your computer and use it in GitHub Desktop.
Better way of getting fractional pi approximations. See the corresponding blog post at http://sleepingcyb.org/2012/12/20/approximating-pi-redux for an explanation.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
from __future__ import division | |
from fractions import gcd | |
from itertools import izip, tee, izip_longest | |
import sys | |
def pi_noncanon_seq(): | |
""" | |
produces a generalized continued fraction from of pi, starting with | |
(0 4) (1 1) (3 4) (5 9) ... | |
""" | |
yield (0, 4) | |
k = 1 | |
while True: | |
yield (2*k - 1, k**2) | |
k += 1 | |
def canonicalize(seq): | |
""" | |
takes a generalized continued fraction and converts it to a simple continued | |
fraction. | |
""" | |
a, b, c, d = 1, 0, 0, 1 | |
for p, q in seq: | |
a, b, c, d = p*a+b, q*a, p*c+d, q*c | |
if c != 0 and d != 0 and a//c == b//d: | |
t = a//c | |
a, b, c, d = c, d, a-t*c, b-t*d | |
yield t | |
def pi_seq(): | |
""" returns pi as a simple continued fraction """ | |
return canonicalize(pi_noncanon_seq()) | |
def condense(seq): | |
""" takes a finite sequence and folds it into a rational number. """ | |
a, b, c, d = 1, 0, 0, 1 | |
for p in seq: | |
a, b, c, d = p*a+b, a, p*c+d, c | |
return (a, c) | |
def seq_cmp(seqa, seqb): | |
""" compares two sequences """ | |
k = 1 | |
for aa, bb in izip_longest(seqa, seqb, fillvalue=None): | |
if aa == None and bb == None: | |
return 0 | |
if aa == None: | |
return k | |
if bb == None: | |
return -k | |
if aa != bb: | |
return cmp(aa, bb) * k | |
k = -k | |
def best_rational_approximations(seq): | |
""" | |
takes a generator and outputs best rational approximations until it runs out | |
of terms. | |
""" | |
seq_part = [next(seq)] # integer part | |
yield (seq_part[0], 1) | |
while True: | |
seq, seq_c = tee(seq) | |
try: | |
fin = next(seq) | |
except StopIteration: | |
break | |
fin_min = fin - (fin // 2) | |
if fin % 2 == 0: | |
if seq_cmp(reversed(seq_part + [fin]), seq_c) <= 0: | |
fin_min += 1 | |
for x in range(fin_min, fin + 1): | |
yield condense(seq_part + [x]) | |
seq_part.append(fin) | |
def head(n, seq): | |
""" generates the first n elements of seq, then stops """ | |
for _, x in izip(range(n), seq): | |
yield x | |
if __name__ == "__main__": | |
lim = int(sys.argv[1]) | |
for p, q in head(lim, best_rational_approximations(pi_seq())): | |
print "%30d/%-30d %.15f" % (p, q, p/q) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment