RE: nth root

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

  1. Tim Roberts

    Tim Roberts Guest

    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
    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


    --
    http://mail.python.org/mailman/listinfo/python-list
    Tim Roberts, Jan 31, 2009
    #1
    1. Advertising

  2. 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
    #2
    1. Advertising

  3. Tim Roberts

    Steve Holden Guest

    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
    #3
  4. 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
    #4
  5. 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
    #5
  6. Tim Roberts

    Dan Goodman Guest

    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
    #6
  7. Tim Roberts

    Mensanator Guest

    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
    #7
  8. Tim Roberts

    Mensanator Guest

    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
    #8
  9. 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
    #9
  10. Tim Roberts

    Dan Goodman Guest

    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
    #10
  11. Tim Roberts

    ajaksu Guest

    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 :D
    ajaksu, Jan 31, 2009
    #11
  12. 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 :D


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

    Mark
    Mark Dickinson, Jan 31, 2009
    #12
  13. Tim Roberts

    Tim Roberts Guest

    "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
    #13
    1. Advertising

Want to reply to this thread or ask your own question?

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. Replies:
    0
    Views:
    1,199
  2. Code4u
    Replies:
    4
    Views:
    2,624
    Stephen Howe
    Jul 13, 2005
  3. Protoman
    Replies:
    14
    Views:
    1,741
    Alf P. Steinbach
    Sep 17, 2005
  4. Tim Roberts

    RE: nth root

    Tim Roberts, Feb 1, 2009, in forum: Python
    Replies:
    10
    Views:
    602
    Mensanator
    Feb 2, 2009
  5. Peter Laurens

    Nth Root in Ruby

    Peter Laurens, Oct 22, 2007, in forum: Ruby
    Replies:
    3
    Views:
    123
    Peter Laurens
    Oct 22, 2007
Loading...

Share This Page