Ticket #10720: trac_10720_power_series_nth_root_2.patch

File trac_10720_power_series_nth_root_2.patch, 4.9 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 1296578728 -3600
    # Node ID 6232d130dc8af54e6a37c6fca6f92bbd4377d65f
    # Parent  120c07be6358d93bcff503363d379c26b8342f2b
    added nth_root for power series
    
    diff -r 120c07be6358 -r 6232d130dc8a 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: R.<t> = ZZ[[]]
     1124            sage: p = 1 + 4*t^2 + 3*t^3 + O(t^4)
     1125            sage: p.nth_root(2)                 
     1126            1 + 2*t^2 + 3/2*t^3 + O(t^4)
     1127            sage: R.<x,y> = QQ[] 
     1128            sage: S.<t> = R[[]]
     1129            sage: p = 8*t^3 + x*t^4 + (x^2+y^2)*t^5 + O(t^6)
     1130            sage: p.nth_root(3)     
     1131            2*t + 1/12*x*t^2 + (23/288*x^2 + 1/12*y^2)*t^3 + O(t^4)
     1132
     1133        TESTS::
     1134
     1135            sage: K.<x> = QQ[[]]
     1136            sage: (x^10/2).nth_root(2)
     1137            Traceback (most recent call last):
     1138            ...
     1139            ValueError: unable to take the 2-th root of 2
     1140            sage: (x^2).nth_root(3)
     1141            Traceback (most recent call last):
     1142            ...
     1143            ValueError: power series does not have a 3-th root since 3 does not divide the valuation 2
     1144
     1145        AUTHORS:
     1146         
     1147        - Mario Pernici
     1148        """
     1149        if n != Integer(n):
     1150            raise ValueError, 'n should be an integer'
     1151        if self.is_zero():
     1152            if n > 0:
     1153                return self._parent(0).O(self.prec()/n)
     1154            else:
     1155                raise 'ValueError', 'undefined result'
     1156        if self == 1:
     1157            return self
     1158        prec = self.prec()
     1159        if prec is infinity and self.degree() > 0:
     1160            prec = self._parent.default_prec()
     1161 
     1162        if n == 0:
     1163            if self:
     1164                return 1
     1165            else:
     1166                raise ValueError, "0^0 expression"
     1167        if n == 1:
     1168            return self
     1169
     1170        if n < 0:
     1171            n = -n
     1172            sign = 1
     1173        else:
     1174            sign = 0
     1175        n = Integer(n)
     1176
     1177        val = self.valuation()
     1178        if val%n != 0:
     1179            raise ValueError, "power series does not have a %s-th root since %s does not divide the valuation %s" %(n,n,val)
     1180        gen = self.parent().gen()
     1181        if val:
     1182            self = self.valuation_zero_part()
     1183            prec -= val
     1184        try:
     1185            first_coeff = ~self[0]
     1186        except ValueError, ZeroDivisionError:
     1187            raise ZeroDivisionError, "leading coefficient must be a unit"
     1188
     1189        A = self.truncate()
     1190        R = A.parent()     # R is the corresponding polynomial ring
     1191
     1192        # compute the initial value first_coeff^(1/n)
     1193        try:
     1194            c = first_coeff.nth_root(n)
     1195        except:
     1196            try:
     1197                c = first_coeff.constant_coefficient()
     1198            except:
     1199                raise ValueError, 'unable to take the %s-th root of %s' %(n,first_coeff)
     1200            if c == first_coeff:
     1201                c = c.nth_root(n)
     1202
     1203        if first_coeff.parent() != R.base_ring():
     1204            R = R.change_ring(first_coeff.parent())
     1205        current = R(c)
     1206
     1207        b_n = 1/n
     1208        a_n = 1+1/n
     1209        for next_prec in sage.misc.misc.newton_method_sizes(prec)[1:]:
     1210            z = current**(n+1) * A.truncate(next_prec)
     1211            current = a_n*current - b_n * z.truncate(next_prec)
     1212
     1213        S = self._parent
     1214        if current.base_ring() is not self.base_ring():
     1215            S1 = S.change_ring(current.base_ring())
     1216            res = S1(current, prec=prec)
     1217        else:
     1218            res = S(current, prec=prec)
     1219
     1220        if sign:
     1221            if val:
     1222                return res*gen**(-val/n)
     1223            return res
     1224        else:
     1225            p = res**-1
     1226            if val:
     1227                return p*gen**(val/n)
     1228            return p
     1229
    11011230    def sqrt(self, prec=None, extend=False, all=False, name=None):
    11021231        r""" The square root function.
    11031232