Ticket #10369: 10369_nffactor.patch

File 10369_nffactor.patch, 7.8 KB (added by jdemeyer, 11 years ago)
  • sage/rings/polynomial/polynomial_element.pyx

    # HG changeset patch
    # User Jeroen Demeyer <jdemeyer@cage.ugent.be>
    # Date 1292018604 -3600
    # Node ID 04ad692114783a6bb3664db8ad6cef90e8be42c4
    # Parent  93c88cef0de236ce3bd64cf117fff39306a943af
    #10369: Remove special case for factoring non-monic polynomials over number fields.
    Also, add doctests that the issues from ticket #10369 are fixed and rename
    "Pari" to "PARI".
    
    diff -r 93c88cef0de2 -r 04ad69211478 sage/rings/polynomial/polynomial_element.pyx
    a b  
    1818    True
    1919"""
    2020
    21 ################################################################################
     21#*****************************************************************************
    2222#       Copyright (C) 2007 William Stein <wstein@gmail.com>
    2323#
    2424#  Distributed under the terms of the GNU General Public License (GPL)
    25 #
     25#  as published by the Free Software Foundation; either version 2 of
     26#  the License, or (at your option) any later version.
    2627#                  http://www.gnu.org/licenses/
    27 ################################################################################
     28#*****************************************************************************
     29
    2830
    2931cdef is_FractionField, is_RealField, is_ComplexField
    3032cdef coerce_binop, generic_power, parent
     
    26502652            sage: pol = t^3 + (-4*a^3 + 2*a)*t^2 - 11/3*a^2*t + 2/3*a^3 - 4/3*a
    26512653            sage: pol.factor()
    26522654            (t - 2*a^3 + a) * (t - 4/3*a^3 + 2/3*a) * (t - 2/3*a^3 + 1/3*a)
     2655
     2656        Test that ticket #10369 is fixed::
     2657
     2658            sage: x = polygen(QQ)
     2659            sage: K.<a> = NumberField(x^6 + x^5 + x^4 + x^3 + x^2 + x + 1)
     2660            sage: R.<t> = PolynomialRing(K)
     2661           
     2662            sage: pol = (-1/7*a^5 - 1/7*a^4 - 1/7*a^3 - 1/7*a^2 - 2/7*a - 1/7)*t^10 + (4/7*a^5 - 2/7*a^4 - 2/7*a^3 - 2/7*a^2 - 2/7*a - 6/7)*t^9 + (90/49*a^5 + 152/49*a^4 + 18/49*a^3 + 24/49*a^2 + 30/49*a + 36/49)*t^8 + (-10/49*a^5 + 10/7*a^4 + 198/49*a^3 - 102/49*a^2 - 60/49*a - 26/49)*t^7 + (40/49*a^5 + 45/49*a^4 + 60/49*a^3 + 277/49*a^2 - 204/49*a - 78/49)*t^6 + (90/49*a^5 + 110/49*a^4 + 2*a^3 + 80/49*a^2 + 46/7*a - 30/7)*t^5 + (30/7*a^5 + 260/49*a^4 + 250/49*a^3 + 232/49*a^2 + 32/7*a + 8)*t^4 + (-184/49*a^5 - 58/49*a^4 - 52/49*a^3 - 66/49*a^2 - 72/49*a - 72/49)*t^3 + (18/49*a^5 - 32/49*a^4 + 10/49*a^3 + 4/49*a^2)*t^2 + (2/49*a^4 - 4/49*a^3 + 2/49*a^2)*t
     2663            sage: pol.factor()
     2664            (-1/7*a^5 - 1/7*a^4 - 1/7*a^3 - 1/7*a^2 - 2/7*a - 1/7) * t * (t - a^5 - a^4 - a^3 - a^2 - a - 1)^4 * (t^5 + (-12/7*a^5 - 10/7*a^4 - 8/7*a^3 - 6/7*a^2 - 4/7*a - 2/7)*t^4 + (12/7*a^5 - 8/7*a^3 + 16/7*a^2 + 2/7*a + 20/7)*t^3 + (-20/7*a^5 - 20/7*a^3 - 20/7*a^2 + 4/7*a - 2)*t^2 + (12/7*a^5 + 12/7*a^3 + 2/7*a + 16/7)*t - 4/7*a^5 - 4/7*a^3 - 4/7*a - 2/7)
     2665
     2666            sage: pol = (1/7*a^2 - 1/7*a)*t^10 + (4/7*a - 6/7)*t^9 + (102/49*a^5 + 99/49*a^4 + 96/49*a^3 + 93/49*a^2 + 90/49*a + 150/49)*t^8 + (-160/49*a^5 - 36/49*a^4 - 48/49*a^3 - 8/7*a^2 - 60/49*a - 60/49)*t^7 + (30/49*a^5 - 55/49*a^4 + 20/49*a^3 + 5/49*a^2)*t^6 + (6/49*a^4 - 12/49*a^3 + 6/49*a^2)*t^5
     2667            sage: pol.factor()
     2668            (1/7*a^2 - 1/7*a) * t^5 * (t^5 + (-40/7*a^5 - 38/7*a^4 - 36/7*a^3 - 34/7*a^2 - 32/7*a - 30/7)*t^4 + (60/7*a^5 - 30/7*a^4 - 18/7*a^3 - 9/7*a^2 - 3/7*a)*t^3 + (60/7*a^4 - 40/7*a^3 - 16/7*a^2 - 4/7*a)*t^2 + (30/7*a^3 - 25/7*a^2 - 5/7*a)*t + 6/7*a^2 - 6/7*a)
     2669
     2670            sage: pol = x^10 + (4/7*a - 6/7)*x^9 + (9/49*a^2 - 3/7*a + 15/49)*x^8 + (8/343*a^3 - 32/343*a^2 + 40/343*a - 20/343)*x^7 + (5/2401*a^4 - 20/2401*a^3 + 40/2401*a^2 - 5/343*a + 15/2401)*x^6 + (-6/16807*a^4 + 12/16807*a^3 - 18/16807*a^2 + 12/16807*a - 6/16807)*x^5
     2671            sage: pol.factor()
     2672            x^5 * (x^5 + (4/7*a - 6/7)*x^4 + (9/49*a^2 - 3/7*a + 15/49)*x^3 + (8/343*a^3 - 32/343*a^2 + 40/343*a - 20/343)*x^2 + (5/2401*a^4 - 20/2401*a^3 + 40/2401*a^2 - 5/343*a + 15/2401)*x - 6/16807*a^4 + 12/16807*a^3 - 18/16807*a^2 + 12/16807*a - 6/16807)
     2673
    26532674        """
    26542675        # PERFORMANCE NOTE:
    26552676        #     In many tests with SMALL degree PARI is substantially
     
    26652686        ## NTL uses the Berlekamp-Zassenhaus method with van Hoeij's improvements.
    26662687        ## But so does Magma since about Jul 2001.
    26672688        ##
    2668         ## But here's the kicker. Pari also uses this algorithm. Even Maple uses
     2689        ## But here's the kicker. PARI also uses this algorithm. Even Maple uses
    26692690        ## it!
    26702691        ##
    26712692        ## NTL's LLL algorithms are extremely well developed (van Hoeij uses
     
    27502771                return Factorization([(self._parent(p), e) for p, e in factors], R(factors.unit()))
    27512772       
    27522773            if R.defining_polynomial().denominator() == 1:
    2753 
    2754                 if (self.leading_coefficient() == 1):
    2755                     unit = None
    2756                     v = [ x._pari_("a") for x in self.list() ]
    2757                 else:
    2758                     unit = self.leading_coefficient()
    2759                     x = self.parent().gen()
    2760                     f1 = self(x/unit)*(unit**(self.degree()-1))
    2761                     ff1 = f1.factor()
    2762                     ff1 = Factorization([(f(unit*x).monic(),e) for f,e in ff1], unit=unit)
    2763                     return ff1
    2764 #                    temp_f = self * 1/unit
    2765 #                    v = [ x._pari_("a") for x in temp_f.list() ]
     2774                v = [ x._pari_("a") for x in self.list() ]
    27662775                f = pari(v).Polrev()
    27672776                Rpari = R.pari_nf()
    27682777                if (Rpari.variable() != "a"):
     
    27702779                    Rpari[0] = Rpari[0]("a")
    27712780                    Rpari[6] = [ x("a") for x in Rpari[6] ]
    27722781                G = list(Rpari.nffactor(f))
    2773                 return self._factor_pari_helper(G, unit=unit)
     2782                # PARI's nffactor() ignores the unit, _factor_pari_helper()
     2783                # adds back the unit of the factorization.
     2784                return self._factor_pari_helper(G)
    27742785
    27752786            else:
    27762787
     
    28702881
    28712882    def _factor_pari_helper(self, G, n=None, unit=None):
    28722883        """
    2873         Fix up and normalize a factorization that came from Pari.
     2884        Fix up and normalize a factorization that came from PARI.
    28742885       
    28752886        TESTS::
    28762887           
     
    29232934       
    29242935        else:
    29252936            # Otherwise we have to adjust for
    2926             # the content ignored by Pari.
     2937            # the content ignored by PARI.
    29272938            content_fix = R.base_ring()(1)
    29282939            for f, e in F:
    29292940                if not f.is_monic():
     
    44314442       
    44324443        If K and L are floating-point (RDF, CDF, RR, or CC), then a
    44334444        floating-point root-finder is used. If L is RDF or CDF then we
    4434         default to using NumPy's roots(); otherwise, we use Pari's
     4445        default to using NumPy's roots(); otherwise, we use PARI's
    44354446        polroots(). This choice can be overridden with
    44364447        algorithm='pari' or algorithm='numpy'. If the algorithm is
    44374448        unspecified and NumPy's roots() algorithm fails, then we fall
     
    44564467        If L is floating-point and K is not, then we attempt to change the
    44574468        polynomial ring to L (using .change_ring()) (or, if L is complex
    44584469        and K is not, to the corresponding real field). Then we use either
    4459         Pari or numpy as specified above.
     4470        PARI or numpy as specified above.
    44604471       
    44614472        For all other cases where K is different than L, we just use
    44624473        .change_ring(L) and proceed as below.
     
    45524563            if algorithm != 'numpy' and algorithm != 'either' and algorithm != 'pari':
    45534564                raise ValueError, "Unknown algorithm '%s'" % algorithm
    45544565
    4555             # We should support GSL, too.  We could also support Pari's
     4566            # We should support GSL, too.  We could also support PARI's
    45564567            # old Newton-iteration algorithm.
    45574568
    45584569            input_arbprec = (is_RealField(K) or