Ticket #10720: trac_10720_power_series_nth_root.patch

sage/rings/power_series_ring_element.pyx

@@ -923 +923 @@
         #        b.append(-b[0]*sum([b[n-i]*a[i] for i in range(1,n+1) if i < len(a)]))
         #return self.parent()(b, prec=prec)
     
+
     def valuation_zero_part(self):
         r"""
         Factor self as as `q^n \cdot (a_0 + a_1 q + \cdots)` with
@@ -1098 +1099 @@
             except TypeError:
                 return False
     
+    def nth_root(self, n):
+        """ Returns an `n^{th}` root of self
+
+        INPUT:
+        
+          - n - An integer number
+
+        ALGORITHM: Newton's method to compute $x = y^{\frac{-1}{n}}$
+        
+        .. math::
+        
+           x_{i+1} = (1 + \frac{1}{n}) x_i - \frac{1}{n} y x_{i}^{n+1})
+
+        EXAMPLES::
+        
+        sage: R.<t> = QQ[]
+        sage: p = 1 + t + t^2 + O(t^3)
+        sage: p.nth_root(2)
+        1 + 1/2*t + 3/8*t^2 + O(t^3)
+        sage: p.nth_root(-3)
+        1 - 1/3*t - 1/9*t^2 + O(t^3)
+        sage: p = 1024*t^10*(1 + 2*t + t^2) + O(t^12)
+        sage: p.nth_root(2)                          
+        32*t^5 + 32*t^6 + O(t^7)
+        sage: p.nth_root(5)
+        4*t^2 + 8/5*t^3 + O(t^4)
+
+        sage: R.<x,y> = QQ[]
+        sage: S.<t> = R[[]]
+        sage: p = 1 + x*t + (x^2+y^2)*t^2 + O(t^3)
+        sage: p.nth_root(7)
+        1 + 1/7*x*t + (4/49*x^2 + 1/7*y^2)*t^2 + O(t^3)
+        """
+
+        if n != Integer(n):
+            raise ValueError, 'n should be an integer'
+        if self.is_zero():
+            if n > 0:
+                return self._parent(0).O(self.prec()/n)
+            else:
+                return 'ValueError', 'undefined result'
+        if self == 1:
+            return self
+        prec = self.prec()
+        if prec is infinity and self.degree() > 0:
+            prec = self._parent.default_prec()
+ 
+        if n == 0:
+            if self:
+                return 1
+            else:
+                raise ValueError, "0^0 expression"
+        if n == 1:
+            return self
+
+        if n < 0:
+            n = -n
+            sign = 1
+        else:
+            sign = 0
+        n = Integer(n)
+
+        val = self.valuation()
+        if val%n != 0:
+            raise NotImplementedError, "n does not divide the valuation"
+        gen = self.parent().gen()
+        if val:
+            self = self.valuation_zero_part()
+            prec -= val
+        try:
+            first_coeff = ~self[0]
+        except ValueError, ZeroDivisionError:
+            raise ZeroDivisionError, "leading coefficient must be a unit"
+
+        A = self.truncate()
+        R = A.parent()     # R is the corresponding polynomial ring
+        if first_coeff == 1:
+            current = R(1)
+        else:
+            #current = R(first_coeff**(1/n))
+            current = R(first_coeff.nth_root(n))
+        a_n = R(1+1/n)
+        b_n = R(1/n)
+        for next_prec in sage.misc.misc.newton_method_sizes(prec)[1:]:
+            z = current**(n+1) * A.truncate(next_prec)
+            current = a_n*current - b_n * z.truncate(next_prec)
+
+        res = self._parent(current, prec=prec)
+        if sign:
+            if val:
+                return res*gen**(-val/n)
+            return res
+        else:
+            p = res**-1
+            if val:
+                return p*gen**(val/n)
+            return p
+
     def sqrt(self, prec=None, extend=False, all=False, name=None):
         r""" The square root function.