Ticket #10720: trac_10720_power_series_nth_root.patch

File trac_10720_power_series_nth_root.patch, 3.8 KB (added by pernici, 12 years ago)
  • sage/rings/power_series_ring_element.pyx

    # HG changeset patch
    # User Mario Pernici <mario.pernici@gmail.com>
    # Date 1296497929 -3600
    # Node ID f6966619b689f57b56b7908134d90cb57baa2331
    # Parent  120c07be6358d93bcff503363d379c26b8342f2b
    nth_root in power series
    
    diff -r 120c07be6358 -r f6966619b689 sage/rings/power_series_ring_element.pyx
    a b  
    923923        #        b.append(-b[0]*sum([b[n-i]*a[i] for i in range(1,n+1) if i < len(a)]))
    924924        #return self.parent()(b, prec=prec)
    925925   
     926
    926927    def valuation_zero_part(self):
    927928        r"""
    928929        Factor self as as `q^n \cdot (a_0 + a_1 q + \cdots)` with
     
    10981099            except TypeError:
    10991100                return False
    11001101   
     1102    def nth_root(self, n):
     1103        """ Returns an `n^{th}` root of self
     1104
     1105        INPUT:
     1106       
     1107          - n - An integer number
     1108
     1109        ALGORITHM: Newton's method to compute $x = y^{\frac{-1}{n}}$
     1110       
     1111        .. math::
     1112       
     1113           x_{i+1} = (1 + \frac{1}{n}) x_i - \frac{1}{n} y x_{i}^{n+1})
     1114
     1115        EXAMPLES::
     1116       
     1117        sage: R.<t> = QQ[]
     1118        sage: p = 1 + t + t^2 + O(t^3)
     1119        sage: p.nth_root(2)
     1120        1 + 1/2*t + 3/8*t^2 + O(t^3)
     1121        sage: p.nth_root(-3)
     1122        1 - 1/3*t - 1/9*t^2 + O(t^3)
     1123        sage: p = 1024*t^10*(1 + 2*t + t^2) + O(t^12)
     1124        sage: p.nth_root(2)                         
     1125        32*t^5 + 32*t^6 + O(t^7)
     1126        sage: p.nth_root(5)
     1127        4*t^2 + 8/5*t^3 + O(t^4)
     1128
     1129        sage: R.<x,y> = QQ[]
     1130        sage: S.<t> = R[[]]
     1131        sage: p = 1 + x*t + (x^2+y^2)*t^2 + O(t^3)
     1132        sage: p.nth_root(7)
     1133        1 + 1/7*x*t + (4/49*x^2 + 1/7*y^2)*t^2 + O(t^3)
     1134        """
     1135
     1136        if n != Integer(n):
     1137            raise ValueError, 'n should be an integer'
     1138        if self.is_zero():
     1139            if n > 0:
     1140                return self._parent(0).O(self.prec()/n)
     1141            else:
     1142                return 'ValueError', 'undefined result'
     1143        if self == 1:
     1144            return self
     1145        prec = self.prec()
     1146        if prec is infinity and self.degree() > 0:
     1147            prec = self._parent.default_prec()
     1148 
     1149        if n == 0:
     1150            if self:
     1151                return 1
     1152            else:
     1153                raise ValueError, "0^0 expression"
     1154        if n == 1:
     1155            return self
     1156
     1157        if n < 0:
     1158            n = -n
     1159            sign = 1
     1160        else:
     1161            sign = 0
     1162        n = Integer(n)
     1163
     1164        val = self.valuation()
     1165        if val%n != 0:
     1166            raise NotImplementedError, "n does not divide the valuation"
     1167        gen = self.parent().gen()
     1168        if val:
     1169            self = self.valuation_zero_part()
     1170            prec -= val
     1171        try:
     1172            first_coeff = ~self[0]
     1173        except ValueError, ZeroDivisionError:
     1174            raise ZeroDivisionError, "leading coefficient must be a unit"
     1175
     1176        A = self.truncate()
     1177        R = A.parent()     # R is the corresponding polynomial ring
     1178        if first_coeff == 1:
     1179            current = R(1)
     1180        else:
     1181            #current = R(first_coeff**(1/n))
     1182            current = R(first_coeff.nth_root(n))
     1183        a_n = R(1+1/n)
     1184        b_n = R(1/n)
     1185        for next_prec in sage.misc.misc.newton_method_sizes(prec)[1:]:
     1186            z = current**(n+1) * A.truncate(next_prec)
     1187            current = a_n*current - b_n * z.truncate(next_prec)
     1188
     1189        res = self._parent(current, prec=prec)
     1190        if sign:
     1191            if val:
     1192                return res*gen**(-val/n)
     1193            return res
     1194        else:
     1195            p = res**-1
     1196            if val:
     1197                return p*gen**(val/n)
     1198            return p
     1199
    11011200    def sqrt(self, prec=None, extend=False, all=False, name=None):
    11021201        r""" The square root function.
    11031202