# 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_n = 1/n a_n = 1+1/n for next_prec in sage.misc.misc.newton_method_sizes(prec)[1:]: z = current**(n+1) * A.truncate(next_prec) z = __pow_trunc(current,n+1,next_prec) * A.truncate(next_prec) current = a_n*current - b_n * z.truncate(next_prec) S = self._parent def make_element_from_parent_v0(parent, *args): return parent(*args) def __pow_trunc(a, n, prec): """ compute a^n with precision prec; helper function for nth_root code adapted from generic_power_c in sage.structure.element.pyx """ if n < 4: # These cases will probably be called often # and don't benefit from the code below if n == 1: return a elif n == 2: return (a*a).truncate(prec) elif n == 3: return ((a*a).truncate(prec)*a).truncate(prec) # check for idempotence, and store the result otherwise aa = a*a if aa == a: return a # since we've computed a^2, let's start squaring there # so, let's keep the least-significant bit around, just # in case. m = n & 1 n = n >> 1 # One multiplication can be saved by starting with # the second-smallest power needed rather than with 1 # we've already squared a, so let's start there. apow = aa while n&1 == 0: apow = (apow*apow).truncate(prec) n = n >> 1 power = apow n = n >> 1 # now multiply that least-significant bit in... if m: power = (power * a).truncate(prec) # and this is straight from the book. while n != 0: apow = (apow*apow).truncate(prec) if n&1 != 0: power = (power*apow).truncate(prec) n = n >> 1 return power