Re: nth root
MRAB
> 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?
>
Have you tried an iterative approach?
def root_13(x):
step = 1
while step ** 13 < x:
step *= 2
root = 0
while step:
if (root + step) ** 13 <= x:
root += step
step //= 2
return root
Probably not that fast in Python, though. numpy or psyco might help.
Tim Roberts
What on earth are numpy and psyco? Do I need to watch the Lord of the Rings?
Tim
Dan Goodman
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
of (10**100)**(1./13) (floats have about 7 and doubles about 16).
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
Tim Roberts
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
________________________________
Dan
--
http://mail.python.org/mailman/listinfo/python-list
Mark Dickinson
> 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
> On Jan 31, 5:43 am, "Tim Roberts" <t.robe...@cqu.edu.au> 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.
>
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/
Mark Dickinson
> [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
> 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
> 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
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
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
Mark Dickinson
> 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
Dan Goodman
> 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
ajaksu
[...]
> 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 :D
Mark Dickinson
> also a great way make some innocent bystander waste his eyesight
> trying to figure out the magic trick :D
Oh, come on! At least I put the two lines next to each other! :-)
Mark
Robert Kern
> Unfortunately, unless I'm doing something wrong, this appears to take 20 times as long... :-)
>
> What on earth are numpy and psyco? Do I need to watch the Lord of the Rings?
No, but Google would help.
--
Robert Kern
"I have come to believe that the whole world is an enigma, a harmless enigma
that is made terrible by our own mad attempt to interpret it as though it had
an underlying truth."
-- Umberto Eco
Tim Roberts
>
>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, ti...@probo.com
Providenza & Boekelheide, Inc.
Tim Roberts
for root in powersp:
nroot = round(bignum**(1.0/root))
if bignum==long(nroot)**root:
.........
which is probably very inefficient, but I can't see anything better.....
Tim
Paul Rubin
First of all, convert the bignum to a float outside of the loop. Second, precompute
the inverses 1.0/2.0, 1.0/3.0, .... Third, do the exponentiation and comparison only
if the float equivalent is very close. I bet almost all the time you're spending is
in bignum to float conversion.
The article
Daniel J. Bernstein (1998). "Detecting perfect powers in essentially
linear time". Mathematics of Computation 67 (223): 1253-1283
(http://cr.yp.to/papers/powers-ams.pdf)
may be of interest.
Paul Rubin
> for root in powersp:
> nroot = round(bignum**(1.0/root))
> if bignum==long(nroot)**root:
> .........
> which is probably very inefficient, but I can't see anything better.....
Actually that will not be accurate enough to know if bignum is, say, a perfect square.
You will have to do the small-power tests at full precision (or use a fancy algorithm).
Larger powers can be done with floats.
Tim Roberts
Yes, very good, on all counts. Many thanks.
Tim
________________________________
From: Paul Rubin [mailto:"http://phr.cx"@NOSPAM.invalid]
Sent: Sun 01-Feb-09 3:53 PM
To: pytho...@python.org
Subject: Re: nth root
The article
may be of interest.
--
http://mail.python.org/mailman/listinfo/python-list
casevh
Take a look at gmpy and the is_power function. I think it will do
exactly what you want.
http://code.google.com/p/gmpy/
casevh
Peter Pearson
You've gotten several promising leads on this thread, but
here's one more thing that might help. Note that if
c == b**13, then c%p == pow( b%p, 13, p ) for any prime p.
So if you have a 100-digit c and a "candidate" 13-th root b,
and you choose p = 101 (for example), a false b has only a 1%
chance of passing the c%p == pow( b%p, 13, p ) test, which
might save some fraction of the cost of whatever more
rigorous test you subsequently apply.
--
To email me, substitute nowhere->spamcop, invalid->net.
Mensanator
And the root function will give you the root AND tell you whether
it was an integral root:
>>> gmpy.root(a,13)
(mpz(3221), 0)
In this case, it wasn't.
>
> http://code.google.com/p/gmpy/
>
> casevh
casevh
exact integral root for any exponent. is_power will give you an answer
to that question but won't tell you what the root or exponent is. Once
you know that the number is a perfect power, you can root to find the
root.
>
> >http://code.google.com/p/gmpy/
>
> > casevh
>
>
Mensanator
casevh
That's the gotcha. :) You still need to test all prime exponents until
you find the correct one. But it is much faster to use is_power to
check whether or not a number has representation as a**b and then try
all the possible exponents than to just try all the possible exponents
on all the numbers.
>
>
>
> > > >http://code.google.com/p/gmpy/
>
> > > > casevh
>
>
Mensanator
Ok. I was under the impression that the OP was only interested
in the 13th root, so that is_power wouldn't necessarily be of any
use to him. If is_power returned true, you would still have to do
root(a,13) to see if it's actually a 13th power. That's executing
two gmpy functions and takes twice as much time as doing root(a,13)
alone.
>
>
>
>
>
> > > > >http://code.google.com/p/gmpy/
>
> > > > > casevh