sagemathgh-37441: Bind to FLINT/NTL API for polynomial modular expone…

…ntiation (new version)

    
This is a rejuvenation of the pull request by @remyoudompheng
sagemath#35320.

All I have done is rebase with develop and modify a few lines to fit
with the current structuring of SageMath. I have also added a few
additional docstrings.

As an example of the performance gain:

### Flint polynomials over prime fields

Old timings:

About 5x faster for flint using `nmod_poly`

```py
sage: R.<x> = GF(5)[]
sage: %time pow(x+1, 5**100, x^5 + 4*x + 3)
CPU times: user 551 µs, sys: 1e+03 ns, total: 552 µs
Wall time: 555 µs
x + 1
```

New timings:

```py
sage: R.<x> = GF(5)[]
sage:  %time pow(x+1, 5**100, x^5 + 4*x + 3)
CPU times: user 175 µs, sys: 1e+03 ns, total: 176 µs
Wall time: 179 µs
x + 1
```

### NTL polynomials over prime fields

About 1000x faster for NTL $\textrm{GF}(p)$
Old timings:

```py
sage: set_random_seed(0)
....: p = 2**255 - 19
....: R.<x> = F[]
....: f = R.random_element(degree=12)
....: g = R.random_element(degree=12)
....: %timeit pow(f, 2**33, g)
92.3 ms ± 947 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
sage: %timeit pow(f, 2**64, g)
1.29 s ± 35.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
```

New timings:

```py
sage: set_random_seed(0)
....: p = 2**255 - 19
....: R.<x> = F[]
....: f = R.random_element(degree=12)
....: g = R.random_element(degree=12)
....: %timeit pow(f, 2**33, g)
651 µs ± 28.9 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
sage:  %timeit pow(f, 2**64, g)
1.17 ms ± 10.7 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops
each)
```

### NTL polynomials over extension fields

About 10x faster for NTL polynomial rings over $\textrm{GF}(p^k)$
Old timings:

```py
sage: set_random_seed(0)
....:
....: p = 2**255 - 19
....: F.<a> = GF(p^3)
....: q = F.order()
....: R.<x> = F[]
....:
....: f = R.random_element(degree=12)
....: g = R.random_element(degree=12)
....:
....: %time _ = pow(f, q, g)
....:
....: p = 29
....: F.<a> = GF(p^100)
....: q = F.order()
....: R.<x> = F[]
....:
....: f = R.random_element(degree=12)
....: g = R.random_element(degree=12)
....:
....: %time _ = pow(f, q, g)
CPU times: user 409 ms, sys: 1.45 ms, total: 410 ms
Wall time: 411 ms
CPU times: user 14.6 s, sys: 202 ms, total: 14.8 s
Wall time: 14.9 s
```

New timings:

```py
sage: set_random_seed(0)
....:
....: p = 2**255 - 19
....: F.<a> = GF(p^3)
....: q = F.order()
....: R.<x> = F[]
....:
....: f = R.random_element(degree=12)
....: g = R.random_element(degree=12)
....:
....: %time _ = pow(f, q, g)
....:
....: p = 29
....: F.<a> = GF(p^100)
....: q = F.order()
....: R.<x> = F[]
....:
....: f = R.random_element(degree=12)
....: g = R.random_element(degree=12)
....:
....: %time _ = pow(f, q, g)
CPU times: user 371 ms, sys: 761 µs, total: 371 ms
Wall time: 371 ms
CPU times: user 1.54 s, sys: 10.5 ms, total: 1.55 s
Wall time: 1.55 s
```

Closes sagemath#35320
    
URL: sagemath#37441
Reported by: Giacomo Pope
Reviewer(s): Giacomo Pope, Travis Scrimshaw

src/sage/libs/ntl/ntl_ZZ_pX.pxd

@@ -1,12 +1,15 @@
 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
     cdef ntl_ZZ_pContext_class c
     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

src/sage/libs/ntl/ntl_ZZ_pX.pyx

@@ -29,6 +29,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
@@ -489,23 +490,78 @@ 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::
 
             sage: c = ntl.ZZ_pContext(20)
             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 (<ntl_ZZ_pX>self)._pow(n)
+        else:
+            return (<ntl_ZZ_pX>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, (<Integer>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
 

src/sage/rings/polynomial/polynomial_modn_dense_ntl.pyx

@@ -49,7 +49,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 *
@@ -213,7 +213,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:
@@ -845,7 +844,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.<x> = PolynomialRing(Integers(100), implementation='NTL')
@@ -1800,6 +1799,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.<x> = 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.<x> = PolynomialRing(Integers(101), implementation='NTL')
+            sage: f = (x-1)^(-5)
+            sage: type(f)
+            <class 'sage.rings.fraction_field_element.FractionFieldElement_1poly_field'>
+            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.<x> = PolynomialRing(Integers(101), implementation='NTL')
+            sage: R(0)^0
+            1
+
+        The value returned from ``0^0`` should belong to our ring::
+
+            sage: R.<x> = PolynomialRing(Integers(101), implementation='NTL')
+            sage: type(R(0)^0) == type(R(0))
+            True
+
+        The modulus can have smaller degree than ``self``::
+
+            sage: R.<x> = PolynomialRing(GF(101), implementation="NTL")
+            sage: pow(x^4 + 1, 100, x^2 + x + 1)
+            100*x + 100
+
+        Canonical coercion should apply::
+
+            sage: R.<x> = 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 = (<Element>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):
         """

src/sage/rings/polynomial/polynomial_zmod_flint.pxd

@@ -1,5 +1,6 @@
 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
@@ -14,4 +15,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=?)

src/sage/rings/polynomial/polynomial_zmod_flint.pyx

@@ -41,7 +41,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
@@ -62,6 +62,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 *
 
@@ -482,6 +483,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.<x> = 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.<x> = PolynomialRing(GF(5), implementation="FLINT")
+            sage: pow(x^4, 6, x^2 + x + 1)
+            1
+
+        TESTS:
+
+        Canonical coercion applies::
+
+            sage: R.<x> = 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 (<Polynomial_zmod_flint>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.<x> = 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, (<Integer>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::

src/sage/rings/polynomial/polynomial_zz_pex.pxd

@@ -1,11 +1,11 @@
 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
 
 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)