Velocity Reviews > RE: nth root

# RE: nth root

Tim Roberts
Guest
Posts: n/a

 01-31-2009
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 [(E-Mail Removed)]
Sent: Sat 31-Jan-09 3:11 PM
To: http://www.velocityreviews.com/forums/(E-Mail Removed)
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

Mark Dickinson
Guest
Posts: n/a

 01-31-2009
On Jan 31, 5:43*am, "Tim Roberts" <(E-Mail Removed)> 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

Steve Holden
Guest
Posts: n/a

 01-31-2009
Mark Dickinson wrote:
> On Jan 31, 5:43 am, "Tim Roberts" <(E-Mail Removed)> 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**po wer", "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/

Mark Dickinson
Guest
Posts: n/a

 01-31-2009
On Jan 31, 1:23*pm, Steve Holden <(E-Mail Removed)> 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
Guest
Posts: n/a

 01-31-2009
On Jan 31, 1:23*pm, Steve Holden <(E-Mail Removed)> 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

Dan Goodman
Guest
Posts: n/a

 01-31-2009
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

Mensanator
Guest
Posts: n/a

 01-31-2009
On Jan 31, 8:05*am, Mark Dickinson <(E-Mail Removed)> wrote:
> On Jan 31, 1:23*pm, Steve Holden <(E-Mail Removed)> 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(402150353421291543309380909399609895399601923 2,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
Guest
Posts: n/a

 01-31-2009
On Jan 31, 10:53*am, Mensanator <(E-Mail Removed)> wrote:
> On Jan 31, 8:05*am, Mark Dickinson <(E-Mail Removed)> wrote:
>
> > On Jan 31, 1:23*pm, Steve Holden <(E-Mail Removed)> 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(402150353421291543309380909399609895399601923 2,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(402150353421291543309380909399609895399 6019232)
>>> 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

Mark Dickinson
Guest
Posts: n/a

 01-31-2009
On Jan 31, 4:48*pm, Dan Goodman <(E-Mail Removed)> 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

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

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

12345

Mark

Dan Goodman
Guest
Posts: n/a

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

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

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

> 12345

Good point! Oops indeed.

Dan