Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
r""" Dense univariate polynomials over `\ZZ`, implemented using NTL.
AUTHORS:
- David Harvey: split off from polynomial_element_generic.py (2007-09) - David Harvey: rewrote to talk to NTL directly, instead of via ntl.pyx (2007-09); a lot of this was based on Joel Mohler's recent rewrite of the NTL wrapper
Sage includes two implementations of dense univariate polynomials over `\ZZ`; this file contains the implementation based on NTL, but there is also an implementation based on FLINT in :mod:`sage.rings.polynomial.polynomial_integer_dense_flint`.
The FLINT implementation is preferred (FLINT's arithmetic operations are generally faster), so it is the default; to use the NTL implementation, you can do::
sage: K.<x> = PolynomialRing(ZZ, implementation='NTL') sage: K Univariate Polynomial Ring in x over Integer Ring (using NTL) """
#***************************************************************************** # Copyright (C) 2007 William Stein <wstein@gmail.com> # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 2 of the License, or # (at your option) any later version. # http://www.gnu.org/licenses/ #*****************************************************************************
from __future__ import absolute_import, print_function
from cysignals.memory cimport sig_free from cysignals.signals cimport sig_on, sig_off from sage.ext.cplusplus cimport ccrepr
include "sage/libs/ntl/decl.pxi"
from sage.rings.polynomial.polynomial_element cimport Polynomial from sage.structure.element cimport ModuleElement, Element
from sage.rings.integer_ring cimport IntegerRing_class
from sage.libs.ntl.ntl_ZZX cimport ntl_ZZX
from sage.rings.integer cimport Integer from sage.rings.real_mpfr cimport RealNumber, RealField_class from sage.rings.real_mpfi cimport RealIntervalFieldElement
from sage.libs.ntl.ZZX cimport *
from sage.rings.polynomial.evaluation cimport ZZX_evaluation_mpfr, ZZX_evaluation_mpfi
cdef class Polynomial_integer_dense_ntl(Polynomial): r""" A dense polynomial over the integers, implemented via NTL. """ cdef Polynomial_integer_dense_ntl _new(self): r""" Quickly creates a new initialized Polynomial_integer_dense_ntl with the correct parent and _is_gen == 0. """
def __init__(self, parent, x=None, check=True, is_gen=False, construct=False): r""" EXAMPLES::
sage: R.<x> = PolynomialRing(ZZ, implementation='NTL') sage: x x
Construct from list::
sage: R([]) 0 sage: R([1, -2, 3]) 3*x^2 - 2*x + 1
Coercions from other rings are attempted automatically::
sage: R([1, -6/3, 3]) 3*x^2 - 2*x + 1 sage: R([1, 5/2, 2]) Traceback (most recent call last): ... TypeError: no conversion of this rational to integer
Construct from constant::
sage: R(3) 3
Coercion from PARI polynomial::
sage: f = R([-1, 2, 5]); f 5*x^2 + 2*x - 1 sage: type(f) <type 'sage.rings.polynomial.polynomial_integer_dense_ntl.Polynomial_integer_dense_ntl'> sage: type(pari(f)) <type 'cypari2.gen.Gen'> sage: type(R(pari(f))) <type 'sage.rings.polynomial.polynomial_integer_dense_ntl.Polynomial_integer_dense_ntl'> sage: R(pari(f)) 5*x^2 + 2*x - 1
Coercion from NTL polynomial::
sage: f = ntl.ZZX([1, 2, 3]) sage: print(R(f)) 3*x^2 + 2*x + 1
Coercion from dictionary::
sage: f = R({2: -4, 3: 47}); f 47*x^3 - 4*x^2
Coercion from fraction field element with trivial denominator::
sage: f = (x^3 - 1) / (x - 1) sage: type(f) <type 'sage.rings.fraction_field_element.FractionFieldElement'> sage: g = R(f); g x^2 + x + 1
NTL polynomials are limited in size to slightly under the word length::
sage: PolynomialRing(ZZ, 'x', implementation='NTL')({2^3: 1}) x^8 sage: import sys sage: PolynomialRing(ZZ, 'x', implementation='NTL')({sys.maxsize>>1: 1}) Traceback (most recent call last): ... OverflowError: Dense NTL integer polynomials have a maximum degree of 268435455 # 32-bit OverflowError: Dense NTL integer polynomials have a maximum degree of 1152921504606846975 # 64-bit """
cdef Py_ssize_t degree cdef Py_ssize_t i cdef ZZ_c y
return # leave initialized to 0 polynomial.
# copy with NTL assignment operator self.__poly = (<Polynomial_integer_dense_ntl>x).__poly return else: # coerce coefficients into Sage integers
# find max degree to allocate only once raise ValueError("Negative monomial degrees not allowed: %s" % i) # now fill them in ZZX_SetCoeff_long(self.__poly, i, a) else: a = ZZ(a)
# copy with NTL assignment operator
isinstance(x.numerator(), Polynomial_integer_dense_ntl): if x.denominator() == 1: # fraction of the form f(x)/1 self.__poly = (<Polynomial_integer_dense_ntl>x.numerator()).__poly return
x = [x] # constant polynomials
raise OverflowError("Dense NTL integer polynomials have a maximum degree of %s" % (NTL_OVFBND-1))
else:
def content(self): r""" Return the greatest common divisor of the coefficients of this polynomial. The sign is the sign of the leading coefficient. The content of the zero polynomial is zero.
EXAMPLES::
sage: R.<x> = PolynomialRing(ZZ, implementation='NTL') sage: (2*x^2 - 4*x^4 + 14*x^7).content() 2 sage: (2*x^2 - 4*x^4 - 14*x^7).content() -2 sage: x.content() 1 sage: R(1).content() 1 sage: R(0).content() 0 """ cdef ZZ_c y
def _eval_mpfr_(self, RealNumber a): r""" Evaluate this polynomial on the real number element ``a``.
This method uses Horner's rule and might not be appropriate for polynomials of large degree.
TESTS::
sage: R.<x> = PolynomialRing(ZZ, implementation='NTL') sage: (x+1)._eval_mpfr_(RR(1.2)) 2.20000000000000 sage: (x^2)._eval_mpfr_(RR(2.2)) 4.84000000000000 sage: R.zero()._eval_mpfr_(RR(2.1)) 0.000000000000000 sage: R.one()._eval_mpfr_(RR(2.1)) 1.00000000000000
sage: p = x^3 - 2*x^2 + x -1 sage: p._eval_mpfr_(RR(1.3)) -0.883000000000000 """
def _eval_mpfi_(self, RealIntervalFieldElement a): r""" Evaluate this polynomial on the real interval ``a``.
This method uses Horner's rule and might not be appropriate for polynomials of large degree.
TESTS::
sage: R.<x> = PolynomialRing(ZZ, implementation='NTL') sage: (x+1)._eval_mpfi_(RIF(1.5)) 2.5000000000000000? sage: (x^2)._eval_mpfi_(RIF(1.333,1.334)) 1.78? sage: R.zero()._eval_mpfi_(RIF(2.1)) 0 sage: R.one()._eval_mpfi_(RIF(2.1)) 1
sage: p = x^3 - x^2 - x - 1 sage: r = p.roots(RIF, multiplicities=False)[0] sage: p._eval_mpfi_(r) 0.?e-27 """
def __reduce__(self): r""" Used for pickling.
EXAMPLES::
sage: R.<x> = PolynomialRing(ZZ, implementation='NTL') sage: loads(dumps(x)) == x True sage: f = 2*x + 3 sage: loads(dumps(f)) == f True """
cdef get_unsafe(self, Py_ssize_t n): """ Return the `n`-th coefficient of ``self``.
EXAMPLES::
sage: R.<x> = PolynomialRing(ZZ, implementation='NTL') sage: f = 2*x^2 - 3 sage: f[0] -3 sage: f[1] 0 sage: f[2] 2 sage: f[3] 0 sage: f[-1] 0 sage: f = 1 + x + 2*x^2 + 3*x^3 + 4*x^4 + 5*x^5 sage: f[:4] 3*x^3 + 2*x^2 + x + 1 sage: f[:100] 5*x^5 + 4*x^4 + 3*x^3 + 2*x^2 + x + 1 """
def _repr(self, name=None, bint latex=False): """ Return string representation of this polynomial.
EXAMPLES::
sage: R.<x> = PolynomialRing(ZZ, 'x', implementation='NTL') sage: (-x+1)^5 -x^5 + 5*x^4 - 10*x^3 + 10*x^2 - 5*x + 1 """ cdef long i else: all.append(" %s %s%s^{%s}" % (sign_str, coeff_str, name, i)) else: else: else:
def _latex_(self, name=None): """ Return the latex representation of this polynomial.
EXAMPLES::
sage: R.<t> = ZZ['t'] sage: latex(t^10-t^2-5*t+1) t^{10} - t^{2} - 5t + 1 sage: latex(cyclotomic_polynomial(10^5)) x^{40000} - x^{30000} + x^{20000} - x^{10000} + 1 """ if name is None: name = self.parent().latex_variable_names()[0] return self._repr(name, latex=True)
cpdef _add_(self, right): r""" Returns self plus right.
EXAMPLES::
sage: R.<x> = PolynomialRing(ZZ, implementation='NTL') sage: f = 2*x + 1 sage: g = -3*x^2 + 6 sage: f + g -3*x^2 + 2*x + 7 """ (<Polynomial_integer_dense_ntl>right).__poly)
cpdef _sub_(self, right): r""" Return self minus right.
EXAMPLES::
sage: R.<x> = PolynomialRing(ZZ, implementation='NTL') sage: f = 2*x + 1 sage: g = -3*x^2 + 6 sage: f - g 3*x^2 + 2*x - 5 """ (<Polynomial_integer_dense_ntl>right).__poly)
cpdef _neg_(self): r""" Returns negative of self.
EXAMPLES::
sage: R.<x> = PolynomialRing(ZZ, implementation='NTL') sage: f = 2*x - 1 sage: -f -2*x + 1 """
r""" Attempts to divide self by right, and return a quotient and remainder.
If right is monic, then it returns ``(q, r)`` where `self = q * right + r` and `deg(r) < deg(right)`.
If right is not monic, then it returns `(q, 0)` where q = self/right if right exactly divides self, otherwise it raises an exception.
EXAMPLES::
sage: R.<x> = PolynomialRing(ZZ, implementation='NTL') sage: f = R(range(10)); g = R([-1, 0, 1]) sage: q, r = f.quo_rem(g) sage: q, r (9*x^7 + 8*x^6 + 16*x^5 + 14*x^4 + 21*x^3 + 18*x^2 + 24*x + 20, 25*x + 20) sage: q*g + r == f True
sage: 0//(2*x) 0
sage: f = x^2 sage: f.quo_rem(0) Traceback (most recent call last): ... ArithmeticError: division by zero polynomial
sage: f = (x^2 + 3) * (2*x - 1) sage: f.quo_rem(2*x - 1) (x^2 + 3, 0)
sage: f = x^2 sage: f.quo_rem(2*x - 1) Traceback (most recent call last): ... ArithmeticError: division not exact in Z[x] (consider coercing to Q[x] first)
TESTS::
sage: z = R(0) sage: z.quo_rem(1) (0, 0) sage: z.quo_rem(x) (0, 0) sage: z.quo_rem(2*x) (0, 0)
"""
cdef ZZX_c *q cdef ZZX_c *r cdef int divisible
# divisor is monic. Just do the division and remainder else: # Non-monic divisor. Check whether it divides exactly. # exactly divisible else: # division failed: clean up and raise exception
r""" Return the GCD of self and right. The leading coefficient need not be 1.
EXAMPLES::
sage: R.<x> = PolynomialRing(ZZ, implementation='NTL') sage: f = (6*x + 47)*(7*x^2 - 2*x + 38) sage: g = (6*x + 47)*(3*x^3 + 2*x + 1) sage: f.gcd(g) 6*x + 47 """ # todo: we're doing an unnecessary copy here
""" Return the LCM of self and right.
EXAMPLES::
sage: R.<x> = PolynomialRing(ZZ, implementation='NTL') sage: f = (6*x + 47)*(7*x^2 - 2*x + 38) sage: g = (6*x + 47)*(3*x^3 + 2*x + 1) sage: h = f.lcm(g); h 126*x^6 + 951*x^5 + 486*x^4 + 6034*x^3 + 585*x^2 + 3706*x + 1786 sage: h == (6*x + 47)*(7*x^2 - 2*x + 38)*(3*x^3 + 2*x + 1) True """
""" This function can't in general return ``(g,s,t)`` as above, since they need not exist. Instead, over the integers, we first multiply `g` by a divisor of the resultant of `a/g` and `b/g`, up to sign, and return ``g, u, v`` such that ``g = s*self + s*right``. But note that this `g` may be a multiple of the gcd.
If ``self`` and ``right`` are coprime as polynomials over the rationals, then ``g`` is guaranteed to be the resultant of self and right, as a constant polynomial.
EXAMPLES::
sage: P.<x> = PolynomialRing(ZZ, implementation='NTL') sage: F = (x^2 + 2)*x^3; G = (x^2+2)*(x-3) sage: g, u, v = F.xgcd(G) sage: g, u, v (27*x^2 + 54, 1, -x^2 - 3*x - 9) sage: u*F + v*G 27*x^2 + 54 sage: x.xgcd(P(0)) (x, 1, 0) sage: f = P(0) sage: f.xgcd(x) (x, 0, 1) sage: F = (x-3)^3; G = (x-15)^2 sage: g, u, v = F.xgcd(G) sage: g, u, v (2985984, -432*x + 8208, 432*x^2 + 864*x + 14256) sage: u*F + v*G 2985984 """ cdef ZZX_c *s cdef ZZX_c *t cdef ZZ_c *r
else:
cpdef _mul_(self, right): r""" Returns self multiplied by right.
EXAMPLES::
sage: R.<x> = PolynomialRing(ZZ, implementation='NTL') sage: (x - 2)*(x^2 - 8*x + 16) x^3 - 10*x^2 + 32*x - 32 """ (<Polynomial_integer_dense_ntl>right).__poly)
cpdef _lmul_(self, Element right): r""" Returns self multiplied by right, where right is a scalar (integer).
EXAMPLES::
sage: R.<x> = PolynomialRing(ZZ, implementation='NTL') sage: x*3 3*x sage: (2*x^2 + 4)*3 6*x^2 + 12 """ cdef ZZ_c _right
cpdef _rmul_(self, Element right): r""" Returns self multiplied by right, where right is a scalar (integer).
EXAMPLES::
sage: R.<x> = PolynomialRing(ZZ, implementation='NTL') sage: 3*x 3*x sage: 3*(2*x^2 + 4) 6*x^2 + 12 """ cdef ZZ_c _right
def __floordiv__(self, right): """ EXAMPLES::
sage: R.<x> = PolynomialRing(ZZ, implementation='NTL') sage: f = R([9,6,1]) ; f x^2 + 6*x + 9 sage: f // x x + 6 sage: f // 3 2*x + 3 sage: g = x^3 ; g x^3 sage: f // g 0 sage: g // f x - 6 """ d = ZZ(right[0]) return self.parent()([c // d for c in self.list()], construct=True) else:
def _unsafe_mutate(self, long n, value): r""" Sets coefficient of `x^n` to value.
This is very unsafe, because Sage polynomials are supposed to be immutable. (Shhhh don't tell anyone!)
EXAMPLES::
sage: R.<x> = PolynomialRing(ZZ, implementation='NTL') sage: f = 2*x^2 + 3; f 2*x^2 + 3 sage: f._unsafe_mutate(1, 42); f 2*x^2 + 42*x + 3 """ raise IndexError("n must be >= 0") cdef ZZ_c y
def real_root_intervals(self): """ Returns isolating intervals for the real roots of this polynomial.
EXAMPLES: We compute the roots of the characteristic polynomial of some Salem numbers::
sage: R.<x> = PolynomialRing(ZZ, implementation='NTL') sage: f = 1 - x^2 - x^3 - x^4 + x^6 sage: f.real_root_intervals() [((1/2, 3/4), 1), ((1, 3/2), 1)] """
## def __copy__(self): ## f = Polynomial_integer_dense(self.parent()) ## f.__poly = self.__poly.copy() ## return f
def degree(self, gen=None): """ Return the degree of this polynomial. The zero polynomial has degree -1.
EXAMPLES::
sage: R.<x> = PolynomialRing(ZZ, implementation='NTL') sage: x.degree() 1 sage: (x^2).degree() 2 sage: R(1).degree() 0 sage: R(0).degree() -1 """
def discriminant(self, proof=True): r""" Return the discriminant of self, which is by definition
.. MATH::
(-1)^{m(m-1)/2} {\mbox{\tt resultant}}(a, a')/lc(a),
where `m = deg(a)`, and `lc(a)` is the leading coefficient of a. If ``proof`` is False (the default is True), then this function may use a randomized strategy that errors with probability no more than `2^{-80}`.
EXAMPLES::
sage: f = ntl.ZZX([1,2,0,3]) sage: f.discriminant() -339 sage: f.discriminant(proof=False) -339 """ cdef ZZ_c* temp = ZZX_discriminant(&self.__poly, proof) cdef Integer x = Integer.__new__(Integer) ZZ_to_mpz(x.value, temp) del temp return x
def __pari__(self, variable=None): """ EXAMPLES::
sage: t = PolynomialRing(ZZ,"t",implementation='NTL').gen() sage: f = t^3 + 3*t - 17 sage: pari(f) t^3 + 3*t - 17 """
def squarefree_decomposition(self): """ Return the square-free decomposition of self. This is a partial factorization of self into square-free, relatively prime polynomials.
This is a wrapper for the NTL function SquareFreeDecomp.
EXAMPLES::
sage: R.<x> = PolynomialRing(ZZ, implementation='NTL') sage: p = 37 * (x-1)^2 * (x-2)^2 * (x-3)^3 * (x-4) sage: p.squarefree_decomposition() (37) * (x - 4) * (x^2 - 3*x + 2)^2 * (x - 3)^3
TESTS:
Verify that :trac:`13053` has been resolved::
sage: R.<x> = PolynomialRing(ZZ, implementation='NTL') sage: f=-x^2 sage: f.squarefree_decomposition() (-1) * x^2
"""
cdef ZZX_c** v cdef long* e cdef long i, n cdef Polynomial_integer_dense_ntl z
def _factor_pari(self): """ Use pari to factor self.
EXAMPLES::
sage: R.<x> = PolynomialRing(ZZ, implementation='NTL') sage: f = R([9,6,1]) ; f x^2 + 6*x + 9 sage: f.factor() (x + 3)^2 sage: f._factor_pari() (x + 3)^2 """
def _factor_ntl(self): """ Use NTL to factor self.
AUTHOR:
- Joel B. Mohler
EXAMPLES::
sage: R.<x> = PolynomialRing(ZZ, implementation='NTL') sage: f = R([9,6,1]) sage: f._factor_ntl() (x + 3)^2 """ cdef Polynomial_integer_dense_ntl fac_py cdef ZZ_c content cdef vec_pair_ZZX_long_c factors cdef long i sig_on() sig_off() fac_py = self._new() ZZX_SetCoeff(fac_py.__poly, 0, content) if ZZX_deg(fac_py.__poly) == 0 and ZZ_to_int(fac_py.__poly.rep.elts())==-1: unit = fac_py else: results.append( (fac_py,1) )
def factor(self): """ This function overrides the generic polynomial factorization to make a somewhat intelligent decision to use Pari or NTL based on some benchmarking.
Note: This function factors the content of the polynomial, which can take very long if it's a really big integer. If you do not need the content factored, divide it out of your polynomial before calling this function.
EXAMPLES::
sage: R.<x>=ZZ[] sage: f=x^4-1 sage: f.factor() (x - 1) * (x + 1) * (x^2 + 1) sage: f=1-x sage: f.factor() (-1) * (x - 1) sage: f.factor().unit() -1 sage: f = -30*x; f.factor() (-1) * 2 * 3 * 5 * x """ cdef int i # it appears that pari has a window from about degrees 30 and 300 in which it beats NTL. else: return c.factor()*g._factor_pari()
def factor_mod(self, p): """ Return the factorization of self modulo the prime p.
INPUT:
- ``p`` -- prime
OUTPUT: factorization of self reduced modulo p.
EXAMPLES::
sage: R.<x> = PolynomialRing(ZZ, 'x', implementation='NTL') sage: f = -3*x*(x-2)*(x-9) + x sage: f.factor_mod(3) x sage: f = -3*x*(x-2)*(x-9) sage: f.factor_mod(3) Traceback (most recent call last): ... ArithmeticError: factorization of 0 is not defined
sage: f = 2*x*(x-2)*(x-9) sage: f.factor_mod(7) (2) * x * (x + 5)^2 """ raise ValueError("p must be prime")
def factor_padic(self, p, prec=10): """ Return `p`-adic factorization of self to given precision.
INPUT:
- ``p`` -- prime
- ``prec`` -- integer; the precision
OUTPUT:
- factorization of ``self`` over the completion at `p`.
EXAMPLES::
sage: R.<x> = PolynomialRing(ZZ, implementation='NTL') sage: f = x^2 + 1 sage: f.factor_padic(5, 4) ((1 + O(5^4))*x + (2 + 5 + 2*5^2 + 5^3 + O(5^4))) * ((1 + O(5^4))*x + (3 + 3*5 + 2*5^2 + 3*5^3 + O(5^4)))
A more difficult example::
sage: f = 100 * (5*x + 1)^2 * (x + 5)^2 sage: f.factor_padic(5, 10) (4 + O(5^10)) * ((5 + O(5^11)))^2 * ((1 + O(5^10))*x + (5 + O(5^10)))^2 * ((5 + O(5^10))*x + (1 + O(5^10)))^2
"""
# Parent field for coefficients and polynomial
# Factor the *exact* polynomial using factorpadic()
cpdef list list(self, bint copy=True): """ Return a new copy of the list of the underlying elements of ``self``.
EXAMPLES::
sage: x = PolynomialRing(ZZ,'x',implementation='NTL').0 sage: f = x^3 + 3*x - 17 sage: f.list() [-17, 3, 0, 1] sage: f = PolynomialRing(ZZ,'x',implementation='NTL')(0) sage: f.list() [] """
""" Returns the resultant of self and other, which must lie in the same polynomial ring.
If proof = False (the default is proof=True), then this function may use a randomized strategy that errors with probability no more than `2^{-80}`.
INPUT:
- other -- a polynomial
OUTPUT:
an element of the base ring of the polynomial ring
EXAMPLES::
sage: x = PolynomialRing(ZZ,'x',implementation='NTL').0 sage: f = x^3 + x + 1; g = x^3 - x - 1 sage: r = f.resultant(g); r -8 sage: r.parent() is ZZ True """ |