Quote:

> 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

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