gh-38476: Include modular composition for polynomial rings over finit…

…e fields

    
Modular composition was very slow for polynomial rings over finite
fields, so I have exposed and used the functions from NTL / Flint.

Compared to the naive method we see very nice speed ups.

```py
sage: R.<x> = GF(2**127 - 1)[]
sage: f = R.random_element(degree=100)
sage: g = R.random_element(degree=100)
sage: h = R.random_element(degree=50)
sage: %timeit f(g) % h
960 ms ± 10.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
sage: %timeit f.modular_composition(g, h)
721 µs ± 10.4 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
sage: f(g) % h == f.modular_composition(g, h)
True
sage:
sage: R.<x> = GF((2**127 - 1, 2))[]
sage: f = R.random_element(degree=100)
sage: g = R.random_element(degree=100)
sage: h = R.random_element(degree=50)
sage: %timeit f(g) % h
1.03 s ± 6.74 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
sage: %timeit f.modular_composition(g, h)
29.1 ms ± 346 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
sage: f(g) % h == f.modular_composition(g, h)
True
sage:
sage: R.<x> = GF(163)[]
sage: f = R.random_element(degree=100)
sage: g = R.random_element(degree=100)
sage: h = R.random_element(degree=50)
sage: %timeit f(g) % h
1.68 ms ± 39.9 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops
each)
sage: %timeit f.modular_composition(g, h)
283 µs ± 5.27 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
sage: f(g) % h == f.modular_composition(g, h)
True
```

Both `./sage -t` and `./sage -tox ...` are erroring on my new build
after 10.4 so I'm going to let the CI catch my typos.
    
URL: #38476
Reported by: Giacomo Pope
Reviewer(s): Giacomo Pope, Lorenz Panny

src/sage/libs/ntl/ZZ_pEX.pxd

@@ -82,6 +82,8 @@ cdef extern from "ntlwrap.h":
     void ZZ_pEX_DivRem_pre "DivRem"(ZZ_pEX_c q, ZZ_pEX_c r, ZZ_pEX_c a, ZZ_pEX_Modulus_c F)
     void ZZ_pEX_div_pre "div"(ZZ_pEX_c q, ZZ_pEX_c a, ZZ_pEX_Modulus_c F)
 
+    void ZZ_pEX_CompMod "CompMod"(ZZ_pEX_c x, ZZ_pEX_c f, ZZ_pEX_c g, ZZ_pEX_Modulus_c F)
+
     void ZZ_pEX_MinPolyMod "MinPolyMod"(ZZ_pEX_c h, ZZ_pEX_c g, ZZ_pEX_c f)
     void ZZ_pEX_MinPolyMod_pre "MinPolyMod"(ZZ_pEX_c h, ZZ_pEX_c g, ZZ_pEX_Modulus_c F)
 

src/sage/libs/ntl/ntl_ZZ_pX.pyx

@@ -1472,6 +1472,46 @@ cdef class ntl_ZZ_pX():
         #ZZ_pX_preallocate_space(&self.x, n)
         sig_off()
 
+    def compose_mod(self, ntl_ZZ_pX other, ntl_ZZ_pX modulus):
+        r"""
+        Compute `f(g) \bmod h`.
+
+        To be precise about the order fo compostion, given ``self``, ``other``
+        and ``modulus`` as `f(x)`, `g(x)` and `h(x)` compute `f(g(x)) \bmod h(x)`.
+
+        INPUT:
+
+        - ``other`` -- a polynomial `g(x)`
+        - ``modulus`` -- a polynomial `h(x)`
+
+        EXAMPLES::
+
+            sage: c = ntl.ZZ_pContext(17)
+            sage: f = ntl.ZZ_pX([1,2,3],c)
+            sage: g = ntl.ZZ_pX([3,2,1],c)
+            sage: h = ntl.ZZ_pX([1,0,1],c)
+            sage: f.compose_mod(g, h)
+            [5 11]
+
+        AUTHORS:
+
+        - Giacomo Pope (2024-08) initial implementation
+        """
+        cdef ntl_ZZ_pX r = self._new()
+        cdef ZZ_pX_Modulus_c mod
+
+        sig_on()
+        # NTL requires that g is reduced
+        other = other % modulus
+
+        # Build the modulus type
+        ZZ_pX_Modulus_build(mod, modulus.x)
+
+        # Compute f(g) % h
+        ZZ_pX_CompMod(r.x, self.x, other.x, mod)
+        sig_off()
+        return r
+
 cdef class ntl_ZZ_pX_Modulus():
     """
     Thin holder for ZZ_pX_Moduli.

src/sage/rings/polynomial/polynomial_gf2x.pyx

@@ -256,6 +256,12 @@ cdef class Polynomial_GF2X(Polynomial_template):
         verbose("Res %5.3f s"%cputime(t),level=1)
         return res
 
+    # Other polynomials have compose_mod as methods following the naming of
+    # NTL/Flint bindings but the above method predates these. We expose
+    # compose_mod here so all polynomial ring elements which support this can
+    # use either name
+    compose_mod = modular_composition
+
     @cached_method
     def is_irreducible(self):
         r"""

src/sage/rings/polynomial/polynomial_modn_dense_ntl.pyx

@@ -424,6 +424,49 @@ cdef class Polynomial_dense_mod_n(Polynomial):
         """
         return small_roots(self, *args, **kwds)
 
+    def compose_mod(self, other, modulus):
+        r"""
+        Compute `f(g) \pmod h`.
+
+        To be precise about the order fo compostion, given ``self``, ``other``
+        and ``modulus`` as `f(x)`, `g(x)` and `h(x)` compute `f(g(x)) \bmod h(x)`.
+
+        INPUT:
+
+        - ``other`` -- a polynomial `g(x)`
+        - ``modulus`` -- a polynomial `h(x)`
+
+        EXAMPLES::
+
+            sage: R.<x> = GF(2**127 - 1)[]
+            sage: f = R.random_element()
+            sage: g = R.random_element()
+            sage: g.compose_mod(g, f) == g(g) % f
+            True
+
+            sage: R.<x> = GF(163)[]
+            sage: f = R([i for i in range(100)])
+            sage: g = R([i**2 for i in range(100)])
+            sage: h = 1 + x + x**5
+            sage: f.compose_mod(g, h)
+            82*x^4 + 56*x^3 + 45*x^2 + 60*x + 127
+            sage: f.compose_mod(g, h) == f(g) % h
+            True
+
+        AUTHORS:
+
+        - Giacomo Pope (2024-08) initial implementation
+        """
+        elt = self.ntl_ZZ_pX()
+        mod = modulus.ntl_ZZ_pX()
+        other = other.ntl_ZZ_pX()
+        res = elt.compose_mod(other, mod)
+        return self.parent()(res, construct=True)
+
+    # compose_mod is the natural name from the NTL bindings, but polynomial_gf2x
+    # has modular_composition as the method name so here we allow both
+    modular_composition = compose_mod
+
 
 def small_roots(self, X=None, beta=1.0, epsilon=None, **kwds):
     r"""

src/sage/rings/polynomial/polynomial_zmod_flint.pyx

@@ -989,3 +989,48 @@ cdef class Polynomial_zmod_flint(Polynomial_template):
         from sage.rings.polynomial.polynomial_ring_constructor import _single_variate
         R = _single_variate(parent.base_ring(), name=name, implementation='NTL')
         return parent(R(self % other).minpoly_mod(R(other)))
+
+    def compose_mod(self, other, modulus):
+        r"""
+        Compute `f(g) \bmod h`.
+
+        To be precise about the order fo compostion, given ``self``, ``other``
+        and ``modulus`` as `f(x)`, `g(x)` and `h(x)` compute `f(g(x)) \bmod h(x)`.
+
+        INPUT:
+
+        - ``other`` -- a polynomial `g(x)`
+        - ``modulus`` -- a polynomial `h(x)`
+
+        EXAMPLES::
+
+            sage: R.<x> = GF(163)[]
+            sage: f = R.random_element()
+            sage: g = R.random_element()
+            sage: g.compose_mod(g, f) == g(g) % f
+            True
+
+            sage: f = R([i for i in range(100)])
+            sage: g = R([i**2 for i in range(100)])
+            sage: h = 1 + x + x**5
+            sage: f.compose_mod(g, h)
+            82*x^4 + 56*x^3 + 45*x^2 + 60*x + 127
+            sage: f.compose_mod(g, h) == f(g) % h
+            True
+
+        AUTHORS:
+
+        - Giacomo Pope (2024-08) initial implementation
+        """
+        cdef Polynomial_zmod_flint res = self._new()
+
+        sig_on()
+        nmod_poly_compose_mod(&res.x, &(<Polynomial_zmod_flint>self).x, &(<Polynomial_zmod_flint>other).x, &(<Polynomial_zmod_flint>modulus).x)
+        sig_off()
+
+        return res
+
+    # compose_mod is the natural name from the Flint bindings, but
+    # polynomial_gf2x has modular_composition as the method name so here we
+    # allow both
+    modular_composition = compose_mod

src/sage/rings/polynomial/polynomial_zz_pex.pyx

@@ -711,3 +711,67 @@ cdef class Polynomial_ZZ_pEX(Polynomial_template):
             ZZ_pEX_PowerMod_ZZ_pre(r.x, y, e_ZZ, mod)
         sig_off()
         return r
+
+    def compose_mod(self, other, modulus):
+        r"""
+        Compute `f(g) \bmod h`.
+
+        To be precise about the order fo compostion, given ``self``, ``other``
+        and ``modulus`` as `f(x)`, `g(x)` and `h(x)` compute `f(g(x)) \bmod h(x)`.
+
+        INPUT:
+
+        - ``other`` -- a polynomial `g(x)`
+        - ``modulus`` -- a polynomial `h(x)`
+
+        EXAMPLES::
+
+            sage: R.<x> = GF(3**6)[]
+            sage: f = R.random_element()
+            sage: g = R.random_element()
+            sage: g.compose_mod(g, f) == g(g) % f
+            True
+
+            sage: F.<z3> = GF(3**6)
+            sage: R.<x> = F[]
+            sage: f = 2*z3^2*x^2 + (z3 + 1)*x + z3^2 + 2
+            sage: g = (z3^2 + 2*z3)*x^2 + (2*z3 + 2)*x + 2*z3^2 + z3 + 2
+            sage: h = (2*z3 + 2)*x^2 + (2*z3^2 + 1)*x + 2*z3^2 + z3 + 2
+            sage: f.compose_mod(g, h)
+            (z3^5 + z3^4 + z3^3 + z3^2 + z3)*x + z3^5 + z3^3 + 2*z3 + 2
+            sage: f.compose_mod(g, h) == f(g) % h
+            True
+
+        AUTHORS:
+
+        - Giacomo Pope (2024-08) initial implementation
+        """
+        self._parent._modulus.restore()
+
+        # Ensure all the parents match
+        if other.parent() is not self._parent:
+            other = self._parent.coerce(other)
+        if modulus.parent() is not self._parent:
+            modulus = self._parent.coerce(modulus)
+
+        # Create the output polynomial
+        cdef Polynomial_ZZ_pEX r
+        r = Polynomial_ZZ_pEX.__new__(Polynomial_ZZ_pEX)
+        celement_construct(&r.x, (<Polynomial_template>self)._cparent)
+        r._parent = (<Polynomial_template>self)._parent
+        r._cparent = (<Polynomial_template>self)._cparent
+
+        # Create ZZ_pEX_Modulus type from modulus input
+        cdef ZZ_pEX_Modulus_c mod
+        ZZ_pEX_Modulus_build(mod, (<Polynomial_ZZ_pEX>modulus).x)
+
+        # Compute f(g) mod h
+        sig_on()
+        ZZ_pEX_CompMod(r.x, (<Polynomial_ZZ_pEX>self).x, (<Polynomial_ZZ_pEX>(other % modulus)).x, mod)
+        sig_off()
+
+        return r
+
+    # compose_mod is the natural name from the NTL bindings, but polynomial_gf2x
+    # has modular_composition as the method name so here we allow both
+    modular_composition = compose_mod