Ticket #11890: 11890.patch

File 11890.patch, 13.6 KB (added by jdemeyer, 9 years ago)
  • sage/libs/pari/gen.pyx

    # HG changeset patch
    # User Jeroen Demeyer <jdemeyer@cage.ugent.be>
    # Date 1317743222 -7200
    # Node ID d6f5bb675950d1e56a1511e2de3c00c3ea72ecb3
    # Parent  d1131c9a3fa1eecabf31b0c93b04a61eedd6abbb
    Add "important" argument to pari_nf() to allow factoring of polynomials over big number fields.
    
    diff --git a/sage/libs/pari/gen.pyx b/sage/libs/pari/gen.pyx
    a b  
    453453        return self.g
    454454   
    455455    def list(self):
    456         if typ(self.g) == t_POL:
    457             raise NotImplementedError, \
    458                 "please report, t_POL.list() is broken, should not be used!"
    459         if typ(self.g) == t_SER:
    460             raise NotImplementedError, \
    461                 "please report, t_SER.list() is broken, should not be used!"
    462         #return list(self.Vecrev())
     456        """
     457        Convert self to a list of PARI gens.
     458
     459        EXAMPLES:
     460
     461        A PARI vector becomes a Sage list::
     462
     463            sage: L = pari("vector(10,i,i^2)").list()
     464            sage: L
     465            [1, 4, 9, 16, 25, 36, 49, 64, 81, 100]
     466            sage: type(L)
     467            <type 'list'>
     468            sage: type(L[0])
     469            <type 'sage.libs.pari.gen.gen'>
     470
     471        For polynomials, list() behaves as for ordinary Sage polynomials::
     472
     473            sage: pol = pari("x^3 + 5/3*x"); pol.list()
     474            [0, 5/3, 0, 1]
     475
     476        For power series, we get all coefficients, including leading and
     477        trailing zeros::
     478
     479            sage: R.<x> = LaurentSeriesRing(QQ)
     480            sage: s = x^2 + O(x^8)
     481            sage: s.list()
     482            [1]
     483            sage: pari(s).list()
     484            [0, 0, 1, 0, 0, 0, 0, 0]
     485
     486        For matrices, we get a list of columns::
     487
     488            sage: M = matrix(ZZ,3,2,[1,4,2,5,3,6]); M
     489            [1 4]
     490            [2 5]
     491            [3 6]
     492            sage: pari(M).list()
     493            [[1, 2, 3]~, [4, 5, 6]~]
     494
     495        For "scalar" types, we get a 1-element list containing ``self``::
     496
     497            sage: pari("42").list()
     498            [42]
     499        """
     500        if typ(self.g) == t_POL or typ(self.g) == t_SER:
     501            return list(self.Vecrev())
    463502        return list(self.Vec())
    464503       
    465504    def __reduce__(self):
     
    75647603    def __call__(self, x):
    75657604        return self.eval(x)
    75667605
     7606    def factornf(self, t):
     7607        """
     7608        Factorization of the polynomial ``self`` over the number field
     7609        defined by the polynomial ``t``.  This does not require that `t`
     7610        is integral, nor that the discriminant of the number field can be
     7611        factored.
     7612
     7613        EXAMPLES::
     7614
     7615            sage: x = polygen(QQ)
     7616            sage: K.<a> = NumberField(x^2 - 1/8)
     7617            sage: pari(x^2 - 2).factornf(K.pari_polynomial("a"))
     7618            [x + Mod(-4*a, 8*a^2 - 1), 1; x + Mod(4*a, 8*a^2 - 1), 1]
     7619        """
     7620        t0GEN(t)
     7621        sig_on()
     7622        return self.new_gen(polfnf(self.g, t0))
     7623
    75677624    def factorpadic(self, p, long r=20, long flag=0):
    75687625        """
    75697626        self.factorpadic(p,r=20,flag=0): p-adic factorization of the
  • sage/rings/number_field/number_field.py

    diff --git a/sage/rings/number_field/number_field.py b/sage/rings/number_field/number_field.py
    a b  
    26132613            self.__pari_polynomial_var = name
    26142614            return self.__pari_polynomial
    26152615
    2616     def pari_nf(self):
     2616    def pari_nf(self, important=True):
    26172617        """
    26182618        PARI number field corresponding to this field.
    26192619       
    26202620        This is the number field constructed using nfinit(). This is the same
    26212621        as the number field got by doing pari(self) or gp(self).
     2622
     2623        INPUT:
     2624
     2625        - ``important`` -- (default: True) bool.  If False, raise a
     2626          ``RuntimeError`` if we need to do a difficult discriminant
     2627          factorization.  Useful when the PARI nf structure is useful
     2628          but not strictly required, such as for factoring polynomials
     2629          over this number field.
    26222630       
    26232631        EXAMPLES::
    26242632       
     
    26462654            ...
    26472655            TypeError: Unable to coerce number field defined by non-integral polynomial to PARI.
    26482656
    2649         We illustrate the maximize_at_primes and assume_disc_small parameters::
    2650        
    2651             sage: p = next_prime(10^40); q = next_prime(p)
    2652 
    2653         The following would take a very long time without the
    2654         maximize_at_primes option above::
     2657        With ``important=False``, we simply bail out if we cannot
     2658        easily factor the discriminant::
     2659
     2660            sage: p = next_prime(10^40); q = next_prime(10^41)
     2661            sage: K.<a> = NumberField(x^2 - p*q)
     2662            sage: K.pari_nf(important=False)
     2663            Traceback (most recent call last):
     2664            ...
     2665            RuntimeError: Unable to factor discriminant with trial division
     2666
     2667        Next, we illustrate the ``maximize_at_primes`` and ``assume_disc_small``
     2668        parameters of the NumberField contructor. The following would take a
     2669        very long time without the ``maximize_at_primes`` option::
    26552670       
    26562671            sage: K.<a> = NumberField(x^2 - p*q, maximize_at_primes=[p])
    26572672            sage: K.pari_nf()
     
    26692684            return self.__pari_nf
    26702685        except AttributeError:
    26712686            f = self.pari_polynomial()
    2672             self.__pari_nf = pari([f, self._pari_integral_basis()]).nfinit()
     2687            self.__pari_nf = pari([f, self._pari_integral_basis(important=important)]).nfinit()
    26732688            return self.__pari_nf
    26742689
    26752690    def pari_zk(self):
     
    40104025        """
    40114026        return self.maximal_order(v=v).basis()
    40124027   
    4013     def _pari_integral_basis(self, v=None):
     4028    def _pari_integral_basis(self, v=None, important=True):
    40144029        """
    40154030        Internal function returning an integral basis of this number field in
    40164031        PARI format.
    40174032       
    40184033        INPUT:
    40194034       
    4020        
    4021         -  ``v`` - None, a prime, or a list of primes. See the
     4035        -  ``v`` -- None, a prime, or a list of primes. See the
    40224036           documentation for self.maximal_order.
    40234037       
     4038        - ``important`` -- (default:True) bool.  If False, raise a
     4039          ``RuntimeError`` if we need to do a difficult discriminant
     4040          factorization.  Useful when the PARI nf structure is useful
     4041          but not strictly required.
    40244042       
    40254043        EXAMPLES::
    40264044       
     
    40464064        v = self._normalize_prime_list(v)
    40474065        try:
    40484066            return self._integral_basis_dict[v]
    4049         except:
     4067        except (AttributeError, KeyError):
    40504068            f = self.pari_polynomial()
    4051             if len(v) == 0:
    4052                 B = f.nfbasis(1 if self._assume_disc_small else 0)
    4053             else:
     4069            if len(v) > 0:
    40544070                m = self._pari_disc_factorization_matrix(v)
    40554071                B = f.nfbasis(fa = m)
     4072            elif self._assume_disc_small:
     4073                B = f.nfbasis(1)
     4074            elif not important:
     4075                # Trial divide the discriminant
     4076                m = self.pari_polynomial().poldisc().abs().factor(limit=0)
     4077                # Since we only need a *squarefree* factorization, we need
     4078                # trial division up to D^(1/3) instead of D^(1/2).
     4079                trialdivlimit = pari(pari._primelimit()**3)
     4080                if all([ p < trialdivlimit or p.isprime() for p in m[0] ]):
     4081                    B = f.nfbasis(fa = m)
     4082                else:
     4083                    raise RuntimeError, "Unable to factor discriminant with trial division"
     4084            else:
     4085                B = f.nfbasis()
    40564086
    40574087            self._integral_basis_dict[v] = B
    40584088            return B
     
    66336663              To:   Number Field in a0 with defining polynomial x^2 - 1/2*c0 over its base field
    66346664              Defn: z |--> a0)
    66356665
    6636         We can relativize over a really large field.  This requires great care
    6637         to not factor or do any operation that would trigger a pari nfinit()
    6638         internally::
     6666        We can relativize over a really large field::
    66396667
    66406668            sage: K.<a> = CyclotomicField(3^3*2^3)
    66416669            sage: R = K.relativize(a^(3^2), 't'); R
     
    67106738        # the generator of the intermediate field.
    67116739
    67126740        # this strange step avoids computing an embedding of the base_field L
    6713         # into self, a computation which triggers a pari nfinit().  Without
    6714         # this, the relativize done over a huge field above is not feasible.
     6741        # into self.
    67156742        base_hom = L.hom([alpha], self)
    67166743        from_M = M.Hom(self)([self.gen()], base_hom=base_hom, check=True)
    67176744
     
    79547981            if p % n == 1:
    79557982                return p
    79567983   
    7957     def _pari_integral_basis(self, v=None):
     7984    def _pari_integral_basis(self, v=None, important=True):
    79587985        """
    79597986        Internal function returning an integral basis of this number field in
    79607987        PARI format.
    79617988       
    79627989        This field is cyclomotic, so this is a trivial computation, since
    7963         the power basis on the generator is an integral basis. Thus the v
    7964         parameter is ignored.
     7990        the power basis on the generator is an integral basis. Thus the ``v``
     7991        and ``important`` parameters are ignored.
    79657992       
    79667993        EXAMPLES::
    79677994       
  • sage/rings/polynomial/polynomial_element.pyx

    diff --git a/sage/rings/polynomial/polynomial_element.pyx b/sage/rings/polynomial/polynomial_element.pyx
    a b  
    26912691            sage: pol.factor()
    26922692            (t - 2*a^3 + a) * (t - 4/3*a^3 + 2/3*a) * (t - 2/3*a^3 + 1/3*a)
    26932693
     2694        Test that this factorization really uses ``nffactor()`` internally::
     2695
     2696            sage: pari.default("debug", 3)
     2697            sage: F = pol.factor()
     2698
     2699            Entering nffactor:
     2700            ...
     2701            sage: pari.default("debug", 0)
     2702
    26942703        Test that ticket #10369 is fixed::
    26952704
    26962705            sage: x = polygen(QQ)
     
    27092718            sage: pol.factor()
    27102719            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)
    27112720
     2721        Factoring over a number field over which we cannot factor the
     2722        discriminant::
     2723
     2724            sage: p = next_prime(10^50); q = next_prime(10^51); n = p*q;
     2725            sage: K.<a> = QuadraticField(p*q)
     2726            sage: R.<x> = PolynomialRing(K)
     2727            sage: factor(x^2 + 1)
     2728            x^2 + 1
     2729            sage: factor( (x - a) * (x + 2*a) )
     2730            (x - a) * (x + 2*a)
    27122731        """
    27132732        # PERFORMANCE NOTE:
    27142733        #     In many tests with SMALL degree PARI is substantially
     
    28042823        elif is_NumberField(R):
    28052824       
    28062825            if R.degree() == 1:
    2807                
    28082826                factors = self.change_ring(QQ).factor()
    28092827                return Factorization([(self._parent(p), e) for p, e in factors], R(factors.unit()))
    28102828       
    2811             if R.defining_polynomial().denominator() == 1:
    2812                 v = [ x._pari_("a") for x in self.list() ]
    2813                 f = pari(v).Polrev()
    2814                 Rpari = R.pari_nf()
    2815                 if (Rpari.variable() != "a"):
    2816                     Rpari = copy.copy(Rpari)
    2817                     Rpari[0] = Rpari[0]("a")
    2818                     Rpari[6] = [ x("a") for x in Rpari[6] ]
     2829            # Use "a" for the number field variable
     2830            v = [ c._pari_("a") for c in self.list() ]
     2831            f = pari(v).Polrev()
     2832            try:
     2833                # Try to compute the PARI nf structure with important=False.
     2834                # This will raise RuntimeError if the computation is too
     2835                # difficult.  It will raise TypeError if the defining
     2836                # polynomial is not integral.
     2837                Rpari = R.pari_nf(important=False).nf_subst("'a")
     2838                # Factor using nffactor()
    28192839                G = list(Rpari.nffactor(f))
    2820                 # PARI's nffactor() ignores the unit, _factor_pari_helper()
    2821                 # adds back the unit of the factorization.
    2822                 return self._factor_pari_helper(G)
    2823 
    2824             else:
    2825 
    2826                 Rdenom = R.defining_polynomial().denominator()
    2827 
    2828                 new_Rpoly = (R.defining_polynomial() * Rdenom).change_variable_name("a")
    2829 
    2830                 Rpari, Rdiff = new_Rpoly._pari_().nfinit(3)
    2831 
    2832                 AZ = QQ['z']
    2833                 Raux = NumberField(AZ(Rpari[0]),'alpha')
    2834 
    2835                 S, gSRaux, fRauxS = Raux.change_generator(Raux(Rdiff))
    2836 
    2837                 phi_RS = R.Hom(S)([S.gen(0)])
    2838                 phi_SR = S.Hom(R)([R.gen(0)])
    2839 
    2840                 unit = self.leading_coefficient()
    2841                 temp_f = self * 1/unit
    2842 
    2843                 v = [ gSRaux(phi_RS(x))._pari_("a") for x in temp_f.list() ]
    2844                 f = pari(v).Polrev()
    2845 
    2846                 pari_factors = Rpari.nffactor(f)
    2847 
    2848                 factors = [ ( self.parent([ phi_SR(fRauxS(Raux(pari_factors[0][i][j])))
    2849                                             for j in range(len(pari_factors[0][i])) ]) ,
    2850                              int(pari_factors[1][i]) )
    2851                             for i in range(pari_factors.nrows()) ]
    2852 
    2853                 return Factorization(factors, unit)
    2854        
     2840            except (RuntimeError, TypeError):
     2841                # Use factornf() which only needs the defining polynomial
     2842                # and which does not require an integral polynomial.
     2843                G = list(f.factornf(R.pari_polynomial("a")))
     2844            # PARI's nffactor() ignores the unit, _factor_pari_helper()
     2845            # adds back the unit of the factorization.
     2846            return self._factor_pari_helper(G)
    28552847
    28562848        elif is_RealField(R):
    28572849            n = pari.set_real_precision(int(3.5*R.prec()) + 1)