source: sage/rings/polynomial/padics/polynomial_padic_capped_relative_dense.py @ 8504:146366eccef3

Revision 8504:146366eccef3, 42.0 KB checked in by David Roe <roed@…>, 5 years ago (diff)

Fixing doctests and merge problems.
Line 
1"""
2p-adic Capped Relative Dense Polynomials
3"""
4
5import sage.rings.polynomial.polynomial_element_generic
6import sage.rings.polynomial.polynomial_element
7import sage.rings.polynomial.polynomial_integer_dense_ntl
8import sage.rings.integer
9import sage.rings.integer_ring
10import sage.rings.padics.misc as misc
11import sage.rings.padics.precision_error as precision_error
12import sage.rings.fraction_field_element as fraction_field_element
13import copy
14
15from sage.libs.all import pari, pari_gen
16from sage.libs.ntl.all import ZZX
17from sage.structure.factorization import Factorization
18from sage.rings.infinity import infinity
19
20min = misc.min
21ZZ = sage.rings.integer_ring.ZZ
22PrecisionError = precision_error.PrecisionError
23Integer = sage.rings.integer.Integer
24Polynomial = sage.rings.polynomial.polynomial_element.Polynomial
25is_Polynomial = sage.rings.polynomial.polynomial_element.is_Polynomial
26Polynomial_generic_domain = sage.rings.polynomial.polynomial_element_generic.Polynomial_generic_domain
27Polynomial_integer_dense = sage.rings.polynomial.polynomial_integer_dense_ntl.Polynomial_integer_dense_ntl
28
29class Polynomial_padic_capped_relative_dense(Polynomial_generic_domain):
30    def __init__(self, parent, x=None, check=True, is_gen=False, construct = False, absprec = infinity, relprec = infinity):
31        """
32        TESTS:
33            sage: K = Qp(13,7)
34            sage: R.<t> = K[]
35            sage: R([K(13), K(1)])
36            (1 + O(13^7))*t + (13 + O(13^8))
37            sage: T.<t> = ZZ[]
38            sage: R(t + 2)
39            (1 + O(13^7))*t + (2 + O(13^7))
40        """
41        Polynomial.__init__(self, parent, is_gen=is_gen)
42        if construct:
43            (self._poly, self._valbase, self._relprecs, self._normalized, self._valaddeds, self._list) = x #the last two of these may be None
44            return
45        elif is_gen:
46            from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
47            self._poly = PolynomialRing(ZZ, parent.variable_name()).gen()
48            self._valbase = 0
49            self._valaddeds = [infinity, 0]
50            self._relprecs = [infinity, parent.base_ring().precision_cap()]
51            self._normalized = True
52            self._list = None
53            return
54
55        #First we list the types that are turned into Polynomials
56        if isinstance(x, ZZX):
57            from polynomial_ring_constructor import PolynomialRing
58            x = Polynomial_integer_dense(PolynomialRing(ZZ, parent.variable_name()), x, construct = True)
59        elif isinstance(x, fraction_field_element.FractionFieldElement) and \
60               x.denominator() == 1:
61            #Currently we ignore precision information in the denominator.  This should be changed eventually
62            x = x.numerator()
63
64        #We now coerce various types into lists of coefficients.  There are fast pathways for some types of polynomials
65        if isinstance(x, Polynomial):
66            if x.parent() is self.parent():
67                if not absprec is infinity or not relprec is infinity:
68                    x._normalize()
69                self._poly = x._poly
70                self._valbase = x._valbase
71                self._valaddeds = x._valaddeds
72                self._relprecs = x._relprecs
73                self._normalized = x._normalized
74                self._list = x._list
75                if not absprec is infinity or not relprec is infinity:
76                    self._adjust_prec_info(absprec, relprec)
77                return
78            elif x.base_ring() is ZZ:
79                self._poly = x
80                self._valbase = Integer(0)
81                p = parent.base_ring().prime()
82                self._relprecs = [c.valuation(p) + parent.base_ring().precision_cap() for c in x.list()]
83                self._comp_valaddeds()
84                self._normalized = len(self._valaddeds) == 0 or (min(self._valaddeds) == 0)
85                self._list = None
86                if not absprec is infinity or not relprec is infinity:
87                    self._adjust_prec_info(absprec, relprec)
88                return
89            else:
90                x = [parent.base_ring()(a) for a in x.list()]
91                check = False
92        elif isinstance(x, dict):
93            zero = parent.base_ring()(0)
94            n = max(x.keys())
95            v = [zero for _ in xrange(n + 1)]
96            for i, z in x.iteritems():
97                v[i] = z
98            x = v
99        elif isinstance(x, pari_gen):
100            x = [parent.base_ring()(w) for w in x.Vecrev()]
101            check = False
102        #The default behavior if we haven't already figured out what the type is is to assume it coerces into the base_ring as a constant polynomial
103        elif not isinstance(x, list):
104            x = [x] # constant polynomial
105       
106        if check:
107            x = [parent.base_ring()(z) for z in x]
108
109        from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
110
111        # Remove this -- for p-adics this is terrible, since it kills any non exact zero.
112        #if len(x) == 1 and not x[0]:
113        #    x = []
114       
115        self._list = x
116        self._valaddeds = [a.valuation() for a in x]
117        self._valbase = sage.rings.padics.misc.min(self._valaddeds)
118        if self._valbase is infinity:
119            self._valaddeds = []
120            self._relprecs = []
121            self._poly = PolynomialRing(ZZ, parent.variable_name())()
122            self._normalized = True
123            if not absprec is infinity or not relprec is infinity:
124                self._adjust_prec_info(absprec, relprec)
125        else:
126            self._valaddeds = [c - self._valbase for c in self._valaddeds]
127            self._relprecs = [a.precision_absolute() - self._valbase for a in x]
128            self._poly = PolynomialRing(ZZ, parent.variable_name())([a >> self._valbase for a in x])
129            self._normalized = True
130            if not absprec is infinity or not relprec is infinity:
131                self._adjust_prec_info(absprec, relprec)
132
133    def _normalize(self):
134        # Currently slow: need to optimize
135        if not self._normalized:
136            if self._valaddeds is None:
137                self._comp_valaddeds()
138            val = sage.rings.padics.misc.min(self._valaddeds)
139            prime_pow = self.base_ring().prime_pow
140            selflist = self._poly.list()
141            if val is infinity:
142                pass
143            elif val != 0:
144                self._relprecs = [max(prec - val,0) for prec in self._relprecs]
145                v = [Integer(0) if (e is infinity) else ((c // prime_pow(val)) % prime_pow(e)) for (c,e) in zip(selflist, self._relprecs)]
146                self._poly = self._poly.parent()(v, check=False)
147                self._valbase += val
148                self._valaddeds = [c - val for c in self._valaddeds]
149            else:
150                self._poly = self._poly.parent()([Integer(0) if (e is infinity) else (c % prime_pow(e)) for (c,e) in zip(selflist, self._relprecs)], check=False)
151            self._normalized = True
152
153    def _reduce_poly(self):
154        selflist = self._poly.list()
155        prime_pow = self.base_ring().prime_pow
156        self._poly = self._poly.parent()([Integer(0) if (e is infinity) else (c % prime_pow(e)) for (c, e) in zip(selflist, self._relprecs)], check=False)
157
158    def __reduce__(self):
159        """
160        For pickling.  This function is here because the relative precisions were getting screwed up for some reason.
161        """
162        return make_padic_poly, (self.parent(), (self._poly, self._valbase, self._relprecs, self._normalized, self._valaddeds, self._list), 0)
163
164    def _comp_list(self):
165        if self.degree() == -1 and self._valbase == infinity:
166            self._list = []
167            return self._list
168        polylist = self._poly.list()
169        polylen = len(polylist)
170        self._list = [self.base_ring()(polylist[i], absprec = self._relprecs[i]) << self._valbase for i in range(polylen)] \
171                     + [self.base_ring()(0, absprec = self._relprecs[i] + self._valbase) for i in range(polylen, len(self._relprecs))]
172        while self._list[-1]._is_exact_zero():
173            self._list.pop()
174
175    def _comp_valaddeds(self):
176        self._valaddeds = []
177        for i in range(self._poly.degree() + 1):
178            tmp = self._poly.list()[i].valuation(self.parent().base_ring().prime())
179            if tmp is infinity or tmp > self._relprecs[i]:
180                self._valaddeds.append(self._relprecs[i])
181            else:
182                self._valaddeds.append(tmp)
183        for i in range(self._poly.degree() + 1, len(self._relprecs)):
184            self._valaddeds.append(self._relprecs[i])
185
186    def _adjust_prec_info(self, absprec=infinity, relprec=infinity):
187        r"""
188        Assumes that self._poly, self._val and self._relprec are set initially and adjusts self._val and self._relprec to the termwise minimum of absprec and relprec.
189        """
190        return
191
192#         min = sage.rings.padics.misc.min
193#         slen = len(self._relprec)
194#         if isinstance(absprec, list):
195#             alen = len(absprec)
196#         elif absprec is infinity:
197#             alen = 0
198#             absprec = []
199#         else:
200#             alen = 1
201#         if isinstance(relprec, list):
202#             rlen = len(relprec)
203#         elif relprec is infinity:
204#             rlen = 0
205#             relprec = []
206#         else:
207#             rlen = 1
208#         preclen = max(slen, rlen, alen)
209#         if not isinstance(absprec, list):
210#             absprec = [absprec] * preclen
211#         if not isinstance(relprec, list):
212#             relprec = [relprec] * preclen
213#         vallist = [c.valuation(self.base_ring().prime()) + self._val for c in self._poly.list()] #######
214#         vmin = min(vallist)
215#         amin = min(absprec)
216#         if amin < vmin:
217#             vmin = amin
218#         if vmin < self._val:
219#             vadjust =
220
221#         if not isinstance(absprec, list):
222#             self._val = min(vallist + [absprec])
223#             absprec = [absprec] * preclen
224#         else:
225#             self._val = padics.misc.min(vallist + absprec)
226#             absprec = absprec + [infinity] * (preclen - len(absprec))
227#         if self._val is infinity:
228#             self._relprec = []
229#             return
230#         if not isinstance(relprec, list):
231#             relprec = [relprec] * preclen
232#         else:
233#             relprec = relprec + [parent.base_ring().precision_cap()] * (preclen - len(relprec))
234#         self._relprec = [min(a, v + r) - self._val for (a, r, v) in zip(absprec, relprec, vallist)]
235#Remember to normalize at the end if self._normalized is true because you need to reduce mod p^n
236
237    def _getprecpoly(self, n):
238        one = Integer(1)
239        return self._poly.parent()([(0 if (c is infinity) else (one << (n * c))) for c in self._relprecs])
240
241    def _getvalpoly(self, n):
242        one = Integer(1)
243        if self._valaddeds is None:
244            self._comp_valaddeds()
245        return self._poly.parent()([(0 if (c is infinity) else (one << (n * c))) for c in self._valaddeds] + \
246                                   [(0 if (c is infinity) else (one << (n * c))) for c in self._relprecs[len(self._valaddeds):]])
247
248    def list(self):
249        """
250        Returns a list of coefficients of self.
251
252        NOTE:
253        The length of the list returned may be greater
254        than expected since it includes any leading zeros
255        that have finite absolute precision.
256
257        EXAMPLES:
258        sage: K = Qp(13,7)
259        sage: R.<t> = K[]       
260        sage: a = 2*t^3 + 169*t - 1
261        sage: a
262        (2 + O(13^7))*t^3 + (13^2 + O(13^9))*t + (12 + 12*13 + 12*13^2 + 12*13^3 + 12*13^4 + 12*13^5 + 12*13^6 + O(13^7))
263        sage: a.list()
264        [12 + 12*13 + 12*13^2 + 12*13^3 + 12*13^4 + 12*13^5 + 12*13^6 + O(13^7),
265         13^2 + O(13^9),
266         0,
267         2 + O(13^7)]
268         """
269
270        if self._list is None:
271            self._comp_list()
272        return list(self._list)
273
274    def _repr(self, name=None):
275        """
276        TESTS:
277            sage: k = Qp(5,10)
278            sage: R.<x> = k[]
279            sage: f = R([k(0,-3), 0, k(0,-1)]); f
280            (O(5^-1))*x^2 + (O(5^-3))
281            sage: f + f
282            (O(5^-1))*x^2 + (O(5^-3))       
283        """
284        # TODO: what is new here (that doesn't come from parent class)?
285        s = " "
286        coeffs = self.list()
287        m = len(coeffs)
288        while m > 0 and coeffs[m-1].valuation() == infinity:
289            m -= 1
290        r = reversed(xrange(m))
291        if name is None:
292            name = self.parent().variable_name()
293        for n in r:
294            x = coeffs[n]
295            if x.valuation() < infinity:
296                if n != m-1:
297                    s += " + "
298                x = "(%s)"%x
299                if n > 1:
300                    var = "*%s^%s"%(name,n)
301                elif n==1:
302                    var = "*%s"%name
303                else:
304                    var = ""
305                s += "%s%s"%(x,var)
306        if s==" ":
307            return "0"
308        return s[1:]
309
310    def content(self):
311        """
312        Returns the content of self.
313
314        EXAMPLES:
315        sage: K = Qp(13,7)
316        sage: R.<t> = K[]       
317        sage: a = 13^7*t^3 + K(169,4)*t - 13^4
318        sage: a.content()
319        13^2 + O(13^9)
320        """
321        if self._normalized:
322            return self._valbase
323        if self._valaddeds is None:
324            self._comp_valaddeds()
325        return self.base_ring()(self.base_ring().prime_pow(min(self._valaddeds) + self._valbase))
326
327    def lift(self):
328        """
329        Returns an integer polynomial congruent to this one modulo the precision of each coefficient.
330
331        NOTE: The lift that is returned will not necessarily be the same for polynomials with
332              the same coefficients (ie same values and precisions): it will depend on how
333              the polynomials are created.
334
335        EXAMPLES:
336        sage: K = Qp(13,7)
337        sage: R.<t> = K[]       
338        sage: a = 13^7*t^3 + K(169,4)*t - 13^4
339        sage: a.lift()
340        62748517*t^3 + 169*t - 28561
341        """
342        return self.base_ring().prime_pow(self._valbase) * self._poly
343
344    def __getitem__(self, n):
345        """
346        Returns the coefficient of x^n
347
348        EXAMPLES:
349        sage: K = Qp(13,7)
350        sage: R.<t> = K[]       
351        sage: a = 13^7*t^3 + K(169,4)*t - 13^4
352        sage: a[1]
353        13^2 + O(13^4)
354        """
355        if n >= len(self._relprecs):
356            return self.base_ring()(0)
357        if not self._list is None:
358            return self._list[n]
359        return self.base_ring()(self.base_ring().prime_pow(self._valbase) * self._poly[n], absprec = self._valbase + self._relprecs[n])
360
361    def __getslice__(self, i, j):
362        """
363        EXAMPLES:
364        sage: K = Qp(13,7)
365        sage: R.<t> = K[]       
366        sage: a = 13^7*t^3 + K(169,4)*t - 13^4
367        sage: a[1:2]
368        (13^2 + O(13^4))*t
369        """
370        if i < 0:
371            i = len(self._relprecs) + i
372            if i < 0:
373                raise IndexError, "list index out of range"
374        if j > len(self._relprecs):
375            j = len(self._relprecs)
376        elif j < 0:
377            j = len(self._relprecs) + j
378            if j < 0:
379                raise IndexError, "list index out of range"
380        if i >= j:
381            return Polynomial_padic_capped_relative_dense(self.parent(), [])
382        else:
383            return Polynomial_padic_capped_relative_dense(self.parent(), (self._poly[i:j], self._valbase, [infinity]*i + self._relprecs[i:j], False, None if self._valaddeds is None else [infinity]*i + self._valaddeds[i:j], None if self._list is None else [self.base_ring()(0)] * i + self._list[i:j]), construct = True)
384
385    def _add_(self, right):
386        """
387        Returns the sum of self and right.
388
389        EXAMPLES:
390        sage: K = Qp(13,7)
391        sage: R.<t> = K[]       
392        sage: a = t^4 + 17*t^2 + 1
393        sage: b = -t^4 + 9*t^2 + 13*t - 1
394        sage: c = a + b; c
395        (O(13^7))*t^4 + (2*13 + O(13^7))*t^2 + (13 + O(13^8))*t + (O(13^7))
396        sage: c.list()
397        [O(13^7), 13 + O(13^8), 2*13 + O(13^7), 0, O(13^7)]
398        """
399        selfpoly = self._poly
400        rightpoly = right._poly
401        if self._valbase > right._valbase:
402            selfpoly = selfpoly * self.base_ring().prime_pow(self._valbase - right._valbase)
403            baseval = right._valbase
404        elif self._valbase < right._valbase:
405            rightpoly = rightpoly * self.base_ring().prime_pow(right._valbase - self._valbase)
406            baseval = self._valbase
407        else:
408            baseval = self._valbase
409        # Currently we don't reduce the coefficients of the answer modulo the appropriate power of p or normalize
410        return Polynomial_padic_capped_relative_dense(self.parent(), \
411                                                      (selfpoly + rightpoly, \
412                                                       baseval, \
413                                                       [min(a + self._valbase - baseval, b + right._valbase - baseval) for (a, b) in
414                                                              zip(_extend_by_infinity(self._relprecs, max(len(self._relprecs), len(right._relprecs))), \
415                                                                  _extend_by_infinity(right._relprecs, max(len(self._relprecs), len(right._relprecs))))], \
416                                                       False, None, None), construct = True)
417
418    def _sub_(self, right):
419        """
420        Returns the sum of self and right.
421
422        EXAMPLES:
423        sage: K = Qp(13,7)
424        sage: R.<t> = K[]       
425        sage: a = t^4 + 17*t^2 + 1
426        sage: b = t^4 - 9*t^2 - 13*t + 1
427        sage: c = a - b; c
428        (O(13^7))*t^4 + (2*13 + O(13^7))*t^2 + (13 + O(13^8))*t + (O(13^7))
429        sage: c.list()
430        [O(13^7), 13 + O(13^8), 2*13 + O(13^7), 0, O(13^7)]
431        """
432        selfpoly = self._poly
433        rightpoly = right._poly
434        if self._valbase > right._valbase:
435            selfpoly = selfpoly * self.base_ring().prime_pow(self._valbase - right._valbase)
436            baseval = right._valbase
437        elif self._valbase < right._valbase:
438            rightpoly = rightpoly * self.base_ring().prime_pow(right._valbase - self._valbase)
439            baseval = self._valbase
440        else:
441            baseval = self._valbase
442        # Currently we don't reduce the coefficients of the answer modulo the appropriate power of p or normalize
443        return Polynomial_padic_capped_relative_dense(self.parent(), \
444                                                      (selfpoly - rightpoly, \
445                                                       baseval, \
446                                                       [min(a + self._valbase - baseval, b + right._valbase - baseval) for (a, b) in
447                                                              zip(_extend_by_infinity(self._relprecs, max(len(self._relprecs), len(right._relprecs))), \
448                                                                  _extend_by_infinity(right._relprecs, max(len(self._relprecs), len(right._relprecs))))], \
449                                                       False, None, None), construct = True)
450
451    def _mul_(self, right):
452        r"""
453        Multiplies self and right.
454
455        ALGORITHM: We use an algorithm thought up by Joe Wetherell to
456        find the precisions of the product.  It works as follows:
457        Suppose $f = \sum_i a_i x^i$ and $g = \sum_j b_j x^j$. Let $N
458        = \max(\deg f, \deg g) + 1$ (in the actual implementation we
459        use $N = 2^{\lfloor \log_2\max(\deg f, \deg g)\rfloor + 1}$).
460        The valuations and absolute precisions of each coefficient
461        contribute to the absolute precision of the kth coefficient of
462        the product in the following way: for each $i + j = k$, you
463        take the valuation of $a_i$ plus the absolute precision of
464        $b_j$, and then take the valuation of $b_j$ plus the absolute
465        precision of $a_i$, take the minimum of those two, and then
466        take the minimum over all $i$, $j$ summing to $k$.
467
468        You can simulate this as follows. Construct new polynomials of
469        degree $N$:
470
471        \begin{align*}
472        A &= \sum_i N^{\mbox{valuation of $a_i$}} x^i \\
473        B &= \sum_j N^{\mbox{absolute precision of $b_j$}} x^j \\
474        C &= \sum_i N^{\mbox{absolute precision of $a_i$}} x^i \\
475        D &= \sum_j N^{\mbox{valuation of $b_j$}} x^j \\
476        \end{align*}
477
478        Now you compute AB and CD. Because you're representing things
479        'N-adically', you don't get any 'overflow', and you can just
480        read off what the precisions of the product are. In fact it
481        tells you more, it tells you exactly how many terms of each
482        combination of valuation modulus contribute to each term of
483        the product (though this feature is not currently exposed in
484        our implementation.
485
486        Since we're working 'N-adically' we can just consider
487        $N^{\infty} = 0$.
488
489        NOTE: The timing of normalization in arithmetic operations
490        may very well change as we do more tests on the relative time
491        requirements of these operations.
492
493        EXAMPLES:
494        sage: K = Qp(13,7)
495        sage: R.<t> = K[]       
496        sage: a = t^4 + 17*t^2 + 1
497        sage: b = -t^4 + 9*t^2 + 13*t - 1
498        sage: c = a + b; c
499        (O(13^7))*t^4 + (2*13 + O(13^7))*t^2 + (13 + O(13^8))*t + (O(13^7))
500        sage: d = R([K(1,4), K(2, 6), K(1, 5)]); d
501        (1 + O(13^5))*t^2 + (2 + O(13^6))*t + (1 + O(13^4))
502        sage: e = c * d; e
503        (O(13^7))*t^6 + (O(13^7))*t^5 + (2*13 + O(13^6))*t^4 + (5*13 + O(13^6))*t^3 + (4*13 + O(13^5))*t^2 + (13 + O(13^5))*t + (O(13^7))
504        sage: e.list()
505        [O(13^7),
506         13 + O(13^5),
507         4*13 + O(13^5),
508         5*13 + O(13^6),
509         2*13 + O(13^6),
510         O(13^7),
511         O(13^7)]
512        """
513        self._normalize()
514        right._normalize()
515        zzpoly = self._poly * right._poly
516        if len(self._relprecs) == 0 or len(right._relprecs) == 0:
517            return self.parent()(0)
518        n = Integer(len(self._relprecs) + len(right._relprecs) - 1).exact_log(2) + 1
519        precpoly1 = self._getprecpoly(n) * right._getvalpoly(n)
520        precpoly2 = self._getvalpoly(n) * right._getprecpoly(n)
521        # These two will be the same length
522        tn = Integer(1) << n
523        preclist = [min(a.valuation(tn), b.valuation(tn)) for (a, b) in zip(precpoly1.list(), precpoly2.list())]
524        answer = Polynomial_padic_capped_relative_dense(self.parent(), (zzpoly, self._valbase + right._valbase, preclist, False, None, None), construct = True)
525        answer._reduce_poly()
526        return answer
527
528    def _lmul_(self, right):
529        return self._rmul_(right)
530
531    def _rmul_(self, left):
532        """
533        Returns self multiplied by a constant
534
535        EXAMPLES:
536        sage: K = Qp(13,7)
537        sage: R.<t> = K[]
538        sage: a = t^4 + K(13,5)*t^2 + 13
539        sage: K(13,7) * a
540        (13 + O(13^7))*t^4 + (13^2 + O(13^6))*t^2 + (13^2 + O(13^8))
541        """
542        if self._valaddeds is None:
543            self._comp_valaddeds()
544        if left != 0:
545            val, unit = left._val_unit()
546            left_rprec = left.precision_relative()
547            relprecs = [min(left_rprec + self._valaddeds[i], self._relprecs[i]) for i in range(len(self._relprecs))]
548        elif left._is_exact_zero():
549            return Polynomial_padic_capped_relative_dense(self.parent(), [])
550        else:
551            return Polynomial_padic_capped_relative_dense(self.parent(), (self._poly.parent()(0), self._valbase + left.valuation(), self._valaddeds, False, self._valaddeds, None), construct = True)
552        return Polynomial_padic_capped_relative_dense(self.parent(), (self._poly._rmul_(unit), self._valbase + val, relprecs, False, self._valaddeds, None), construct = True)
553
554    def _neg_(self):
555        """
556        Returns the negation of self.
557
558        EXAMPLES:
559        sage: K = Qp(13,2)
560        sage: R.<t> = K[]
561        sage: a = t^4 + 13*t^2 + 4
562        sage: -a
563        (12 + 12*13 + O(13^2))*t^4 + (12*13 + 12*13^2 + O(13^3))*t^2 + (9 + 12*13 + O(13^2))
564        """
565        return Polynomial_padic_capped_relative_dense(self.parent(), (-self._poly, self._valbase, self._relprecs, False, self._valaddeds, None), construct = True)
566
567    def lshift_coeffs(self, shift, no_list = False):
568        """
569        Returns a new polynomials whose coefficients are multiplied by p^shift.
570
571        EXAMPLES:
572        sage: K = Qp(13, 4)
573        sage: R.<t> = K[]
574        sage: a = t + 52
575        sage: a.lshift_coeffs(3)
576        (13^3 + O(13^7))*t + (4*13^4 + O(13^8))
577        """
578        if shift < 0:
579            return self.rshift_coeffs(-shift, no_list)
580        if no_list or self._list is None:
581            return Polynomial_padic_capped_relative_dense(self.parent(), (self._poly, self._valbase + shift, self._relprecs, False, self._valaddeds, None), construct = True)
582        else:
583            return Polynomial_padic_capped_relative_dense(self.parent(), (self._poly, self._valbase + shift, self._relprecs, False, self._valaddeds, [c.__lshift__(shift) for c in self._list]), construct = True)
584
585    def rshift_coeffs(self, shift, no_list = False):
586        """
587        Returns a new polynomial whose coefficients are p-adiclly
588        shifted to the right by shift.
589
590        NOTES: Type Qp(5)(0).__rshift__? for more information.
591
592        EXAMPLES:
593        sage: K = Zp(13, 4)
594        sage: R.<t> = K[]
595        sage: a = t^2 + K(13,3)*t + 169; a
596        (1 + O(13^4))*t^2 + (13 + O(13^3))*t + (13^2 + O(13^6))
597        sage: b = a.rshift_coeffs(1); b
598        (O(13^3))*t^2 + (1 + O(13^2))*t + (13 + O(13^5))
599        sage: b.list()
600        [13 + O(13^5), 1 + O(13^2), O(13^3)]
601        sage: b = a.rshift_coeffs(2); b
602        (O(13^2))*t^2 + (O(13))*t + (1 + O(13^4))
603        sage: b.list()
604        [1 + O(13^4), O(13), O(13^2)]
605        """
606        if shift < 0:
607            return self.lshift_coeffs(-shift, no_list) # We can't just absorb this into the next if statement because we allow rshift to preserve _normalized
608        if self.base_ring().is_field() or shift <= self._valbase:
609            if no_list or self._list is None:
610                return Polynomial_padic_capped_relative_dense(self.parent(), (self._poly, self._valbase - shift, self._relprecs, self._normalized, self._valaddeds, None), construct = True)
611            else:
612                return Polynomial_padic_capped_relative_dense(self.parent(), (self._poly, self._valbase - shift, self._relprecs, self._normalized, self._valaddeds, [c.__rshift__(shift) for c in self._list]), construct = True)
613        else:
614            shift = shift - self._valbase
615            fdiv = self.base_ring().prime_pow(shift)
616            return Polynomial_padic_capped_relative_dense(self.parent(), (self._poly // fdiv, 0, [0 if a <= shift else a - shift for a in self._relprecs], False, None, None), construct = True)
617
618    #def __floordiv__(self, right):
619    #    if is_Polynomial(right) and right.is_constant() and right[0] in self.base_ring():
620    #        d = self.base_ring()(right[0])
621    #    elif (right in self.base_ring()):
622    #        d = self.base_ring()(right)
623    #    else:
624    #        raise NotImplementedError
625    #    return self._rmul_(self.base_ring()(~d.unit_part())).rshift_coeffs(d.valuation())
626
627    def _unsafe_mutate(self, n, value):
628        """
629        It's a really bad idea to use this function for p-adic polynomials.  There are speed issues, and it may not be bug-free currently.
630        """
631        n = int(n)
632        value = self.base_ring()(value)
633        if self.is_gen():
634            raise ValueError, "cannot modify generator"
635        if n < 0:
636            raise IndexError, "n must be >= 0"
637        if self._valbase is infinity:
638            if value._is_exact_zero():
639                return
640            self._valbase = value.valuation()
641            if value != 0:
642                self._poly._unsafe_mutate(self, n, value.unit_part().lift())
643                self._relprecs = [infinity] * n + [value.precision_relative()]
644            else:
645                self._relprecs = [infinity] * n + [0]
646            self._valaddeds = [infinity] * n + [0]
647            zero = self.base_ring()(0)
648            self._list = [zero] * n + [value]
649            self._normalized = True
650        elif value.valuation() >= self._valbase:
651            # _valbase and _normalized stay the same
652            if value != 0:
653                self._poly._unsafe_mutate(self, n, (value.__rshift__(self._valbase)).lift())
654            else:
655                self._poly._unsafe_mutate(self, n, 0)
656            if n < len(self._relprecs):
657                self._relprecs[n] = value.precision_absolute() - self._valbase
658                if not self._valaddeds is None:
659                    self._valaddeds[n] = value.valuation() - self._valbase
660                if not self._list is None:
661                    self._list[n] = value
662            else:
663                self._relprecs.extend([infinity] * (n - len(self._relprecs)) + [value.precision_absolute() - self._valbase])
664                if not self._valaddeds is None:
665                    self._valaddeds.extend([infinity] * (n - len(self._relprecs)) + [value.valuation() - self._valbase])
666                if not self._list is None:
667                    zero = self.base_ring()(0)
668                    self._list.extend([zero] * (n - len(self._relprecs)) + [value])
669        else:
670            basediff = self._valbase - value.valuation()
671            self._valbase = value.valuation()
672            if not self._valaddeds is None:
673                self._valaddeds = [c + basediff for c in self._valaddeds]
674            self._poly = self._poly * self.base_ring().prime_pow(basediff)
675            if value != 0:
676                self._poly._unsafe_mutate(self, n, value.unit_part().lift())
677            else:
678                self._poly._unsafe_mutate(self, n, 0)
679            if n < len(self._relprecs):
680                self._relprecs[n] = value.precision_relative()
681            else:
682                self._relprecs.extend([infinity] * (n - len(self._relprecs)) + [value.precision_relative()])
683            self._normalized = False
684            if not self._list is None:
685                if n < len(self._list):
686                    self._list[n] = value
687                else:
688                    zero = self._base_ring()(0)
689                    self._list.extend([zero] * (n - len(self._list)) + [value])
690
691    def _pari_(self, variable = None):
692        if variable is None:
693            variable = self.parent().variable_name()
694        return pari(self.list()).Polrev(variable)
695
696    def copy(self):
697        return Polynomial_padic_capped_relative_dense(self.parent(), (self._poly.copy(), self._valbase, copy.copy(self._relprecs), self._normalized, copy.copy(self._valaddeds), copy.copy(self._list)), construct = True)
698
699    def degree(self):
700        """
701        Returns the degree of self, i.e., the largest $n$ so that the
702        coefficient of $x^n$ does not compare equal to $0$.
703
704        EXAMPLES:
705            sage: K = Qp(3,10)
706            sage: x = O(3^5)
707            sage: li =[3^i * x for i in range(0,5)]; li
708            [O(3^5), O(3^6), O(3^7), O(3^8), O(3^9)]
709            sage: R.<T> = K[]
710            sage: f = R(li); f
711            (O(3^9))*T^4 + (O(3^8))*T^3 + (O(3^7))*T^2 + (O(3^6))*T + (O(3^5))
712            sage: f.degree()
713            -1
714        """
715        self._normalize()
716        return self._poly.degree()
717
718    def prec_degree(self):
719        """
720        Returns the largest $n$ so that precision information is
721        stored about the coefficient of $x^n$.
722
723        Always greater than or equal to degree.
724        """
725        return len(self._relprecs) - 1
726
727    def precision_absolute(self, n = None):
728        """
729        Returns absolute precision information about self.
730
731        INPUT:
732        self -- a p-adic polynomial
733        n -- None or an integer (default None). 
734
735        OUTPUT:
736        If n == None, returns a list of absolute precisions of coefficients.  Otherwise,
737        returns the absolute precision of the coefficient of x^n.
738        """
739        if n is None:
740            return [c + self._valbase for c in self._relprecs]
741        return self._relprecs[n] + self._valbase
742
743    def precision_relative(self, n = None):
744        """
745        Returns relative precision information about self.
746
747        INPUT:
748        self -- a p-adic polynomial
749        n -- None or an integer (default None). 
750
751        OUTPUT:
752        If n == None, returns a list of relative precisions of coefficients.  Otherwise,
753        returns the relative precision of the coefficient of x^n.
754        """
755        if n is None:
756            self._normalize()
757            return copy.copy(self._relprecs)
758        n = int(n)
759        if n < 0 or n >= len(self._relprecs) or self._relprecs[n] is infinity:
760            return Integer(0)
761        if self._valaddeds is None:   
762            return self._relprecs[n] - self._poly[n].valuation(self.base_ring().prime())
763        else:
764            return self._relprecs[n] - self._valaddeds[n]
765
766    def valuation_of_coefficient(self, n = None):
767        """
768        Returns valuation information about self's coefficients.
769
770        INPUT:
771        self -- a p-adic polynomial
772        n -- None or an integer (default None). 
773
774        OUTPUT:
775        If n == None, returns a list of valuations of coefficients.  Otherwise,
776        returns the valuation of the coefficient of x^n.
777        """
778        if n is None:
779            self._normalize()
780            return [c + self._valbase for c in self._valadded]
781        n = int(n)
782        if n < 0 or n >= len(self._relprecs):
783            return infinity
784        if self._valaddeds is None:
785            return self._valbase + self._poly[n].valuation(self.base_ring().prime())
786        else:
787            return self._valbase + self._valaddeds[n]
788
789    def valuation(self, val_of_var = None):
790        """
791        Returns the valuation of self
792
793        INPUT:
794        self -- a p-adic polynomial
795        val_of_var -- None or a rational (default None). 
796
797        OUTPUT:
798        If val_of_var == None, returns the largest power of the variable dividing self.  Otherwise,
799        returns the valuation of self where the variable is assigned valuation val_of_var
800        """
801        if val_of_var is None:
802            return self._poly.valuation()
803        if self._valaddeds is None:
804            self._comp_valaddeds()
805        return self._valbase + min([self._valaddeds[i] + val_of_var * i for i in range(len(self._valaddeds))])
806
807    def reverse(self, n = None):
808        """
809        Returns a new polynomial whose coefficients are the reversed coefficients of self, where self is considered as a polynomial of degree n.
810
811        If n is None, defaults to the degree of self.
812        If n is smaller than the degree of self, some coefficients will be discarded.
813
814        EXAMPLES:
815        sage: K = Qp(13,7)
816        sage: R.<t> = K[]
817        sage: f = t^3 + 4*t; f
818        (1 + O(13^7))*t^3 + (4 + O(13^7))*t
819        sage: f.reverse()
820        (4 + O(13^7))*t^2 + (1 + O(13^7))
821        sage: f.reverse(3)
822        (4 + O(13^7))*t^2 + (1 + O(13^7))
823        sage: f.reverse(2)
824        (4 + O(13^7))*t
825        sage: f.reverse(4)
826        (4 + O(13^7))*t^3 + (1 + O(13^7))*t
827        sage: f.reverse(6)
828        (4 + O(13^7))*t^5 + (1 + O(13^7))*t^3
829        """
830        if n is None:
831            n = self._poly.degree()
832        zzlist = self._poly.list()[:(n+1)] + [0] * (n - self._poly.degree())
833        zzlist.reverse()
834        relprec = self._relprecs[:(n+1)] + [infinity] * (n - self.prec_degree())
835        relprec.reverse()
836        if self._valaddeds is None:
837            valadded = None
838        else:
839            valadded = self._valaddeds[:(n+1)] + [infinity] * (n - self.prec_degree())
840            valadded.reverse()
841        if self._list is None:
842            L = None
843        else:
844            L = self._list[:(n+1)] + [self.base_ring()(0)] * (n - self.prec_degree())
845            L.reverse()
846        return Polynomial_padic_capped_relative_dense(self.parent(), (self._poly.parent()(zzlist), self._valbase, relprec, self._normalized, valadded, L), construct = True)
847
848    def rescale(self, a):
849        r"""
850        Return f(a*X)
851
852        NOTE:  Need to write this function for integer polynomials before this works.
853
854        EXAMPLES:
855            sage: K = Zp(13, 5)
856            sage: R.<t> = K[]
857            sage: f = t^3 + K(13, 3) * t
858            sage: f.rescale(2)    # todo: not tested -- in fact, is broken!
859        """
860        negval = False
861        try:
862            a = self.base_ring()(a)
863        except ValueError, msg:
864            if msg == "element has negative valuation.":
865                negval = True
866            else:
867                raise ValueError, msg
868        if negval:
869            return self.parent().base_extend(self.base_ring().fraction_field())(self).rescale(a)
870        if self.base_ring().is_field() and a.valuation() < 0:
871            D = self.prec_degree()
872            return a**D * self.reverse(D).rescale(~a).reverse(D)
873        aval = a.valuation()
874        arprec = a.precision_relative()
875        if self._valaddeds is None:
876            self._comp_valaddeds()
877        valadded = [self._valaddeds[i] + aval * i for i in range(len(self._valaddeds))]
878        relprec = [infinity if (self._relprecs[i] is infinity) else (min(self._relprecs[i] - self._valaddeds[i], arprec) + aval * i + self._valaddeds[i]) for i in range(len(self._relprecs))]
879        relprec[0] = self._relprecs[0]
880        if a == 0:
881            zzpoly = self._poly.parent()(0)
882        else:
883            zzpoly = self._poly.rescale(Integer(a))
884        return Polynomial_padic_capped_relative_dense(self.parent(), (zzpoly, self._valbase, relprec, False, valadded, None), construct = True)
885
886    def quo_rem(self, right):
887        return self._quo_rem_naive(right)
888
889    def _quo_rem_naive(self, right):
890        """
891        An implementation of quo_rem that doesn't have good run-time or precision characteristics.
892        """
893        K = self.base_ring().fraction_field()
894        f = self.base_extend(K)
895        g = right.base_extend(K)
896        if g == 0:
897            raise ZeroDivisionError, "cannot divide by a polynomial indistinguishable from 0"
898        x = f.parent().gen()
899        quo = f.parent()(0)
900        while f.degree() >= g.degree():
901            a = f.leading_coefficient() / g.leading_coefficient()
902            quo = quo + a * (x ** (f.degree() - g.degree()))
903            f = f - a * (x ** (f.degree() - g.degree())) * g
904        return (quo, f)
905
906    #def gcd(self, right):
907    #    raise NotImplementedError
908
909    #def lcm(self, right):
910    #    raise NotImplementedError
911
912    def xgcd(self, right):
913        return self._xgcd(right)
914
915    #def discriminant(self):
916    #    raise NotImplementedError
917
918    def disc(self):
919        return self.discriminant()
920
921    #def resultant(self):
922    #    raise NotImplementedError
923
924    def newton_slopes(self):
925        """
926        Returns a list of the Newton slopes of this polynomial.  These are the valuations of the roots of this polynomial.
927
928        EXAMPLES:
929        sage: K = Qp(13)
930        sage: R.<t> = K[]
931        sage: f = t^4 + 13^5*t^2 + 4*13^2*t - 13^7
932        sage: f.newton_polygon()
933        [(0, 7), (1, 2), (4, 0)]
934        sage: f.newton_slopes()
935        [5, 2/3, 2/3, 2/3]
936        """
937        polygon = self.newton_polygon()
938        if polygon == []:
939            return []
940        answer = [infinity] * polygon[0][0]
941        for m in range(1, len(polygon)):
942            dx = polygon[m][0] - polygon[m - 1][0]
943            dy = polygon[m][1] - polygon[m - 1][1]
944            answer.extend([-dy / dx] * dx)
945        return answer       
946
947    def newton_polygon(self):
948        r"""
949        Returns a list of vertices of the Newton polygon of this polynomial.
950
951        NOTES:
952        The vertices are listed so that the first coordinates are strictly increasing, up to the polynomial's degree (not the limit of available precision information).  Also note that if some coefficients have very low precision an error is raised.
953
954        EXAMPLES:
955        sage: K = Qp(13)
956        sage: R.<t> = K[]
957        sage: f = t^4 + 13^5*t^2 + 4*13^2*t - 13^7
958        sage: f.newton_polygon()
959        [(0, 7), (1, 2), (4, 0)]
960        """
961        if self._poly == 0:
962            return []
963        for x in range(len(self._relprecs)):
964            if not self._relprecs[x] is infinity:
965                break
966        if self._valaddeds is None:
967            self._comp_valaddeds()
968        if self._poly[x] == 0:
969            raise PrecisionError, "first term with non-infinite valuation must have determined valuation"
970        yrel = self._valaddeds[x]
971        answer = [(x, self._valbase + yrel)]
972        xf = self._poly.degree()
973        if xf == x:
974            return answer
975        yfrel = self._valaddeds[xf]
976        curslope = (yfrel - yrel) / (xf - x)
977        for i in range(x + 1, xf):
978            yrel += curslope
979            if self._valaddeds[i] < yrel:
980                if self._relprecs[i] == self._valaddeds[i]:
981                    raise PrecisionError, "not enough precision known in coefficient %s to compute newton polygon"%i
982                yrel = self._valaddeds[i]
983                answer.append((i, self._valbase + yrel))
984                curslope = (yfrel - yrel) / (xf - i)
985        answer.append((xf, self._valbase + self._valaddeds[xf]))
986        return answer
987
988    def hensel_lift(self, a):
989        raise NotImplementedError
990
991    def factor_mod(self):
992        r"""
993        Returns the factorization of self modulo p.
994        """
995        self._normalize()
996        if self._valbase < 0:
997            raise ValueError, "Polynomial does not have integral coefficients"
998        elif self._valbase > 0:
999            raise ValueError, "Factorization of the zero polynomial not defined"
1000        elif min(self._relprecs) <= 0:
1001            raise PrecisionError, "Polynomial is not known to high enough precision"
1002        return self._poly.factor_mod(self.base_ring().prime())
1003
1004    def factor(self):
1005        # This will eventually be improved.
1006        if self == 0:
1007            raise ValueError, "Factorization of the zero polynomial not defined"
1008        from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
1009        from sage.rings.padics.factory import ZpCA
1010        base = self.base_ring()
1011        #print self.list()
1012        m = min([x.precision_absolute() for x in self.list()])
1013        #print m
1014        R = ZpCA(base.prime(), prec = m)
1015        S = PolynomialRing(R, self.parent().variable_name())
1016        F = S(self).factor()
1017        return Factorization([(self.parent()(a), b) for (a, b) in F], base(F.unit()))
1018
1019def _extend_by_infinity(L, n):
1020    return L + [infinity] * (n - len(L))
1021
1022def make_padic_poly(parent, x, version):
1023    if version == 0:
1024        return parent(x, construct = True)
1025    else:
1026        raise ValueError, "unknown pickling version"
Note: See TracBrowser for help on using the repository browser.