# RE: nth root

Discussion in 'Python' started by Tim Roberts, Jan 31, 2009.

1. ### Tim RobertsGuest

Dan,

Thanks - you're probably right - just my intuition said to me that rather than calculating that the 13th root of 4021503534212915433093809093996098953996019232
is 3221.2904208350265....
there must be a quicker way of finding out its between 3221 and 3222....

.....but perhaps not.

Tim

________________________________

From: Dan Goodman [mailto:]
Sent: Sat 31-Jan-09 3:11 PM
To:
Subject: Re: nth root

Takes less than 1 sec here to do (10**100)**(1./13) a million times, and
only about half as long to do (1e100)**(1./13), or about 14 times as
long as to do .2**2. Doesn't look like one could hope for it to be that
much quicker as you need 9 sig figs of accuracy to get the integer part

Dan

Tim wrote:
> In PythonWin I'm running a program to find the 13th root (say) of
> millions of hundred-digit numbers. I'm using
> n = 13
> root = base**(1.0/n)
> which correctly computes the root to a large number of decimal
> places, but therefore takes a long time. All I need is the integer
> component. Is there a quicker way?
>
>
>
> ------------------------------------------------------------------------
>
> --
> http://mail.python.org/mailman/listinfo/python-list

--
http://mail.python.org/mailman/listinfo/python-list

Tim Roberts, Jan 31, 2009

2. ### Mark DickinsonGuest

On Jan 31, 5:43 am, "Tim Roberts" <> wrote:
> Dan,
>
> Thanks - you're probably right - just my intuition said to me that rather than calculating that the 13th root of 4021503534212915433093809093996098953996019232
> is 3221.2904208350265....
> there must be a quicker way of finding out its between 3221 and 3222....
>
> ....but perhaps not.

I don't think you'll find anything much quicker than n**(1./13)
(though I hope
that if you're doing this millions of time then you're precomputing
the 1./13
rather than redoing the division every single time.

What happens behind the scenes here is that your integer is
immediately
converted to a float, then the system math library is used for the
power operation. The integer -> float conversion is probably quite
significant, timewise.

I'd also be a bit worried about accuracy. Is it important to you
that the
integer part of the result is *exactly* right, or is it okay if
(n**13)**(1./13) sometimes comes out as slightly less than n, or if
(n**13-1)**(1./13) sometimes comes out as n?

Mark

Mark Dickinson, Jan 31, 2009

3. ### Steve HoldenGuest

Mark Dickinson wrote:
> On Jan 31, 5:43 am, "Tim Roberts" <> wrote:
>> Dan,
>>
>> Thanks - you're probably right - just my intuition said to me that rather than calculating that the 13th root of 4021503534212915433093809093996098953996019232
>> is 3221.2904208350265....
>> there must be a quicker way of finding out its between 3221 and 3222....
>>
>> ....but perhaps not.

>
> I don't think you'll find anything much quicker than n**(1./13)
> (though I hope
> that if you're doing this millions of time then you're precomputing
> the 1./13
> rather than redoing the division every single time.
>

Compared with the computation involved in the power computation I think
you'll find this makes a negligible difference in timing. But that's
just mu gut instinct, and we both know that a benchmark is the only way
to be certain, right? It just seems like a possibly premature
optimization to me. [sigh. I had to start this, didn't i?]

>>> t1 = timeit.Timer("x =

4021503534212915433093809093996098953996019232**(1.0/13)")
>>> t2 = timeit.Timer("x =

4021503534212915433093809093996098953996019232**power", "power=1.0/13")
>>> t1.timeit()

1.4860000610351562
>>> t2.timeit()

1.3789999485015869
>>>

Hmm, well, I suppose an 9% speed gain might be worth it.

> What happens behind the scenes here is that your integer is
> immediately
> converted to a float, then the system math library is used for the
> power operation. The integer -> float conversion is probably quite
> significant, timewise.
>

I bow to your superior intuition!

> I'd also be a bit worried about accuracy. Is it important to you
> that the
> integer part of the result is *exactly* right, or is it okay if
> (n**13)**(1./13) sometimes comes out as slightly less than n, or if
> (n**13-1)**(1./13) sometimes comes out as n?
>

Much more significant points, given the limited precision of the doubles
Python will be using. Could gmpy do this better, I wonder?

regards
Steve
--
Steve Holden +1 571 484 6266 +1 800 494 3119
Holden Web LLC http://www.holdenweb.com/

Steve Holden, Jan 31, 2009
4. ### Mark DickinsonGuest

On Jan 31, 1:23 pm, Steve Holden <> wrote:
> [Mark]
> > power operation.  The integer -> float conversion is probably quite
> > significant, timewise.

>
> I bow to your superior intuition!

Here's another timing that shows the significance of the int -> float
conversion: (non-debug build of the trunk)

>>> t1 = timeit.Timer("x = n**power", "n = 4021503534212915433093809093996098953996019232; power = 1./13")
>>> t2 = timeit.Timer("x = n**power", "n = 4021503534212915433093809093996098953996019232.; power = 1./13")
>>> t1.timeit()

0.34778499603271484
>>> t2.timeit()

0.26025009155273438

I've got a patch posted to the tracker somewhere that improves
the accuracy of long->float conversions, while also speeding them
up a tiny bit (but not a lot...).

Mark

Mark Dickinson, Jan 31, 2009
5. ### Mark DickinsonGuest

On Jan 31, 1:23 pm, Steve Holden <> wrote:
> Much more significant points, given the limited precision of the doubles
> Python will be using. Could gmpy do this better, I wonder?

Almost certainly, if exact results are wanted! At least, GMP has
an mpz_root function; I don't know offhand whether gmpy makes it
accessible from Python.

Mark

Mark Dickinson, Jan 31, 2009
6. ### Dan GoodmanGuest

Mark Dickinson wrote:
> I'd also be a bit worried about accuracy. Is it important to you
> that the
> integer part of the result is *exactly* right, or is it okay if
> (n**13)**(1./13) sometimes comes out as slightly less than n, or if
> (n**13-1)**(1./13) sometimes comes out as n?

I don't think accuracy is too big a problem here actually (at least for
13th roots). I just tested it with several hundred thousand random 100
digit numbers and it never made a mistake. The precision of double ought
to easily guarantee a correct result. If you let x=int(1e100**(1./13))
then ((x+1)**13-x**13)/x**13=2.6e-7 so you only need around the first 8
or 9 digits of the 100 digit number to compute the 13th root exactly
(well within the accuracy of a double).

OTOH, suppose you were doing cube roots instead then you would need the
first 35 digits of the 100 digit number and this is more accurate than a
double. So for example int(1e100**(1./3)) is a long way from being the
integer part of the true cube root (it's between 10**18 and 10**19 away).

Dan

Dan Goodman, Jan 31, 2009
7. ### MensanatorGuest

On Jan 31, 8:05 am, Mark Dickinson <> wrote:
> On Jan 31, 1:23 pm, Steve Holden <> wrote:
>
> > Much more significant points, given the limited precision of the doubles
> > Python will be using. Could gmpy do this better, I wonder?

>
> Almost certainly, if exact results are wanted!  At least, GMP has
> an mpz_root function; I don't know offhand whether gmpy makes it
> accessible from Python.
>
> Mark

What am I doing wrong here?

IDLE 2.6b1
>>> import timeit
>>> from gmpy import root
>>> root(4021503534212915433093809093996098953996019232,13)

(mpz(3221), 0)
>>> t1 = timeit.Timer("x = root(a,r)", "a = 4021503534212915433093809093996098953996019232; r = 13")
>>> t1.timeit()

Traceback (most recent call last):
File "<pyshell#5>", line 1, in <module>
t1.timeit()
File "C:\Python26\lib\timeit.py", line 193, in timeit
timing = self.inner(it, self.timer)
File "<timeit-src>", line 6, in inner
NameError: global name 'root' is not defined

Mensanator, Jan 31, 2009
8. ### MensanatorGuest

On Jan 31, 10:53 am, Mensanator <> wrote:
> On Jan 31, 8:05 am, Mark Dickinson <> wrote:
>
> > On Jan 31, 1:23 pm, Steve Holden <> wrote:

>
> > > Much more significant points, given the limited precision of the doubles
> > > Python will be using. Could gmpy do this better, I wonder?

>
> > Almost certainly, if exact results are wanted!  At least, GMP has
> > an mpz_root function; I don't know offhand whether gmpy makes it
> > accessible from Python.

>
> > Mark

>
> What am I doing wrong here?
>
> IDLE 2.6b1
>
> >>> import timeit
> >>> from gmpy import root
> >>> root(4021503534212915433093809093996098953996019232,13)

> (mpz(3221), 0)
> >>> t1 = timeit.Timer("x = root(a,r)", "a = 4021503534212915433093809093996098953996019232; r = 13")
> >>> t1.timeit()

>
> Traceback (most recent call last):
>   File "<pyshell#5>", line 1, in <module>
>     t1.timeit()
>   File "C:\Python26\lib\timeit.py", line 193, in timeit
>     timing = self.inner(it, self.timer)
>   File "<timeit-src>", line 6, in inner
> NameError: global name 'root' is not defined

Never mind, I figured it out.

>>> import gmpy
>>> a=gmpy.mpz(4021503534212915433093809093996098953996019232)
>>> r=13
>>> t1 = timeit.Timer("x = root(a,r)", "from gmpy import root; from __main__ import a,r")
>>> t1.timeit()

4.7018698921850728

For comparison:

>>> t2 = timeit.Timer("x = n**power", "n = 4021503534212915433093809093996098953996019232.; power = 1./13")
>>> t2.timeit()

0.43993394115364026

Mensanator, Jan 31, 2009
9. ### Mark DickinsonGuest

On Jan 31, 4:48 pm, Dan Goodman <> wrote:
> I don't think accuracy is too big a problem here actually (at least for
> 13th roots). I just tested it with several hundred thousand random 100
> digit numbers and it never made a mistake.

Well, random numbers is one thing. But how about the following:

>>> n = 12345**13
>>> n

154662214940914131102165197707101295849230845947265625L
>>> int(n ** (1./13)) # should be 12345; okay

12345
>>> int((n-1) ** (1./13)) # should be 12344; oops!

12345

Mark

Mark Dickinson, Jan 31, 2009
10. ### Dan GoodmanGuest

Mark Dickinson wrote:
> Well, random numbers is one thing. But how about the following:
>
>>>> n = 12345**13
>>>> n

> 154662214940914131102165197707101295849230845947265625L
>>>> int(n ** (1./13)) # should be 12345; okay

> 12345
>>>> int((n-1) ** (1./13)) # should be 12344; oops!

> 12345

Good point! Oops indeed.

Dan

Dan Goodman, Jan 31, 2009
11. ### ajaksuGuest

On Jan 31, 12:03 pm, Mark Dickinson <> wrote:
[...]
> t1 = timeit.Timer("x = n**power", "n = 4021503534212915433093809093996098953996019232; power = 1./13")
> t2 = timeit.Timer("x = n**power", "n = 4021503534212915433093809093996098953996019232.; power = 1./13")

And by using a float literal instead of "float
(402150353421291543309...)", (BTW, here -^), it not only is faster but
also a great way make some innocent bystander waste his eyesight
trying to figure out the magic trick

ajaksu, Jan 31, 2009
12. ### Mark DickinsonGuest

On Jan 31, 7:04 pm, ajaksu <> wrote:
> also a great way make some innocent bystander waste his eyesight
> trying to figure out the magic trick

Oh, come on! At least I put the two lines next to each other!

Mark

Mark Dickinson, Jan 31, 2009
13. ### Tim RobertsGuest

"Tim Roberts" <> wrote:
>
>Thanks - you're probably right - just my intuition said to me that rather than calculating that the 13th root of 4021503534212915433093809093996098953996019232
>is 3221.2904208350265....
>there must be a quicker way of finding out its between 3221 and 3222....
>
>....but perhaps not.

Also, remember that the number you computed there is not really the 13th
root of 4021503534212915433093809093996098953996019232. When you convert
it to float to do the exponentiation, you're losing everything after the
15th significant digit. Of course, if all you're looking for is the
nearest integer, that's not very relevant.

Do you really need the absolute number? You could take log(x)/13 and work
with the log of the results. I suspect (without trying it) that's faster
than the exponentiation.

From one Tim Roberts to another.
--
Tim Roberts,
Providenza & Boekelheide, Inc.

Tim Roberts, Feb 1, 2009