cube root

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
16 messages Options
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

cube root

col speed
Hi there, just a quick one.
Is there a way to obtain cube roots in python?
I have been trying:
math.pow(n,1.0/3)
This works with some, but with 64, for example I get:
>>> pow(64,1.0/3)
3.9999999999999996
However:
>>> 4**3
64
Any ideas?
Thanks
Colin

_______________________________________________
Tutor maillist  -  [hidden email]
http://mail.python.org/mailman/listinfo/tutor
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: cube root

Brett Wilkins
Hey Colin,

What you're running into here is the limited accuracy of floating point
values...
You'll likely find this happens a lot in python.. CPython, at least. (I
know, as I do)
I'm not sure, as I've never used it... but perhaps Numeric/Numpy handle
this kinda stuff better (for accuracy's sake)


Cheers,
--Brett

col speed wrote:

> Hi there, just a quick one.
> Is there a way to obtain cube roots in python?
> I have been trying:
> math.pow(n,1.0/3)
> This works with some, but with 64, for example I get:
> >>> pow(64,1.0/3)
> 3.9999999999999996
> However:
> >>> 4**3
> 64
> Any ideas?
> Thanks
> Colin
> ------------------------------------------------------------------------
>
> _______________________________________________
> Tutor maillist  -  [hidden email]
> http://mail.python.org/mailman/listinfo/tutor
>  

_______________________________________________
Tutor maillist  -  [hidden email]
http://mail.python.org/mailman/listinfo/tutor
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: cube root

Senthil Kumaran-6
On Mon, Jan 19, 2009 at 12:11 PM, Brett Wilkins <[hidden email]> wrote:
> Hey Colin,
> What you're running into here is the limited accuracy of floating point
> values...

Here the Python Documentation mentioning about the inherent
limitations in dealing with floating point numbers:

http://docs.python.org/tutorial/floatingpoint.html


--
Senthil
_______________________________________________
Tutor maillist  -  [hidden email]
http://mail.python.org/mailman/listinfo/tutor
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: cube root

Alan Gauld
In reply to this post by Brett Wilkins

"Brett Wilkins" <[hidden email]> wrote

> What you're running into here is the limited accuracy of floating
> point
> values... You'll likely find this happens a lot in python.. CPython,
> at least.

In fact it happens with virtually every programming language
available.
The problem traces back to the way that computers represent floating
point numbers in hardware. The integrated circuits cannot hold
infinitely
long numbers so they need to use a condensed representation and
that inherently introduces small inaccuracies. (Think of 1/3 in
decimal - an infinite series of 1.3333..., the same is true in binary,
some numbers just cannot be represented accurately)

In the specific case being discussed you can disguise the problem
quite easily by displaying the result with a smaller number of decimal
places using round() or string formatting:

>>> 64**(1.0/3)
3.9999999999999996
>>> print 64**(1.0/3)
4.0
>>> round(64**(1.0/3))
4.0
>>> "%5.3f" % 64**(1.0/3)
'4.000'
>>>

When you need to compare floating point numbers you also
need to be careful since:

>>> 4.0 == 64**(1.0/3)     ## is false!

You need to introduce an error range, typically:

>>> e = 0.00000001
>>> 4.0-e < 64**(1.0/3) < 4.0+e     #  is true

You can also use the decimal module for exact decimal
arithmetic but it introduces some other complications.

HTH,


--
Alan Gauld
Author of the Learn to Program web site
http://www.freenetpages.co.uk/hp/alan.gauld


_______________________________________________
Tutor maillist  -  [hidden email]
http://mail.python.org/mailman/listinfo/tutor
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: cube root

Vicent-3
In reply to this post by Brett Wilkins




On Mon, Jan 19, 2009 at 07:41, Brett Wilkins <[hidden email]> wrote:

What you're running into here is the limited accuracy of floating point
values...
You'll likely find this happens a lot in python.. CPython, at least. (I
know, as I do)
I'm not sure, as I've never used it... but perhaps Numeric/Numpy handle
this kinda stuff better (for accuracy's sake)


Maybe this issue can be overcome by using symbolic notation when possible (maybe by using SymPy)? Has anyone any experience with this?

--
Vicent

_______________________________________________
Tutor maillist  -  [hidden email]
http://mail.python.org/mailman/listinfo/tutor
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: cube root

Paul McGuire
In reply to this post by col speed
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  -  [hidden email]
http://mail.python.org/mailman/listinfo/tutor
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: cube root

Brett Wilkins
In reply to this post by col speed
The only language I've run into so far (I haven't used many, mind) that
doesn't have this issue is Scheme :)
(Just learning it at the moment.)

Cheers,
--Brett

P.S. Forgive me if this email doesn't sort properly, sending through
webmail, as I don't have a relaying SMTP available to me currently.

Alan Gauld wrote:

>
> "Brett Wilkins" <[hidden email]> wrote
>
>> What you're running into here is the limited accuracy of floating point
>> values... You'll likely find this happens a lot in python.. CPython,
>> at least.
>
> In fact it happens with virtually every programming language available.
> The problem traces back to the way that computers represent floating
> point numbers in hardware. The integrated circuits cannot hold infinitely
> long numbers so they need to use a condensed representation and
> that inherently introduces small inaccuracies. (Think of 1/3 in
> decimal - an infinite series of 1.3333..., the same is true in binary,
> some numbers just cannot be represented accurately)
>
> In the specific case being discussed you can disguise the problem
> quite easily by displaying the result with a smaller number of decimal
> places using round() or string formatting:
>
>>>> 64**(1.0/3)
> 3.9999999999999996
>>>> print 64**(1.0/3)
> 4.0
>>>> round(64**(1.0/3))
> 4.0
>>>> "%5.3f" % 64**(1.0/3)
> '4.000'
>>>>
>
> When you need to compare floating point numbers you also
> need to be careful since:
>
>>>> 4.0 == 64**(1.0/3)     ## is false!
>
> You need to introduce an error range, typically:
>
>>>> e = 0.00000001
>>>> 4.0-e < 64**(1.0/3) < 4.0+e     #  is true
>
> You can also use the decimal module for exact decimal
> arithmetic but it introduces some other complications.
>
> HTH,
>
>



_______________________________________________
Tutor maillist  -  [hidden email]
http://mail.python.org/mailman/listinfo/tutor
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: cube root

Kent Johnson
On Mon, Jan 19, 2009 at 6:11 AM, Brett Wilkins <[hidden email]> wrote:
> The only language I've run into so far (I haven't used many, mind) that
> doesn't have this issue is Scheme :)

It doesn't have an issue with cube roots or with floating point
inaccuracies in general? If the latter, I would like to know how they
do that...

Kent
_______________________________________________
Tutor maillist  -  [hidden email]
http://mail.python.org/mailman/listinfo/tutor
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: cube root

spir
In reply to this post by Alan Gauld
Do you know any common algorithm to convert decimal (in the sense of fractional) decimals (in the sense of base 10 numbers) into binaries?

123.456 --> 1111011.bbb...
and/or
123456 * 10**(-3) --> bbb... * 2**(-bbb...)

How do python/C achieve that?

denis

------
la vida e estranya
_______________________________________________
Tutor maillist  -  [hidden email]
http://mail.python.org/mailman/listinfo/tutor
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: cube root

Andre Engels
In reply to this post by Brett Wilkins
On Mon, Jan 19, 2009 at 12:11 PM, Brett Wilkins <[hidden email]> wrote:
> The only language I've run into so far (I haven't used many, mind) that
> doesn't have this issue is Scheme :)
> (Just learning it at the moment.)

It doesn't? That would surprise me. The only one that I know to do
this kind of thing correctly is Matlab.

--
André Engels, [hidden email]
_______________________________________________
Tutor maillist  -  [hidden email]
http://mail.python.org/mailman/listinfo/tutor
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: cube root

Brett Wilkins
In reply to this post by col speed
Hmm, Well I thought it was both, but the latter seems untrue (now that I
test a bit more)

(expt 64 (/ 1 3)) gives the value 4, but turning any of those into
floating point numbers seems to give me the infamous 3.999999999999996
thing all over again.

I was originally thinking that scheme would handle exponents as floating
point numbers, the way that C does (ie 2e20 is a fp number...), and I
tested 2e20+1-2e20, which stumps most languages (python just gives 0.0,
gcc 4.0.1 gives the same).

Oh well, sorry to get your hopes up!
 Btw, I'm using MIT/GNU scheme, in case that makes a difference.

Cheers
--Brett

Kent Johnson wrote:
> On Mon, Jan 19, 2009 at 6:11 AM, Brett Wilkins <[hidden email]> wrote:
>> The only language I've run into so far (I haven't used many, mind) that
>> doesn't have this issue is Scheme :)
>
> It doesn't have an issue with cube roots or with floating point
> inaccuracies in general? If the latter, I would like to know how they
> do that...
>
> Kent


_______________________________________________
Tutor maillist  -  [hidden email]
http://mail.python.org/mailman/listinfo/tutor
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: cube root

Andre Engels
In reply to this post by spir
On Mon, Jan 19, 2009 at 1:13 PM, spir <[hidden email]> wrote:
> Do you know any common algorithm to convert decimal (in the sense of fractional) decimals (in the sense of base 10 numbers) into binaries?
>
> 123.456                 --> 1111011.bbb...
> and/or
> 123456 * 10**(-3)       --> bbb... * 2**(-bbb...)
>
> How do python/C achieve that?

There's probably more efficient methods, but a simple method to
convert a decimal fraction to a binary would be the following
(untested):

def tobinary(n,precision=12)
    # n is a number between 0 and 1 that should be converted,
precision is the number of binary digits to use.
    assert 0.0 <= n < 1.0
    outcome = "0."
    compare = 0.5
    for i in xrange(precision):
        if n > compare:
            outcome += "1"
            n -= compare
            if n == 0.0: break
        else:
            outcome += "0"
        compare /= 2
    return outcome

--
André Engels, [hidden email]
_______________________________________________
Tutor maillist  -  [hidden email]
http://mail.python.org/mailman/listinfo/tutor
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: cube root

Ken Oliver
In reply to this post by col speed


-----Original Message-----

>From: Andre Engels <[hidden email]>
>Sent: Jan 19, 2009 7:22 AM
>To: spir <[hidden email]>
>Cc: [hidden email]
>Subject: Re: [Tutor] cube root
>
>On Mon, Jan 19, 2009 at 1:13 PM, spir <[hidden email]> wrote:
>> Do you know any common algorithm to convert decimal (in the sense of fractional) decimals (in the sense of base 10 numbers) into binaries?
>>
>> 123.456                 --> 1111011.bbb...
>> and/or
>> 123456 * 10**(-3)       --> bbb... * 2**(-bbb...)
>>
>> How do python/C achieve that?
>
>There's probably more efficient methods, but a simple method to
>convert a decimal fraction to a binary would be the following
>(untested):
>
>def tobinary(n,precision=12)
>    # n is a number between 0 and 1 that should be converted,
>precision is the number of binary digits to use.
>    assert 0.0 <= n < 1.0
>    outcome = "0."
>    compare = 0.5
>    for i in xrange(precision):
>        if n > compare:
>            outcome += "1"
>            n -= compare
>            if n == 0.0: break
>        else:
>            outcome += "0"
>        compare /= 2
>    return outcome
>
>--
>André Engels, [hidden email]

I hope my memory serves.  To convert decimal fraction into binary fraction, do the following repeatedly to the decimal until desired accuracy is achieved.

Multiply by 2, if result is less than 1, next digit is 0 else next digit is 1 and you drop the whole number part of the result. Continue...

.456 * 2 = .912    (first digit is 0)
.912 * 2 = 1.824   (next digit is 1)
.824 * 2 = 1.648   (next digit is 1)
.648 * 2 = 1.296   (next digit is 1)
.296 * 2 = .592    (next digit is 0)

0.456 ~ 0.01110
_______________________________________________
Tutor maillist  -  [hidden email]
http://mail.python.org/mailman/listinfo/tutor
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: cube root

col speed
In reply to this post by Vicent-3
Wow! I seem to have caused a great deal of comments!
I actually am looking to see if a number is a "perfect cube". I will try out all the suggestions.
Thanks a lot.
Colin
P.S.
I have a small programme that changes decimal to binary (I don't know if it can be changed to work with fractions or not), maybe it will help.

def dec2bin(dec):
    bin = ''
    while dec > 0:
        bin = str(dec & 1) + bin
        dec >>= 1
    return bin

2009/1/19 Vicent <[hidden email]>




On Mon, Jan 19, 2009 at 07:41, Brett Wilkins <[hidden email]> wrote:

What you're running into here is the limited accuracy of floating point
values...
You'll likely find this happens a lot in python.. CPython, at least. (I
know, as I do)
I'm not sure, as I've never used it... but perhaps Numeric/Numpy handle
this kinda stuff better (for accuracy's sake)


Maybe this issue can be overcome by using symbolic notation when possible (maybe by using SymPy)? Has anyone any experience with this?

--
Vicent


_______________________________________________
Tutor maillist  -  [hidden email]
http://mail.python.org/mailman/listinfo/tutor
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: cube root

Chris Fuller-5
On Monday 19 January 2009 18:56, col speed wrote:
> Wow! I seem to have caused a great deal of comments!
> I actually am looking to see if a number is a "perfect cube". I will try
> out all the suggestions.
> Thanks a lot.
> Colin

The reliable way to do this is to store a list of cubes.  If the number you
want to check is not in your list, then compute more cubes until you exceed
that number.  Unless you are dealing with absurdly huge numbers, this won't
be very onerous on computing time.

Cheers
_______________________________________________
Tutor maillist  -  [hidden email]
http://mail.python.org/mailman/listinfo/tutor
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: cube root

XXChris123
This post has NOT been accepted by the mailing list yet.
In reply to this post by col speed
# Find the  cube root of a  perfect  cube
#
# exhaustive enumeration
# brute  force for fun try cube root of 1957816251
x = int(raw_input('Enter an interger:'))
ans = 0
# while ans*ans*ans < abs(x):
while ans**3 < abs(x):

    ans = ans + 1
    print 'current guess =' , ans

print ' last guess = ' , ans
print ' ans*ans*ans = ' , ans*ans*ans
if ans*ans*ans != abs(x):
        print x , 'is not  a perfect  cube'
elif x < 0:
    print 'You entered a negative number'
    ans = - ans
print 'cube root of '+ str(x) + ' is ' + str(ans)



####if you run this program it ask s for input
Loading...