diff --git a/src/sage/libs/ntl/ntl_ZZ_pX.pxd b/src/sage/libs/ntl/ntl_ZZ_pX.pxd index 62b98f25e21..6003cad7958 100644 --- a/src/sage/libs/ntl/ntl_ZZ_pX.pxd +++ b/src/sage/libs/ntl/ntl_ZZ_pX.pxd @@ -1,6 +1,7 @@ # sage_setup: distribution = sagemath-ntl from sage.libs.ntl.ZZ_pX cimport * from sage.libs.ntl.ntl_ZZ_pContext cimport ntl_ZZ_pContext_class +from sage.rings.integer cimport Integer cdef class ntl_ZZ_pX(): cdef ZZ_pX_c x @@ -8,6 +9,8 @@ cdef class ntl_ZZ_pX(): cdef void setitem_from_int(ntl_ZZ_pX self, long i, int value) noexcept cdef int getitem_as_int(ntl_ZZ_pX self, long i) noexcept cdef ntl_ZZ_pX _new(self) + cdef ntl_ZZ_pX _pow(ntl_ZZ_pX self, long exp) + cdef ntl_ZZ_pX _powmod(ntl_ZZ_pX self, Integer exp, ntl_ZZ_pX modulus) cdef class ntl_ZZ_pX_Modulus(): cdef ZZ_pX_Modulus_c x diff --git a/src/sage/libs/ntl/ntl_ZZ_pX.pyx b/src/sage/libs/ntl/ntl_ZZ_pX.pyx index 5de5ae27c92..883d40f4489 100644 --- a/src/sage/libs/ntl/ntl_ZZ_pX.pyx +++ b/src/sage/libs/ntl/ntl_ZZ_pX.pyx @@ -30,6 +30,7 @@ include 'decl.pxi' from cpython.object cimport Py_EQ, Py_NE from sage.cpython.string cimport char_to_str from sage.rings.integer cimport Integer +from sage.libs.ntl.convert cimport PyLong_to_ZZ from sage.libs.ntl.ntl_ZZ cimport ntl_ZZ from sage.libs.ntl.ntl_ZZ_p cimport ntl_ZZ_p from sage.libs.ntl.ntl_ZZ_pContext cimport ntl_ZZ_pContext_class @@ -490,9 +491,12 @@ cdef class ntl_ZZ_pX(): sig_off() return r - def __pow__(ntl_ZZ_pX self, long n, ignored): + def __pow__(self, n, modulus): """ - Return the n-th nonnegative power of self. + Return the ``n``-th nonnegative power of ``self``. + + If ``modulus`` is not ``None``, the exponentiation is performed + modulo the polynomial ``modulus``. EXAMPLES:: @@ -500,13 +504,65 @@ cdef class ntl_ZZ_pX(): sage: g = ntl.ZZ_pX([-1,0,1],c) sage: g**10 [1 0 10 0 5 0 0 0 10 0 8 0 10 0 0 0 5 0 10 0 1] + + sage: x = ntl.ZZ_pX([0,1],c) + sage: x**10 + [0 0 0 0 0 0 0 0 0 0 1] + + Modular exponentiation:: + + sage: c = ntl.ZZ_pContext(20) + sage: f = ntl.ZZ_pX([1,0,1],c) + sage: m = ntl.ZZ_pX([1,0,0,0,0,1],c) + sage: pow(f, 123**45, m) + [1 19 3 0 3] + + Modular exponentiation of ``x``:: + + sage: f = ntl.ZZ_pX([0,1],c) + sage: f.is_x() + True + sage: m = ntl.ZZ_pX([1,1,0,0,0,1],c) + sage: pow(f, 123**45, m) + [15 5 5 11] + """ + if modulus is None: + return (self)._pow(n) + else: + return (self)._powmod(Integer(n), modulus) + + cdef ntl_ZZ_pX _pow(ntl_ZZ_pX self, long n): + """ + Compute the ``n``-th power of ``self``. """ if n < 0: raise NotImplementedError + cdef long ln = n #self.c.restore_c() # restored in _new() cdef ntl_ZZ_pX r = self._new() + if self.is_x() and n >= 1: + ZZ_pX_LeftShift(r.x, self.x, n-1) + else: + sig_on() + ZZ_pX_power(r.x, self.x, ln) + sig_off() + return r + + cdef ntl_ZZ_pX _powmod(ntl_ZZ_pX self, Integer n, ntl_ZZ_pX modulus): + r""" + Compute the ``n``-th power of ``self`` modulo a polynomial. + """ + cdef ntl_ZZ_pX r = self._new() + cdef ZZ_c n_ZZ + cdef ZZ_pX_Modulus_c mod + is_x = self.is_x() sig_on() - ZZ_pX_power(r.x, self.x, n) + mpz_to_ZZ(&n_ZZ, (n).value) + ZZ_pX_Modulus_build(mod, modulus.x) + if is_x: + ZZ_pX_PowerXMod_pre(r.x, n_ZZ, mod) + else: + ZZ_pX_PowerMod_pre(r.x, self.x, n_ZZ, mod) sig_off() return r diff --git a/src/sage/rings/polynomial/polynomial_modn_dense_ntl.pyx b/src/sage/rings/polynomial/polynomial_modn_dense_ntl.pyx index 454d70e3406..ae4fb8b9615 100644 --- a/src/sage/rings/polynomial/polynomial_modn_dense_ntl.pyx +++ b/src/sage/rings/polynomial/polynomial_modn_dense_ntl.pyx @@ -50,7 +50,7 @@ from sage.rings.infinity import infinity from sage.interfaces.singular import singular as singular_default -from sage.structure.element import coerce_binop +from sage.structure.element import canonical_coercion, coerce_binop, have_same_parent from sage.libs.ntl.types cimport NTL_SP_BOUND from sage.libs.ntl.ZZ_p cimport * @@ -214,7 +214,6 @@ cdef class Polynomial_dense_mod_n(Polynomial): def _pow(self, n): n = int(n) - if self.degree() <= 0: return self.parent()(self[0]**n) if n < 0: @@ -846,7 +845,7 @@ cdef class Polynomial_dense_modn_ntl_zz(Polynomial_dense_mod_n): return r def __pow__(Polynomial_dense_modn_ntl_zz self, ee, modulus): - """ + r""" TESTS:: sage: R. = PolynomialRing(Integers(100), implementation='NTL') @@ -1801,6 +1800,87 @@ cdef class Polynomial_dense_mod_p(Polynomial_dense_mod_n): A dense polynomial over the integers modulo `p`, where `p` is prime. """ + def __pow__(self, n, modulus): + """ + Exponentiation of ``self``. + + If ``modulus`` is not ``None``, the exponentiation is performed + modulo the polynomial ``modulus``. + + EXAMPLES:: + + sage: R. = PolynomialRing(Integers(101), implementation='NTL') + sage: (x-1)^5 + x^5 + 96*x^4 + 10*x^3 + 91*x^2 + 5*x + 100 + sage: pow(x-1, 15, x^3+x+1) + 55*x^2 + 6*x + 46 + + TESTS: + + Negative powers work but use the generic + implementation of fraction fields:: + + sage: R. = PolynomialRing(Integers(101), implementation='NTL') + sage: f = (x-1)^(-5) + sage: type(f) + + sage: (f + 2).numerator() + 2*x^5 + 91*x^4 + 20*x^3 + 81*x^2 + 10*x + 100 + + We define ``0^0`` to be unity, :issue:`13895`:: + + sage: R. = PolynomialRing(Integers(101), implementation='NTL') + sage: R(0)^0 + 1 + + The value returned from ``0^0`` should belong to our ring:: + + sage: R. = PolynomialRing(Integers(101), implementation='NTL') + sage: type(R(0)^0) == type(R(0)) + True + + The modulus can have smaller degree than ``self``:: + + sage: R. = PolynomialRing(GF(101), implementation="NTL") + sage: pow(x^4 + 1, 100, x^2 + x + 1) + 100*x + 100 + + Canonical coercion should apply:: + + sage: R. = PolynomialRing(GF(101), implementation="FLINT") + sage: x_ZZ = ZZ["x"].gen() + sage: pow(x+1, 100, 2) + 0 + sage: pow(x+1, 100, x_ZZ^2 + x_ZZ + 1) + 100*x + 100 + sage: pow(x+1, int(100), x_ZZ^2 + x_ZZ + 1) + 100*x + 100 + sage: xx = polygen(GF(97)) + sage: _ = pow(x + 1, 3, xx^3 + xx + 1) + Traceback (most recent call last): + ... + TypeError: no common canonical parent for objects with parents: ... + """ + n = Integer(n) + parent = (self)._parent + + if modulus is not None: + # Similar to coerce_binop + if not have_same_parent(self, modulus): + a, m = canonical_coercion(self, modulus) + if a is not self: + return pow(a, n, m) + modulus = m + self = self % modulus + return parent(pow(self.ntl_ZZ_pX(), n, modulus.ntl_ZZ_pX()), construct=True) + else: + if n < 0: + return ~(self**(-n)) + elif self.degree() <= 0: + return parent(self[0]**n) + else: + return parent(self.ntl_ZZ_pX()**n, construct=True) + @coerce_binop def gcd(self, right): """ diff --git a/src/sage/rings/polynomial/polynomial_zmod_flint.pxd b/src/sage/rings/polynomial/polynomial_zmod_flint.pxd index bfee3e8a037..df0db4e4c18 100644 --- a/src/sage/rings/polynomial/polynomial_zmod_flint.pxd +++ b/src/sage/rings/polynomial/polynomial_zmod_flint.pxd @@ -1,6 +1,7 @@ # sage_setup: distribution = sagemath-flint from sage.libs.flint.types cimport nmod_poly_t, nmod_poly_struct, fmpz_poly_t from sage.structure.parent cimport Parent +from sage.rings.integer cimport Integer from sage.rings.polynomial.polynomial_integer_dense_flint cimport Polynomial_integer_dense_flint ctypedef nmod_poly_struct celement @@ -15,4 +16,5 @@ cdef class Polynomial_zmod_flint(Polynomial_template): cdef int _set_list(self, x) except -1 cdef int _set_fmpz_poly(self, fmpz_poly_t) except -1 cpdef Polynomial _mul_trunc_opposite(self, Polynomial_zmod_flint other, length) + cdef Polynomial _powmod_bigexp(self, Integer exp, Polynomial_zmod_flint modulus) cpdef rational_reconstruction(self, m, n_deg=?, d_deg=?) diff --git a/src/sage/rings/polynomial/polynomial_zmod_flint.pyx b/src/sage/rings/polynomial/polynomial_zmod_flint.pyx index ad4ac79fa5b..374ba90942e 100644 --- a/src/sage/rings/polynomial/polynomial_zmod_flint.pyx +++ b/src/sage/rings/polynomial/polynomial_zmod_flint.pyx @@ -42,7 +42,7 @@ AUTHORS: from sage.libs.ntl.ntl_lzz_pX import ntl_zz_pX from sage.structure.element cimport parent -from sage.structure.element import coerce_binop +from sage.structure.element import coerce_binop, canonical_coercion, have_same_parent from sage.rings.polynomial.polynomial_integer_dense_flint cimport Polynomial_integer_dense_flint from sage.misc.superseded import deprecated_function_alias @@ -63,6 +63,7 @@ include "sage/libs/flint/nmod_poly_linkage.pxi" # and then the interface include "polynomial_template.pxi" +from sage.libs.flint.fmpz cimport * from sage.libs.flint.fmpz_poly cimport * from sage.libs.flint.nmod_poly cimport * @@ -483,6 +484,95 @@ cdef class Polynomial_zmod_flint(Polynomial_template): _mul_short_opposite = _mul_trunc_opposite + def __pow__(self, exp, modulus): + r""" + Exponentiation of ``self``. + + If ``modulus`` is not ``None``, the exponentiation is performed + modulo the polynomial ``modulus``. + + EXAMPLES:: + + sage: R. = GF(5)[] + sage: pow(x+1, 5**50, x^5 + 4*x + 3) + x + 1 + sage: pow(x+1, 5**64, x^5 + 4*x + 3) + x + 4 + sage: pow(x, 5**64, x^5 + 4*x + 3) + x + 3 + + The modulus can have smaller degree than ``self``:: + + sage: R. = PolynomialRing(GF(5), implementation="FLINT") + sage: pow(x^4, 6, x^2 + x + 1) + 1 + + TESTS: + + Canonical coercion applies:: + + sage: R. = PolynomialRing(GF(5), implementation="FLINT") + sage: x_ZZ = ZZ["x"].gen() + sage: pow(x+1, 25, 2) + 0 + sage: pow(x+1, 4, x_ZZ^2 + x_ZZ + 1) + 4*x + 4 + sage: pow(x+1, int(4), x_ZZ^2 + x_ZZ + 1) + 4*x + 4 + sage: xx = polygen(GF(97)) + sage: pow(x + 1, 3, xx^3 + xx + 1) + Traceback (most recent call last): + ... + TypeError: no common canonical parent for objects with parents: ... + """ + exp = Integer(exp) + if modulus is not None: + # Similar to coerce_binop + if not have_same_parent(self, modulus): + a, m = canonical_coercion(self, modulus) + if a is not self: + return pow(a, exp, m) + modulus = m + self = self % modulus + if exp > 0 and exp.bit_length() >= 32: + return (self)._powmod_bigexp(exp, modulus) + + return Polynomial_template.__pow__(self, exp, modulus) + + cdef Polynomial _powmod_bigexp(Polynomial_zmod_flint self, Integer exp, Polynomial_zmod_flint m): + r""" + Modular exponentiation with a large integer exponent. + + It is assumed that checks and coercions have been already performed on arguments. + + TESTS:: + + sage: R. = PolynomialRing(GF(5), implementation="FLINT") + sage: f = x+1 + sage: pow(f, 5**50, x^5 + 4*x + 3) # indirect doctest + x + 1 + sage: pow(x, 5**64, x^5 + 4*x + 3) # indirect doctest + x + 3 + """ + cdef Polynomial_zmod_flint ans = self._new() + # Preconditioning is useful for large exponents: the inverse + # power series helps computing fast quotients. + cdef Polynomial_zmod_flint minv = self._new() + cdef fmpz_t exp_fmpz + + fmpz_init(exp_fmpz) + fmpz_set_mpz(exp_fmpz, (exp).value) + nmod_poly_reverse(&minv.x, &m.x, nmod_poly_length(&m.x)) + nmod_poly_inv_series(&minv.x, &minv.x, nmod_poly_length(&m.x)) + + if self == self._parent.gen(): + nmod_poly_powmod_x_fmpz_preinv(&ans.x, exp_fmpz, &m.x, &minv.x) + else: + nmod_poly_powmod_fmpz_binexp_preinv(&ans.x, &self.x, exp_fmpz, &m.x, &minv.x) + + fmpz_clear(exp_fmpz) + return ans + cpdef Polynomial _power_trunc(self, unsigned long n, long prec): r""" TESTS:: diff --git a/src/sage/rings/polynomial/polynomial_zz_pex.pxd b/src/sage/rings/polynomial/polynomial_zz_pex.pxd index 3107ef76a39..8b0479efd1a 100644 --- a/src/sage/rings/polynomial/polynomial_zz_pex.pxd +++ b/src/sage/rings/polynomial/polynomial_zz_pex.pxd @@ -1,6 +1,7 @@ # sage_setup: distribution = sagemath-ntl from sage.libs.ntl.ZZ_pEX cimport ZZ_pEX_c from sage.libs.ntl.ntl_ZZ_pEContext cimport ZZ_pEContext_ptrs +from sage.rings.integer cimport Integer ctypedef ZZ_pEX_c celement ctypedef ZZ_pEContext_ptrs *cparent @@ -8,5 +9,4 @@ ctypedef ZZ_pEContext_ptrs *cparent include "polynomial_template_header.pxi" cdef class Polynomial_ZZ_pEX(Polynomial_template): - pass - + cdef _powmod_bigexp(Polynomial_ZZ_pEX self, Integer exp, Polynomial_ZZ_pEX modulus) diff --git a/src/sage/rings/polynomial/polynomial_zz_pex.pyx b/src/sage/rings/polynomial/polynomial_zz_pex.pyx index 006953f182f..48ec6745208 100644 --- a/src/sage/rings/polynomial/polynomial_zz_pex.pyx +++ b/src/sage/rings/polynomial/polynomial_zz_pex.pyx @@ -21,7 +21,9 @@ from sage.libs.ntl.ntl_ZZ_pEContext cimport ntl_ZZ_pEContext_class from sage.libs.ntl.ZZ_pE cimport ZZ_pE_to_ZZ_pX from sage.libs.ntl.ZZ_pX cimport ZZ_pX_deg, ZZ_pX_coeff from sage.libs.ntl.ZZ_p cimport ZZ_p_rep -from sage.libs.ntl.convert cimport ZZ_to_mpz +from sage.libs.ntl.convert cimport ZZ_to_mpz, mpz_to_ZZ + +from sage.structure.element import have_same_parent, canonical_coercion # We need to define this stuff before including the templating stuff # to make sure the function get_cparent is found since it is used in @@ -617,3 +619,97 @@ cdef class Polynomial_ZZ_pEX(Polynomial_template): ZZ_pEX_InvTrunc(r.x, self.x, prec) sig_off() return r + + def __pow__(self, exp, modulus): + r""" + Exponentiation of ``self``. + + If ``modulus`` is not ``None``, the exponentiation is performed + modulo the polynomial ``modulus``. + + EXAMPLES:: + + sage: K. = GF(101^2, 'a', modulus=[1,1,1]) + sage: R. = PolynomialRing(K, implementation="NTL") + sage: pow(x, 100) + x^100 + sage: pow(x + 3, 5) + x^5 + 15*x^4 + 90*x^3 + 68*x^2 + x + 41 + + If modulus is not ``None``, performs modular exponentiation:: + + sage: K. = GF(101^2, 'a', modulus=[1,1,1]) + sage: R. = PolynomialRing(K, implementation="NTL") + sage: pow(x, 100, x^2 + x + a) + (19*a + 64)*x + 30*a + 2 + sage: pow(x, 100 * 101**200, x^2 + x + a) + (19*a + 64)*x + 30*a + 2 + + The modulus can have smaller degree than ``self``:: + + sage: K. = GF(101^2, 'a', modulus=[1,1,1]) + sage: R. = PolynomialRing(K, implementation="NTL") + sage: pow(x^4, 25, x^2 + x + a) + (19*a + 64)*x + 30*a + 2 + + TESTS: + + Canonical coercion should apply:: + + sage: xx = GF(101)["x"].gen() + sage: pow(x+1, 25, 2) + 0 + sage: pow(x + a, 101**2, xx^3 + xx + 1) + 4*x^2 + 44*x + a + 70 + sage: pow(x + a, int(101**2), xx^3 + xx + 1) + 4*x^2 + 44*x + a + 70 + sage: xx = polygen(GF(97)) + sage: _ = pow(x + a, 101**2, xx^3 + xx + 1) + Traceback (most recent call last): + ... + TypeError: no common canonical parent for objects with parents: ... + """ + exp = Integer(exp) + if modulus is not None: + # Handle when modulus is zero + if modulus.is_zero(): + raise ZeroDivisionError("modulus must be nonzero") + + # Similar to coerce_binop + if not have_same_parent(self, modulus): + a, m = canonical_coercion(self, modulus) + if a is not self: + return pow(a, exp, m) + modulus = m + self = self % modulus + if exp > 0 and exp.bit_length() >= 32: + return (self)._powmod_bigexp(Integer(exp), modulus) + return Polynomial_template.__pow__(self, exp, modulus) + + cdef _powmod_bigexp(Polynomial_ZZ_pEX self, Integer exp, Polynomial_ZZ_pEX modulus): + """ + Modular exponentiation for large exponents. + """ + self._parent._modulus.restore() + cdef Polynomial_ZZ_pEX r + cdef ZZ_c e_ZZ + cdef ZZ_pEX_c y + cdef ZZ_pEX_Modulus_c mod + + mpz_to_ZZ(&e_ZZ, exp.value) + r = Polynomial_ZZ_pEX.__new__(Polynomial_ZZ_pEX) + celement_construct(&r.x, (self)._cparent) + r._parent = (self)._parent + r._cparent = (self)._cparent + ZZ_pEX_Modulus_build(mod, modulus.x) + + sig_on() + if ZZ_pEX_IsX(self.x): + ZZ_pEX_PowerXMod_ZZ_pre(r.x, e_ZZ, mod) + elif ZZ_pEX_deg(self.x) < ZZ_pEX_deg(modulus.x): + ZZ_pEX_PowerMod_ZZ_pre(r.x, self.x, e_ZZ, mod) + else: + ZZ_pEX_rem_pre(y, self.x, mod) + ZZ_pEX_PowerMod_ZZ_pre(r.x, y, e_ZZ, mod) + sig_off() + return r