Wondering about Domingo's Rational Mean book
Kirby Urner
urner at alumni.princeton.edu
Wed May 31 17:45:36 CEST 2000
david_ullrich at my-deja.com wrote:
> There was a long thread right here on comp.lang.python
>(let's see, is that where I am, I think so...) on continued
>fractions recently. A deja search for "continued fraction"
>works.
>
>DU
Thanks David!
I did the search and indeed retrieved the source code I
was looking for, generously provided by Michael Hudson
<mwh21 at cam.ac.uk> in post:
http://www.deja.com/getdoc.xp?AN=624458140&fmt=text
========================================================
epsilon = 2.2204460492503131e-16
def cfrac(x,err=0.0):
result = []
err = err + epsilon
while 1:
a,b = divmod(x, 1.0)
result.append( int(a) )
if not b:
return result
err = err / (b * (b + err))
if err > b:
return result
x = 1 / b
def converge(partials):
p_1, q_1, p, q = 0, 1, 1, 0
result = []
for a in partials:
p_1, q_1, p, q = p,q, a*p+p_1, a*q+q_1
result.append( (p,q) )
return result
========================================================
Usage:
>>> from roots import cfrac,converge
>>> val = 2.**(1./3.)-1
>>> val
0.25992104989487319
>>> converge(cfrac(val))
[(0, 1), (1, 3), (1, 4), (6, 23), (7, 27), (13, 50), (59, 227),
(72, 277), (131, 504), (1120, 4309), (1251, 4813), (18634, 71691),
(19885, 76504), (217484, 836731), (454853, 1749966),
(672337, 2586697),(3144201, 12096754)]
My Python code for a different method, called RP2 at Domingo's
website, gets me the following (numerator, denominator) pairs
(Note: I've stripped off the long integer Ls for readability):
[(1, 5), (5, 19), (19, 73), (73, 281), (281, 1081), (1081, 4159),
(4159, 16001), (16001, 61561), (61561, 236845), (236845, 911219),
(911219, 3505753), (3505753, 13487761), (13487761, 51891761),
(51891761, 199644319), (199644319, 768096001),
(768096001, 2955112721), (2955112721, 11369270485)]
Note how the latter series uses the denominator of the previous
for the numerator of the next. I'll do a quick comparison of
the error (absolute difference from the floating point target
value of 2**(1./3.)-1 in the form of a table:
METHOD 1 (standard CF) METHOD 2 (rational mean)
0.259921049895 0.0599210498949
0.0734122834385 0.00323684484197
0.00992104989487 0.000352922707867
0.000948515322518 0.000134573026546
0.000661790635614 2.34459423146e-005
7.89501051268e-005 2.80031564742e-006
9.15562174542e-006 2.05026694233e-007
6.7479390618e-006 4.01913080594e-009
4.1497423825e-007 4.48542819553e-009
4.54868859245e-008 9.17280640333e-010
2.73094219461e-009 1.23254961792e-010
1.67198754841e-010 1.10377817997e-011
1.51283430228e-011 2.66620059364e-013
4.93438623295e-013 1.35114142097e-013
1.89515070304e-013 3.44169137634e-014
3.14193115969e-014 5.21804821574e-015
5.55111512313e-016 4.99600361081e-016
Note that METHOD 2 (Domingo's RP2) is actually converging
more quickly, right from the start. However, it's also
true that "standard CF" is reaching its target with fewer
digits (shorter numerators and denominators). Let's
compare "total number of digits" for the two methods:
Number of total digits in numerator,denominator:
METHOD1 METHOD2
2 2
2 3
2 4
3 5
3 7
4 8
5 9
5 10
6 11
8 12
8 13
10 15
10 16
12 17
13 18
13 19
15 21
So it's something of a trade off. You'll get more accurate
fractions with fewer iterations using Method2, but you'll
need fewer digits for comparable accuracy using Method1.
Note that my code for Method2 (so far I haven't shared it,
mostly because it's kinda messy) isn't so general in that
it doesn't accept floating point values (like math.pi)
and aim at them. Rather, my RP2 code finds the nth root
of k, with n and k both whole numbers. So at this point
I can't easily give a converging series of fractions for
math.pi, as Michael does in his post.
Kirby
More information about the Python-list
mailing list