Wow!  Everybody jumped on the "floating point inaccuracy" topic, I'm
surprised no one tried to find out what the OP was trying to do with his
cube root solver in the first place.  Of course, the first-cut approach to
solving the cube root is to raise to the 1/3 power, but this is not the only
possible approach.  Assuming that the OP wanted to try to find cube roots
for integers, or raise an exception if the given number is not a perfect
cube, I came up with this:

def cube_root(n):
    "A modified Newton's Method solver for integral cube roots."
    if int(n) != n:
        raise ValueError("must provide an integer")
    if n in (-1,0,1): return n
    offset = (1,-1)[n > 0]
    x = n/2
    seen = set()
    steps = 0
    while 1:
        y = x**3
        if y == n:
            #~ print "Found %d ^ 1/3 = %d in %d steps" % (n,x,steps)
            return x
        dydx = 3.0*x*x
        x += int((n-y)/dydx)+offset
        x = x or 1
        if x in seen:
            raise ValueError("%d is not a perfect cube" % n)
        seen.add(x)
        steps += 1

This will correctly give 4 for 64, 5 for 125, etc., not 3.99999999999.

Then I went googling for "python cube root", and found this curious
solution:

def root3rd(x):
    y, y1 = None, 2
    while y!=y1:
        y = y1
        y3 = y**3
        d = (2*y3+x)
        y1 = (y*(y3+2*x)+d//2)//d
    return y

This only works for perfect cubes, but could be made more robust with this
test just before returning y:

    if y*y*y != x:
        raise ValueError("%d is not a perfect cube" % x)

I was intrigued at this solution - timing it versus mine showed it to be 5x
faster.  I tried to see some sort of Newton's method implementation, but
instead its more of a binary search.  This was confirmed by adding this line
inside the while loop:
        print y,y1

And to solve for the cube root of 1860867 (123**3), I get this output:
None 2
2 4
4 8
8 16
16 32
32 62
62 105
105 123

with the value 123 returned as the answer.  Note that the guess (variable
y1) is initially doubled and redoubled, until the solution is neared.  I
tried initializing y1 to a different guess, the original number, and got
this:

None 1860867
1860867 930434
930434 465217
465217 232609
232609 116305
116305 58153
58153 29077
29077 14539
14539 7270
7270 3635
3635 1818
1818 909
909 456
456 235
235 141
141 123

Now the guess is halved each time until nearing the solution, then again
converging to the solution 123.

This algorithm is also robust in that it can do more than just cube roots.
By changing this line:
    y3 = y**3
to:
    y3 = y**4
We get a 4th root solver!  Unfortunately, this breaks down at the 5th root,
as I found solutions would not converge.  At this point, I lost interest in
the general applicability of this algorithm, but it is still an interesting
approach - does anyone know the theoretical basis for it?

I guess I was just disappointed that taking cube roots of so small a number
as 64 quickly got us into floating-point roundoff errors, and if the OP is
doing something like looking for perfect cubes, then there are alternatives
to the one-size-fits-almost-all logarithmic implementation of the '**'
operator.

-- Paul

_______________________________________________
Tutor maillist  -  Tutor@python.org
http://mail.python.org/mailman/listinfo/tutor

Reply via email to