Ticket #10720: trac_10720_power_series_nth_root_3.patch

File trac_10720_power_series_nth_root_3.patch, 2.3 KB (added by mario pernici, 12 years ago)
  • sage/rings/power_series_ring_element.pyx

    # HG changeset patch
    # User Mario Pernici <mario.pernici@gmail.com>
    # Date 1296746495 -3600
    # Node ID 185dae69862db184c2c9bc70c1c2e7cb0d071b35
    # Parent  6232d130dc8af54e6a37c6fca6f92bbd4377d65f
    much faster nth_root for n large using __pow_trunc
    
    diff -r 6232d130dc8a -r 185dae69862d sage/rings/power_series_ring_element.pyx
    a b  
    12071207        b_n = 1/n
    12081208        a_n = 1+1/n
    12091209        for next_prec in sage.misc.misc.newton_method_sizes(prec)[1:]:
    1210             z = current**(n+1) * A.truncate(next_prec)
     1210            z = __pow_trunc(current,n+1,next_prec) * A.truncate(next_prec)
    12111211            current = a_n*current - b_n * z.truncate(next_prec)
    12121212
    12131213        S = self._parent
     
    20032003
    20042004def make_element_from_parent_v0(parent, *args):
    20052005    return parent(*args)
     2006
     2007def __pow_trunc(a, n, prec):
     2008    """
     2009    compute a^n with precision prec; helper function for nth_root
     2010
     2011    code adapted from generic_power_c in sage.structure.element.pyx
     2012    """
     2013    if n < 4:
     2014        # These cases will probably be called often
     2015        # and don't benefit from the code below
     2016        if n == 1:
     2017            return a
     2018        elif n == 2:
     2019            return (a*a).truncate(prec)
     2020        elif n == 3:
     2021            return ((a*a).truncate(prec)*a).truncate(prec)
     2022
     2023    # check for idempotence, and store the result otherwise
     2024    aa = a*a
     2025    if aa == a:
     2026        return a
     2027
     2028    # since we've computed a^2, let's start squaring there
     2029    # so, let's keep the least-significant bit around, just
     2030    # in case.
     2031    m = n & 1
     2032    n = n >> 1
     2033
     2034    # One multiplication can be saved by starting with
     2035    # the second-smallest power needed rather than with 1
     2036    # we've already squared a, so let's start there.
     2037    apow = aa
     2038    while n&1 == 0:
     2039        apow = (apow*apow).truncate(prec)
     2040        n = n >> 1
     2041    power = apow
     2042    n = n >> 1
     2043
     2044    # now multiply that least-significant bit in...
     2045    if m:
     2046        power = (power * a).truncate(prec)
     2047
     2048    # and this is straight from the book.
     2049    while n != 0:
     2050        apow = (apow*apow).truncate(prec)
     2051        if n&1 != 0:
     2052            power = (power*apow).truncate(prec)
     2053        n = n >> 1
     2054
     2055    return power
     2056