Wondering about Domingo's Rational Mean book

[Kirby Urner]

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