source: sage/rings/polynomial/polynomial_element_generic.py @ 6878:5fda14d7b30f

Revision 6878:5fda14d7b30f, 36.1 KB checked in by William Stein <wstein@…>, 6 years ago (diff)

merge
Line 
1"""
2Univariate Polynomials
3
4AUTHORS:
5    -- William Stein: first version
6    -- Martin Albrecht: Added singular coercion.
7    -- David Harvey: split off polynomial_integer_dense_ntl.pyx (2007-09)
8    -- Robert Bradshaw: split off polynomial_modn_dense_ntl.pyx (2007-09)
9
10TESTS:
11
12We test coercion in a particularly complicated situation:
13    sage: W.<w>=QQ['w']
14    sage: WZ.<z>=W['z']
15    sage: m = matrix(WZ,2,2,[1,z,z,z^2])
16    sage: a = m.charpoly()
17    sage: R.<x> = WZ[]
18    sage: R(a)
19    x^2 + ((-1)*z^2 + -1)*x
20"""
21
22################################################################################
23#       Copyright (C) 2007 William Stein <wstein@gmail.com>
24#
25#  Distributed under the terms of the GNU General Public License (GPL)
26#
27#                  http://www.gnu.org/licenses/
28################################################################################
29
30import copy
31
32from sage.rings.polynomial.polynomial_element import Polynomial, is_Polynomial, Polynomial_generic_dense
33from sage.structure.element import (IntegralDomainElement, EuclideanDomainElement,
34                                    PrincipalIdealDomainElement)
35
36from sage.rings.polynomial.polynomial_singular_interface import Polynomial_singular_repr
37
38from sage.libs.all import pari, pari_gen
39from sage.libs.ntl.all import ZZ as ntl_ZZ, ZZX, zero_ZZX
40from sage.structure.factorization import Factorization
41
42from sage.rings.infinity import infinity
43from sage.rings.rational_field import QQ
44from sage.rings.integer_ring import ZZ
45import sage.rings.integer as integer
46import sage.rings.complex_field as complex_field
47import sage.rings.arith as arith
48
49import sage.rings.fraction_field_element as fraction_field_element
50import sage.rings.polynomial.polynomial_ring
51
52               
53class Polynomial_generic_sparse(Polynomial):
54    """
55    A generic sparse polynomial.
56
57    EXAMPLES:
58        sage: R.<x> = PolynomialRing(PolynomialRing(QQ, 'y'), sparse=True)
59        sage: f = x^3 - x + 17
60        sage: type(f)
61        <class 'sage.rings.polynomial.polynomial_element_generic.Polynomial_generic_sparse'>
62        sage: loads(f.dumps()) == f
63        True
64    """
65    def __init__(self, parent, x=None, check=True, is_gen=False, construct=False):
66        Polynomial.__init__(self, parent, is_gen=is_gen)
67        if x is None:
68            self.__coeffs = {}
69            return
70        R = parent.base_ring()
71        if isinstance(x, Polynomial):
72            if x.parent() == self.parent():
73                x = dict(x.dict())
74            elif x.parent() == R:
75                x = {0:x}
76            else:
77                w = {}
78                for n, c in x.dict().iteritems():
79                    w[n] = R(c)
80                #raise TypeError, "Cannot coerce %s into %s."%(x, parent)
81        elif isinstance(x, list):
82            y = {}
83            for i in xrange(len(x)):
84                if x[i] != 0:
85                    y[i] = x[i]
86            x = y
87        elif not isinstance(x, dict):
88            x = {0:x}   # constant polynomials
89        elif isinstance(x, pari_gen):
90            x = [R(w) for w in x.Vecrev()]
91            check = True
92        if check:
93            self.__coeffs = {}
94            for i, z in x.iteritems():
95                self.__coeffs[i] = R(z)
96        else:
97            self.__coeffs = x
98        if check:
99            self.__normalize()
100
101    def dict(self):
102        """
103        Return a new copy of the dict of the underlying
104        elements of self.
105
106        EXAMPLES:
107            sage: R.<w> = PolynomialRing(Integers(8), sparse=True)
108            sage: f = 5 + w^1997 - w^10000; f
109            7*w^10000 + w^1997 + 5
110            sage: d = f.dict(); d
111            {0: 5, 10000: 7, 1997: 1}
112            sage: d[0] = 10
113            sage: f.dict()
114            {0: 5, 10000: 7, 1997: 1}           
115        """
116        return dict(self.__coeffs)
117
118    def valuation(self):
119        """
120        EXAMPLES:
121            sage: R.<w> = PolynomialRing(GF(9,'a'), sparse=True)
122            sage: f = w^1997 - w^10000
123            sage: f.valuation()
124            1997
125            sage: R(19).valuation()
126            0
127            sage: R(0).valuation()
128            +Infinity       
129        """
130        c = self.__coeffs.keys()
131        if len(c) == 0:
132            return infinity
133        return ZZ(min(self.__coeffs.keys()))
134
135    def derivative(self):
136        """
137        EXAMPLES:
138            sage: R.<w> = PolynomialRing(ZZ, sparse=True)
139            sage: f = R(range(9)); f
140            8*w^8 + 7*w^7 + 6*w^6 + 5*w^5 + 4*w^4 + 3*w^3 + 2*w^2 + w
141            sage: f.derivative()
142            64*w^7 + 49*w^6 + 36*w^5 + 25*w^4 + 16*w^3 + 9*w^2 + 4*w + 1
143        """
144        d = {}
145        for n, c in self.__coeffs.iteritems():
146            d[n-1] = n*c
147        if d.has_key(-1):
148            del d[-1]
149        return self.polynomial(d)
150
151    def _dict_unsafe(self):
152        """
153        Return unsafe access to the underlying dictionary of coefficients.
154
155        ** DO NOT use this, unless you really really know what you are doing. **
156
157        EXAMPLES:
158            sage: R.<w> = PolynomialRing(ZZ, sparse=True)
159            sage: f = w^15 - w*3; f
160            w^15 - 3*w
161            sage: d = f._dict_unsafe(); d
162            {1: -3, 15: 1}
163            sage: d[1] = 10; f
164            w^15 + 10*w       
165        """
166        return self.__coeffs
167
168    def _repr(self, name=None):
169        r"""
170        EXAMPLES:
171            sage: R.<w> = PolynomialRing(CDF, sparse=True)
172            sage: f = CDF(1,2) + w^5 - CDF(pi)*w + CDF(e)
173            sage: f._repr()
174            '1.0*w^5 + (-3.14159265359)*w + 3.71828182846 + 2.0*I'
175            sage: f._repr(name='z')
176            '1.0*z^5 + (-3.14159265359)*z + 3.71828182846 + 2.0*I'       
177
178        AUTHOR:
179            -- David Harvey (2006-08-05), based on Polynomial._repr()
180        """
181        s = " "
182        m = self.degree() + 1
183        if name is None:
184            name = self.parent().variable_name()
185        atomic_repr = self.parent().base_ring().is_atomic_repr()
186        coeffs = list(self.__coeffs.iteritems())
187        coeffs.sort()
188        for (n, x) in reversed(coeffs):
189            if x != 0:
190                if n != m-1:
191                    s += " + "
192                x = repr(x)
193                if not atomic_repr and n > 0 and (x.find("+") != -1 or x.find("-") != -1):
194                    x = "(%s)"%x
195                if n > 1:
196                    var = "*%s^%s"%(name,n)
197                elif n==1:
198                    var = "*%s"%name
199                else:
200                    var = ""
201                s += "%s%s"%(x,var)
202        if atomic_repr:
203            s = s.replace(" + -", " - ")
204        s = s.replace(" 1*"," ")
205        s = s.replace(" -1*", " -")
206        if s==" ":
207            return "0"
208        return s[1:]
209
210    def __normalize(self):
211        x = self.__coeffs
212        zero = self.base_ring()(0)
213        D = [n for n, z in x.iteritems() if not z]
214        for n in D:
215            del x[n]
216       
217    def __getitem__(self,n):
218        """
219        Return the n-th coefficient of this polynomial.
220
221        Negative indexes are allowed and always return 0 (so you can
222        view the polynomial as embedding Laurent series).
223
224        EXAMPLES:
225            sage: R.<w> = PolynomialRing(RDF, sparse=True)
226            sage: e = RDF(e)
227            sage: f = sum(e^n*w^n for n in range(4)); f
228            20.0855369232*w^3 + 7.38905609893*w^2 + 2.71828182846*w + 1.0
229            sage: f[1]
230            2.71828182846
231            sage: f[5]
232            0.0
233            sage: f[-1]
234            0.0           
235        """
236        if not self.__coeffs.has_key(n):
237            return self.base_ring()(0)
238        return self.__coeffs[n]
239
240    def __getslice__(self, i, j):
241        """
242        EXAMPLES:
243            sage: R.<x> = PolynomialRing(RealField(19), sparse=True)
244            sage: f = (2-3.5*x)^3; f
245            -42.875*x^3 + 73.500*x^2 - 42.000*x + 8.0000
246            sage: f[1:3]
247            73.500*x^2 - 42.000*x
248            sage: f[:2]
249            -42.000*x + 8.0000
250            sage: f[2:]
251            -42.875*x^3 + 73.500*x^2       
252        """
253        if i < 0:
254            i = 0
255        v = {}
256        x = self.__coeffs
257        for k in x.keys():
258            if i <= k and k < j:
259                v[k] = x[k]
260        P = self.parent()
261        return P(v)
262
263    def _unsafe_mutate(self, n, value):
264        r"""
265        Change the coefficient of $x^n$ to value.
266
267        ** NEVER USE THIS ** -- unless you really know what you are doing.
268
269        EXAMPLES:
270            sage: R.<z> = PolynomialRing(CC, sparse=True)
271            sage: f = z^2 + CC.0; f
272            1.00000000000000*z^2 + 1.00000000000000*I
273            sage: f._unsafe_mutate(0, 10)
274            sage: f
275            1.00000000000000*z^2 + 10.0000000000000
276
277        Much more nasty:
278            sage: z._unsafe_mutate(1, 0)
279            sage: z
280            0       
281        """
282        n = int(n)
283        value = self.base_ring()(value)
284        x = self.__coeffs
285        if n < 0:
286            raise IndexError, "polynomial coefficient index must be nonnegative"
287        if value == 0:
288            if x.has_key(n):
289                del x[n]
290        else:
291            x[n] = value
292           
293    def list(self):
294        """
295        Return a new copy of the list of the underlying
296        elements of self.
297
298        EXAMPLES:
299            sage: R.<z> = PolynomialRing(Integers(100), sparse=True)
300            sage: f = 13*z^5 + 15*z^2 + 17*z
301            sage: f.list()
302            [0, 17, 15, 0, 0, 13]       
303        """
304        zero = self.base_ring()(0)
305        v = [zero] * (self.degree()+1)
306        for n, x in self.__coeffs.iteritems():
307            v[n] = x
308        return v
309
310    #def _pari_(self, variable=None):
311    #    if variable is None:
312    #        return self.__pari
313    #    else:
314    #        return self.__pari.subst('x',variable)
315
316    def degree(self):
317        """
318        Return the degree of this sparse polynomial.
319       
320        EXAMPLES:
321            sage: R.<z> = PolynomialRing(ZZ, sparse=True)
322            sage: f = 13*z^50000 + 15*z^2 + 17*z
323            sage: f.degree()
324            50000       
325        """
326        v = self.__coeffs.keys()
327        if len(v) == 0:
328            return -1
329        return max(v)
330
331    def _add_(self, right):
332        r"""
333        EXAMPLES:
334            sage: R.<x> = PolynomialRing(Integers(), sparse=True)
335            sage: (x^100000 + 2*x^50000) + (4*x^75000 - 2*x^50000 + 3*x)
336            x^100000 + 4*x^75000 + 3*x
337
338        AUTHOR:
339            -- David Harvey (2006-08-05)
340        """
341        output = dict(self.__coeffs)
342
343        for (index, coeff) in right.__coeffs.iteritems():
344            if index in output:
345                output[index] += coeff
346            else:
347                output[index] = coeff
348
349        output = self.polynomial(output, check=False)
350        output.__normalize()
351        return output
352
353    def _mul_(self, right):
354        r"""
355        EXAMPLES:
356            sage: R.<x> = PolynomialRing(ZZ, sparse=True)
357            sage: (x^100000 - x^50000) * (x^100000 + x^50000)
358             x^200000 - x^100000
359            sage: (x^100000 - x^50000) * R(0)
360             0
361
362        AUTHOR:
363            -- David Harvey (2006-08-05)
364        """
365        output = {}
366
367        for (index1, coeff1) in self.__coeffs.iteritems():
368            for (index2, coeff2) in right.__coeffs.iteritems():
369                product = coeff1 * coeff2
370                index = index1 + index2
371                if index in output:
372                    output[index] += product
373                else:
374                    output[index] = product
375
376        output = self.polynomial(output, check=False)
377        output.__normalize()
378        return output
379
380    def shift(self, n):
381        r"""
382        Returns this polynomial multiplied by the power $x^n$. If $n$ is negative,
383        terms below $x^n$ will be discarded. Does not change this polynomial.
384
385        EXAMPLES:
386            sage: R.<x> = PolynomialRing(ZZ, sparse=True)
387            sage: p = x^100000 + 2*x + 4
388            sage: type(p)
389            <class 'sage.rings.polynomial.polynomial_element_generic.Polynomial_generic_sparse'>
390            sage: p.shift(0)
391             x^100000 + 2*x + 4
392            sage: p.shift(-1)
393             x^99999 + 2
394            sage: p.shift(-100002)
395             0
396            sage: p.shift(2)
397             x^100002 + 2*x^3 + 4*x^2
398
399        AUTHOR:
400            -- David Harvey (2006-08-06)
401        """
402        n = int(n)
403        if n == 0:
404            return self
405        if n > 0:
406            output = {}
407            for (index, coeff) in self.__coeffs.iteritems():
408                output[index + n] = coeff
409            return self.polynomial(output, check=False)
410        if n < 0:
411            output = {}
412            for (index, coeff) in self.__coeffs.iteritems():
413                if index + n >= 0:
414                    output[index + n] = coeff
415            return self.polynomial(output, check=False)
416
417
418class Polynomial_generic_domain(Polynomial, IntegralDomainElement):
419    def __init__(self, parent, is_gen=False, construct=False):
420        Polynomial.__init__(self, parent, is_gen=is_gen)
421       
422    def is_unit(self):
423        r"""
424        Return True if this polynomial is a unit.
425
426        EXERCISE (Atiyah-McDonald, Ch 1): Let $A[x]$ be a polynomial
427        ring in one variable.  Then $f=\sum a_i x^i \in A[x]$ is a
428        unit if and only if $a_0$ is a unit and $a_1,\ldots, a_n$ are
429        nilpotent.
430
431        EXAMPLES:
432            sage: R.<z> = PolynomialRing(ZZ, sparse=True)
433            sage: (2 + z^3).is_unit()
434            False
435            sage: f = -1 + 3*z^3; f
436            3*z^3 - 1
437            sage: f.is_unit()
438            False
439            sage: R(-3).is_unit()
440            False
441            sage: R(-1).is_unit()
442            True
443            sage: R(0).is_unit()
444            False       
445        """
446        if self.degree() > 0:
447            return False
448        return self[0].is_unit()
449
450class Polynomial_generic_field(Polynomial_singular_repr,
451                               Polynomial_generic_domain,
452                               EuclideanDomainElement):
453
454    def quo_rem(self, other):
455        """
456        Returns a tuple (quotient, remainder) where
457            self = quotient*other + remainder.
458
459        EXAMPLES:
460            sage: R.<y> = PolynomialRing(QQ)
461            sage: K.<t> = NumberField(y^2 - 2)
462            sage: P.<x> = PolynomialRing(K)
463            sage: x.quo_rem(K(1))
464            (x, 0)
465            sage: x.xgcd(K(1))
466            (1, 0, 1)
467        """
468        other = self.parent()(other)
469        if other.is_zero(): 
470            raise ZeroDivisionError, "other must be nonzero"
471           
472        # This is algorithm 3.1.1 in Cohen GTM 138
473        A = self
474        B = other
475        R = A
476        Q = self.polynomial(0)
477        X = self.parent().gen()
478        while R.degree() >= B.degree():
479            aaa = (R.leading_coefficient()/B.leading_coefficient())
480            bbb = X**(R.degree()-B.degree())
481            S = aaa * bbb
482            Q += S
483            R -= S*B           
484        return (Q, R)
485
486    def _gcd(self, other):
487        """
488        Return the GCD of self and other, as a monic polynomial.
489        """
490        g = EuclideanDomainElement._gcd(self, other)
491        c = g.leading_coefficient()
492        if c.is_unit():
493            return (1/c)*g
494        return g
495       
496       
497class Polynomial_generic_sparse_field(Polynomial_generic_sparse, Polynomial_generic_field):
498    """
499    EXAMPLES:
500        sage: R.<x> = PolynomialRing(RR, sparse=True)
501        sage: f = x^3 - x + 17
502        sage: type(f)
503        <class 'sage.rings.polynomial.polynomial_element_generic.Polynomial_generic_sparse_field'>
504        sage: loads(f.dumps()) == f
505        True
506    """
507    def __init__(self, parent, x=None, check=True, is_gen = False, construct=False):
508        Polynomial_generic_sparse.__init__(self, parent, x, check, is_gen)
509
510
511class Polynomial_generic_dense_field(Polynomial_generic_dense, Polynomial_generic_field):
512    def __init__(self, parent, x=None, check=True, is_gen = False, construct=False):
513        Polynomial_generic_dense.__init__(self, parent, x, check, is_gen)
514
515
516class Polynomial_rational_dense(Polynomial_generic_field):
517    """
518    A dense polynomial over the rational numbers.
519    """
520    def __init__(self, parent, x=None, check=True, is_gen=False, construct=False):
521        Polynomial.__init__(self, parent, is_gen=is_gen)
522
523        if construct:
524            self.__poly = x
525            return
526
527        self.__poly = pari([]).Polrev()
528
529        if x is None:
530            return         # leave initialized to 0 polynomial.
531
532
533        if fraction_field_element.is_FractionFieldElement(x):
534            if not x.denominator().is_unit():
535                raise TypeError, "denominator must be a unit"
536            elif x.denominator() != 1:
537                x = x.numerator() * x.denominator().inverse_of_unit()
538            else:
539                x = x.numerator()
540       
541        if isinstance(x, Polynomial):
542            if x.parent() == self.parent():
543                self.__poly = x.__poly.copy()
544                return
545            else:
546                x = [QQ(a) for a in x.list()]
547                check = False
548
549        elif isinstance(x, dict):
550            x = self._dict_to_list(x, QQ(0))
551           
552           
553        elif isinstance(x, pari_gen):
554            f = x.Polrev()
555            self.__poly = f
556            assert self.__poly.type() == "t_POL"           
557            return
558           
559        elif not isinstance(x, list):
560            x = [x]   # constant polynomials
561
562        if check:
563            x = [QQ(z) for z in x]
564
565        self.__list = list(x)
566        while len(self.__list) > 0 and not self.__list[-1]:
567            del self.__list[-1]
568       
569        # NOTE: It is much faster to convert to string and let pari's parser at it,
570        # which is why we pass str(x) in.
571        self.__poly = pari(str(x)).Polrev()
572        assert self.__poly.type() == "t_POL"
573
574    def _repr(self, name=None):
575        if name is None:
576            name = self.parent().variable_name()
577        return str(self.__poly).replace("x", name)
578
579    def _repr_(self):
580        return self._repr()
581
582    def __reduce__(self):
583        return Polynomial_rational_dense, \
584               (self.parent(), self.list(), False, self.is_gen())
585       
586    def __getitem__(self, n):
587        return QQ(self.__poly[n])
588
589    def __getslice__(self, i, j):
590        if i < 0:
591            i = 0
592        v = [QQ(x) for x in self.__poly[i:j]]
593        P = self.parent()
594        return P([0]*int(i) + v)
595
596    def _pow(self, n):
597        if self.degree() <= 0:
598            return self.parent()(self[0]**n)
599        if n < 0:
600            return (~self)**(-n)
601        return Polynomial_rational_dense(self.parent(), self.__poly**n, construct=True)
602     
603    def _add_(self, right):
604        return Polynomial_rational_dense(self.parent(), 
605                                         self.__poly + right.__poly, construct=True)
606
607    def is_irreducible(self):
608        """
609        EXAMPLES:
610            sage: R.<x> = QQ[]
611            sage: (x^2 + 2).is_irreducible()
612            True
613            sage: (x^2 - 1).is_irreducible()
614            False
615        """
616        try:
617            return self.__poly.polisirreducible()
618        except NotImplementedError:
619            F = self.__poly.factor()
620            if len(F) > 1 or F[0][1] > 1:
621                return False
622            return True
623
624    def galois_group(self, pari_group=False, use_kash=False):
625        r"""
626        Return the Galois group of f as a permutation group.
627
628        INPUT:
629            self -- an irreducible polynomial
630           
631            pari_group -- bool (default: False); if True instead return
632                          the Galois group as a PARI group.  This has
633                          a useful label in it, and may be slightly faster
634                          since it doesn't require looking up a group in
635                          Gap.  To get a permutation group from a PARI
636                          group P, type PermutationGroup(P).
637                         
638            use_kash --   bool (default: False); if True use KASH's Galois
639                          command instead of using the PARI C library.
640                          An attempt is always made to use KASH if the
641                          degree of the polynomial is >= 12.
642
643        ALGORITHM: The Galois group is computed using PARI in C
644        library mode, or possibly kash if available.
645
646        \note{ The PARI documentation contains the following warning:
647        The method used is that of resolvent polynomials and is
648        sensitive to the current precision. The precision is updated
649        internally but, in very rare cases, a wrong result may be
650        returned if the initial precision was not sufficient.}
651
652        EXAMPLES:
653            sage: R.<x> = PolynomialRing(QQ)
654            sage: f = x^4 - 17*x^3 - 2*x + 1
655            sage: G = f.galois_group(); G            # uses optional database_gap package
656            Transitive group number 5 of degree 4
657            sage: G.gens()                           # uses optional database_gap package
658            ((1,2,3,4), (1,2))
659            sage: G.order()                          # uses optional database_gap package
660            24
661
662        It is potentially useful to instead obtain the corresponding
663        PARI group, which is little more than a $4$-tuple.  See the
664        PARI manual for the exact details.  (Note that the third
665        entry in the tuple is in the new standard ordering.)
666            sage: f = x^4 - 17*x^3 - 2*x + 1
667            sage: G = f.galois_group(pari_group=True); G
668            PARI group [24, -1, 5, "S4"] of degree 4
669            sage: PermutationGroup(G)                # uses optional database_gap package
670            Transitive group number 5 of degree 4
671
672        You can use KASH to compute Galois groups as well.  The
673        avantage is that KASH can compute Galois groups of fields up
674        to degree 23, whereas PARI only goes to degree 11.  (In my
675        not-so-thorough experiments PARI is faster than KASH.)
676
677            sage: f = x^4 - 17*x^3 - 2*x + 1
678            sage: f.galois_group(use_kash=true)      # requires optional KASH
679            Transitive group number 5 of degree 4
680       
681        """
682        from sage.groups.all import PariGroup, PermutationGroup, TransitiveGroup
683        if not self.is_irreducible():
684            raise ValueError, "polynomial must be irreducible"
685        if self.degree() > 11 or use_kash:
686            # TODO -- maybe use KASH if available or print message that user should install KASH?
687            try:
688                from sage.interfaces.all import kash
689                kash.eval('X := PolynomialRing(RationalField()).1')
690                s = self._repr(name='X')
691                G = kash('Galois(%s)'%s)
692                d = int(kash.eval('%s.ext1'%G.name()))
693                n = int(kash.eval('%s.ext2'%G.name()))
694                return TransitiveGroup(d, n)
695            except RuntimeError:
696                raise NotImplementedError, "Sorry, computation of Galois groups of fields of degree bigger than 11 is not yet implemented.  Try installing the optional free (closed source) KASH package, which supports up to degree $23$."
697        G = self.__poly.polgalois()
698        H = PariGroup(G, self.degree())
699        if pari_group:
700            return H
701        else:
702            return PermutationGroup(H)
703
704    def quo_rem(self, right):
705        """
706        Returns a tuple (quotient, remainder) where
707            self = quotient*right + remainder.
708
709        EXAMPLES:
710            sage: R.<x> = QQ[]
711            sage: f = x^5 + 17*x + 3
712            sage: g = x^3 - 19
713            sage: q,r = f.quo_rem(g)
714            sage: q*g + r
715            x^5 + 17*x + 3       
716        """
717        if not isinstance(right, Polynomial_rational_dense):
718            right = self.parent()(right)
719        if right.parent() != self.parent():
720            raise TypeError
721        v = self.__poly.divrem(right.__poly)
722        return Polynomial_rational_dense(self.parent(), v[0], construct=True), \
723               Polynomial_rational_dense(self.parent(), v[1], construct=True)
724       
725                                   
726    def _mul_(self, right):
727        """
728        EXAMPLES:
729            sage: R.<x> = QQ[]
730            sage: (x - QQ('2/3'))*(x^2 - 8*x + 16)
731            x^3 - 26/3*x^2 + 64/3*x - 32/3
732        """
733        return self.parent()(self.__poly * right.__poly, construct=True)       
734
735    def _sub_(self, right):
736        """
737        EXAMPLES:
738            sage: R.<x> = QQ[]
739            sage: x^5 + 17*x^3 + x+ 3 - (x^3 - 19)
740            x^5 + 16*x^3 + x + 22
741        """
742        return self.parent()(self.__poly - right.__poly, construct=True)
743                               
744    def _unsafe_mutate(self, n, value):
745        try:
746            del self.__list
747        except AttributeError:
748            pass
749        n = int(n)
750        if n < 0:
751            raise IndexError, "n must be >= 0"
752        if n <= self.__poly.poldegree():
753            self.__poly[n] = QQ(value)
754        else:
755            self.__poly = self.__poly + pari('(%s)*x^%s'%(QQ(value),n))
756        if hasattr(self, "__list"):
757            del self.__list
758     
759    def complex_roots(self, flag=0):
760        """
761        Returns the complex roots of this polynomial.
762        INPUT:
763            flag -- optional, and can be
764                    0: (default), uses Schonhage's method modified by Gourdon,
765                    1: uses a modified Newton method.
766        OUTPUT:
767            list of complex roots of this polynomial, counted with multiplicities.
768           
769        NOTE: Calls the pari function polroots.
770
771        EXAMPLE:
772        We compute the roots of the characteristic polynomial of some Salem numbers:
773            sage: R.<x> = PolynomialRing(QQ)
774            sage: f = 1 - x^2 - x^3 - x^4 + x^6
775            sage: f.complex_roots()[0]
776            0.713639173536901
777        """
778        R = self.__poly.polroots(flag)
779        C = complex_field.ComplexField()
780        return [C(a) for a in R]
781
782    def real_root_intervals(self):
783        """
784        Returns isolating intervals for the real roots of this polynomial.
785
786        EXAMPLE:
787            sage: R.<x> = PolynomialRing(QQ)
788            sage: f = (x - 1/2) * (x - 3/4) * (x - 3/2)
789            sage: f.real_root_intervals()
790            [((243/512, 1215/2048), 1), ((729/1024, 1701/2048), 1), ((243/256, 1011/512), 1)]
791        """
792
793        from sage.rings.polynomial.real_roots import real_roots
794
795        return real_roots(self)
796
797    def copy(self):
798        """
799        Return a copy of this polynomial.
800        """
801        f = Polynomial_rational_dense(self.parent())
802        f.__poly = self.__poly.copy()
803        return f
804       
805    def degree(self):
806        """
807        Return the degree of this polynomial.  The zero polynomial
808        has degree -1.
809
810        EXAMPLES:
811            sage: R.<x> = QQ[]
812            sage: (x^5 + 17*x^3 + x+ 3).degree()
813            5
814            sage: R(0).degree()
815            -1
816            sage: type(x.degree())
817            <type 'sage.rings.integer.Integer'>
818        """
819        return ZZ(max(self.__poly.poldegree(), -1))
820
821    def discriminant(self):
822        """
823        EXAMPLES:
824            sage: _.<x> = PolynomialRing(QQ)
825            sage: f = x^3 + 3*x - 17
826            sage: f.discriminant()
827            -7911
828        """
829        return QQ(self.__poly.poldisc())
830
831    def disc(self):
832        """
833        Same as discriminant().
834        """
835        return self.discriminant()
836     
837    def factor_mod(self, p):
838        """
839        Return the factorization of self modulo the prime p.
840
841        INPUT:
842            p -- prime
843
844        OUTPUT:
845            factorization of self reduced modulo p.
846
847        EXAMPLES:
848            sage: R.<x> = QQ[]
849            sage: (x^5 + 17*x^3 + x+ 3).factor_mod(3)
850            x * (x^2 + 1)^2
851            sage: (x^5 + 2).factor_mod(5)
852            (x + 2)^5       
853        """
854        import sage.rings.finite_field as finite_field
855        p = integer.Integer(p)
856        if not p.is_prime():
857            raise ValueError, "p must be prime"
858        G = self._pari_().factormod(p)
859        K = finite_field.FiniteField(p)
860        R = sage.rings.polynomial.polynomial_ring.PolynomialRing(K, names=self.parent().variable_name())
861        return R(1)._factor_pari_helper(G, unit=R(self).leading_coefficient())
862
863    def factor_padic(self, p, prec=10):
864        """
865        Return p-adic factorization of self to given precision.
866
867        INPUT:
868            p -- prime
869            prec -- integer; the precision
870
871        OUTPUT:
872            factorization of self viewed as a polynomial over the p-adics
873
874        EXAMPLES:
875            sage: R.<x> = QQ[]
876            sage: f = x^3 - 2
877            sage.: f.factor_padic(2)
878            (1 + O(2^10))*x^3 + O(2^10)*x^2 + O(2^10)*x + (2 + 2^2 + 2^3 + 2^4 + 2^5 + 2^6 + 2^7 + 2^8 + 2^9 + O(2^10))
879            sage.: f.factor_padic(3)
880            (1 + O(3^10))*x^3 + O(3^10)*x^2 + O(3^10)*x + (1 + 2*3 + 2*3^2 + 2*3^3 + 2*3^4 + 2*3^5 + 2*3^6 + 2*3^7 + 2*3^8 + 2*3^9 + O(3^10))
881            sage.: f.factor_padic(5)
882            ((1 + O(5^10))*x + (2 + 4*5 + 2*5^2 + 2*5^3 + 5^4 + 3*5^5 + 4*5^7 + 2*5^8 + 5^9 + O(5^10))) * ((1 + O(5^10))*x^2 + (3 + 2*5^2 + 2*5^3 + 3*5^4 + 5^5 + 4*5^6 + 2*5^8 + 3*5^9 + O(5^10))*x + (4 + 5 + 2*5^2 + 4*5^3 + 4*5^4 + 3*5^5 + 3*5^6 + 4*5^7 + 4*5^9 + O(5^10)))
883        """
884        from sage.rings.padics.factory import Qp
885        p = integer.Integer(p)
886        if not p.is_prime():
887            raise ValueError, "p must be prime"
888        prec = integer.Integer(prec)
889        if prec <= 0:
890            raise ValueError, "prec must be positive"
891        K = Qp(p, prec, type='capped-rel')
892        R = sage.rings.polynomial.polynomial_ring.PolynomialRing(K, names=self.parent().variable_name())
893        return R(self).factor(absprec = prec)
894
895    def list(self):
896        """
897        Return a new copy of the list of the underlying
898        elements of self.
899
900        EXAMPLES:
901            sage: _.<x> = PolynomialRing(QQ)
902            sage: f = x^3 + 3*x - 17/13; f
903            x^3 + 3*x - 17/13
904            sage: v = f.list(); v
905            [-17/13, 3, 0, 1]
906            sage: v[0] = 0
907            sage: f
908            x^3 + 3*x - 17/13
909            sage: f.list()
910            [-17/13, 3, 0, 1]           
911        """
912        return [QQ(x) for x in self.__poly.Vecrev()]
913
914##     def partial_fraction(self, g):
915##         """
916##         Return partial fraction decomposition of self/g, where g
917##         has the same parent as self.
918##         """
919##         g = self.parent()(g)
920##         from sage.interfaces.maxima import maxima
921##         h = maxima(self)/maxima(g)
922##         k = h.partfrac(self.parent().variable())
923
924    def rescale(self, a):
925        """
926        Return f(a*X).
927        """
928        b = 1
929        c = []
930        for i in range(self.degree()+1):
931            c.append(b*self[i])
932            b *= a
933        return self.parent()(c)
934
935    def hensel_lift(self, p, e):
936        """
937        Assuming that self factors modulo $p$ into distinct factors,
938        computes the Hensel lifts of these factors modulo $p^e$.  We
939        assume that $p$ has integer coefficients.
940        """
941        p = integer.Integer(p)
942        if not p.is_prime():
943            raise ValueError, "p must be prime"
944        e = integer.Integer(e)
945        if e < 1:
946            raise ValueError, "e must be at least 1"
947        F = self.factor_mod(p)
948        y = []
949        for g, n in F:
950            if n > 1:
951                raise ArithmeticError, "The polynomial must be square free modulo p."
952            y.append(g)
953        H = self._pari_().polhensellift(y, p, e)
954        R = integer_mod_ring.IntegerModRing(p**e)
955        S = sage.rings.polynomial.polynomial_ring.PolynomialRing(R, self.parent().variable_name())
956        return [S(eval(str(m.Vec().Polrev().Vec()))) for m in H]
957
958class Polynomial_padic_generic_dense(Polynomial_generic_dense, Polynomial_generic_domain):
959    def __init__(self, parent, x=None, check=True, is_gen = False, construct=False, absprec=None):
960        Polynomial_generic_dense.__init__(self, parent, x, check, is_gen, absprec=absprec)
961
962    def _mul_(self, right):
963        return self._mul_generic(right)
964
965    def __pow__(self, right):
966        #computing f^p in this way loses precision
967        return self._pow(right)
968
969    def clear_zeros(self):
970        """
971        This function replaces coefficients of the polynomial that evaluate as equal to 0 with the zero element of the base ring that has the maximum possible precision.
972       
973        WARNING: this function mutates the underlying polynomial.
974        """
975        coeffs = self._Polynomial_generic_dense__coeffs
976        for n in range(len(coeffs)):
977            if not coeffs[n]:
978                self._Polynomial_generic_dense__coeffs[n] = self.base_ring()(0)
979
980    def _repr(self, name=None):
981        r"""
982        EXAMPLES:
983            sage: R.<w> = PolynomialRing(Zp(5, prec=5, type = 'capped-abs', print_mode = 'val-unit'))
984            sage: f = 24 + R(4/3)*w + w^4
985            sage: f._repr()
986            '(1 + O(5^5))*w^4 + (1043 + O(5^5))*w + (24 + O(5^5))'
987            sage: f._repr(name='z')
988            '(1 + O(5^5))*z^4 + (1043 + O(5^5))*z + (24 + O(5^5))'
989
990        AUTHOR:
991            -- David Roe (2007-03-03), based on Polynomial_generic_dense._repr()
992        """
993        s = " "
994        n = m = self.degree()
995        if name is None:
996            name = self.parent().variable_name()
997        atomic_repr = self.parent().base_ring().is_atomic_repr()
998        coeffs = self.list()
999        for x in reversed(coeffs):
1000            if x.valuation() != infinity:
1001                if n != m:
1002                    s += " + "
1003                x = repr(x)
1004                x = "(%s)"%x
1005                if n > 1:
1006                    var = "*%s^%s"%(name,n)
1007                elif n==1:
1008                    var = "*%s"%name
1009                else:
1010                    var = ""
1011                s += "%s%s"%(x,var)
1012            n -= 1
1013        if s==" ":
1014            return "0"
1015        return s[1:]
1016
1017
1018    def factor(self, absprec = None):
1019        if self == 0:
1020            raise ValueError, "Factorization of 0 not defined"
1021        if absprec is None:
1022            absprec = min([x.precision_absolute() for x in self.list()])
1023        else:
1024            absprec = integer.Integer(absprec)
1025        if absprec <= 0:
1026            raise ValueError, "absprec must be positive"
1027        G = self._pari_().factorpadic(self.base_ring().prime(), absprec)
1028        pols = G[0]
1029        exps = G[1]
1030        F = []
1031        R = self.parent()
1032        for i in xrange(len(pols)):
1033            f = R(pols[i], absprec = absprec)
1034            e = int(exps[i])
1035            F.append((f,e))
1036           
1037        if R.base_ring().is_field():
1038            # When the base ring is a field we normalize
1039            # the irreducible factors so they have leading
1040            # coefficient 1.
1041            for i in range(len(F)):
1042                cur = F[i][0].leading_coefficient()
1043                if cur != 1:
1044                    F[i] = (F[i][0].monic(), F[i][1])
1045            return Factorization(F, self.leading_coefficient())
1046        else:
1047            # When the base ring is not a field, we normalize
1048            # the irreducible factors so that the leading term
1049            # is a power of p.  We also ensure that the gcd of
1050            # the coefficients of each term is 1.
1051            c = self.leading_coefficient().valuation()
1052            u = self.base_ring()(1)
1053            for i in range(len(F)):
1054                upart = F[i][0].leading_coefficient().unit_part()
1055                lval = F[i][0].leading_coefficient().valuation()
1056                if upart != 1:
1057                    F[i] = (F[i][0] // upart, F[i][1])
1058                    u *= upart ** F[i][1]
1059                c -= lval
1060            if c != 0:
1061                F.append((self.parent()(self.base_ring().prime_pow(c)), 1))
1062            return Factorization(F, u)
1063
1064class Polynomial_padic_ring_dense(Polynomial_padic_generic_dense):
1065    def content(self):
1066        if self == 0:
1067            return self.base_ring()(0)
1068        else:
1069            return self.base_ring()(self.base_ring().prime_pow(min([x.valuation() for x in self.list()])))
1070
1071class Polynomial_padic_field_dense(Polynomial_padic_generic_dense, Polynomial_generic_dense_field):
1072    def content(self):
1073        if self != 0:
1074            return self.base_ring()(1)
1075        else:
1076            return self.base_ring()(0)
1077
1078    def _xgcd(self, other):
1079        H = Polynomial_generic_dense_field._xgcd(self, other)
1080        c = ~H[0].leading_coefficient()
1081        return c * H[0], c * H[1], c * H[2]
1082
1083
1084class Polynomial_padic_ring_lazy_dense(Polynomial_padic_ring_dense):
1085    pass
1086
1087class Polynomial_padic_field_lazy_dense(Polynomial_padic_field_dense):
1088    pass
Note: See TracBrowser for help on using the repository browser.