Ticket #10720: trac_10720_power_series_nth_root_3.patch

sage/rings/power_series_ring_element.pyx

@@ -1207 +1207 @@
         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
@@ -2003 +2003 @@
 
 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
+