Ticket #10720: trac_10720_power_series_nth_root_2.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: R.<t> = ZZ[[]]
+            sage: p = 1 + 4*t^2 + 3*t^3 + O(t^4)
+            sage: p.nth_root(2)                 
+            1 + 2*t^2 + 3/2*t^3 + O(t^4)
+            sage: R.<x,y> = QQ[]  
+            sage: S.<t> = R[[]]
+            sage: p = 8*t^3 + x*t^4 + (x^2+y^2)*t^5 + O(t^6)
+            sage: p.nth_root(3)     
+            2*t + 1/12*x*t^2 + (23/288*x^2 + 1/12*y^2)*t^3 + O(t^4)
+
+        TESTS::
+
+            sage: K.<x> = QQ[[]]
+            sage: (x^10/2).nth_root(2)
+            Traceback (most recent call last):
+            ...
+            ValueError: unable to take the 2-th root of 2 
+            sage: (x^2).nth_root(3)
+            Traceback (most recent call last):
+            ...
+            ValueError: power series does not have a 3-th root since 3 does not divide the valuation 2
+
+        AUTHORS:
+          
+        - Mario Pernici
+        """
+        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:
+                raise '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 ValueError, "power series does not have a %s-th root since %s does not divide the valuation %s" %(n,n,val)
+        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
+
+        # compute the initial value first_coeff^(1/n)
+        try:
+            c = first_coeff.nth_root(n)
+        except:
+            try:
+                c = first_coeff.constant_coefficient()
+            except:
+                raise ValueError, 'unable to take the %s-th root of %s' %(n,first_coeff)
+            if c == first_coeff:
+                c = c.nth_root(n)
+
+        if first_coeff.parent() != R.base_ring():
+            R = R.change_ring(first_coeff.parent())
+        current = R(c)
+
+        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)
+            current = a_n*current - b_n * z.truncate(next_prec)
+
+        S = self._parent
+        if current.base_ring() is not self.base_ring():
+            S1 = S.change_ring(current.base_ring())
+            res = S1(current, prec=prec)
+        else:
+            res = S(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.