Middle products for fmpz_mod_poly, mpn_mod_poly

Author: fredrik-johansson

doc/source/fmpz_mod_poly.rst

@@ -495,6 +495,20 @@ Multiplication
     Sets ``res`` to the lowest `n` coefficients of the product of
     ``poly1`` and ``poly2``.
 
+.. function:: void _fmpz_mod_poly_mulmid(fmpz * res, const fmpz * poly1, slong len1, const fmpz * poly2, slong len2, slong nlo, slong nhi, const fmpz_mod_ctx_t ctx)
+
+    Sets ``(res, nhi - nlo)`` to the coefficients at indices `[nlo, nhi)`
+    in the full product of ``(poly1, len1)`` and ``(poly2, len2)``.
+    Assumes that ``len1`` and ``len2`` are positive and that
+    `0 \le nlo < nhi \le len1 + len2 - 1`.
+    Does not support aliasing between the inputs and the output.
+
+.. function:: void fmpz_mod_poly_mulmid(fmpz_mod_poly_t res, const fmpz_mod_poly_t poly1, const fmpz_mod_poly_t poly2, slong nlo, slong nhi, const fmpz_mod_ctx_t ctx)
+
+    Sets ``res`` to the polynomial formed by the coefficients at indices `[nlo, nhi)`
+    in the product of ``poly1`` and ``poly2``. Equivalently, compute
+    `[(poly1 \cdot poly2) \bmod x^{nhi}] / x^{nlo}`.
+
 .. function:: void _fmpz_mod_poly_sqr(fmpz * res, const fmpz * poly, slong len, const fmpz_mod_ctx_t ctx)
 
     Sets ``res`` to the square of ``poly``.

doc/source/mpn_mod.rst

@@ -227,34 +227,40 @@ Multiplication
 All multiplication algorithms optimize for squaring.
 
 .. function:: int _mpn_mod_poly_mullow_classical(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong len, gr_ctx_t ctx)
+              int _mpn_mod_poly_mulmid_classical(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong nlo, slong nhi, gr_ctx_t ctx)
 
     Polynomial multiplication using the schoolbook algorithm.
 
 .. function:: int _mpn_mod_poly_mullow_KS(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong len, gr_ctx_t ctx)
+              int _mpn_mod_poly_mulmid_KS(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong nlo, slong nhi, gr_ctx_t ctx)
 
     Polynomial multiplication using Kronecker substitution (bit packing).
 
 .. function:: int _mpn_mod_poly_mullow_karatsuba(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong len, slong cutoff, gr_ctx_t ctx)
+              int _mpn_mod_poly_mulmid_karatsuba(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong nlo, slong nhi, slong cutoff, gr_ctx_t ctx)
 
     Polynomial multiplication using the Karatsuba algorithm,
     implemented without intermediate modular reductions.
     This algorithm calls itself recursively, switching to
     basecase multiplication (also without intermediate reductions)
     when either *len1* or *len2* is smaller than *cutoff*.
 
-    Currently a full product is computed internally regardless of *len*;
-    truncation only skips the modular reductions.
+    Currently a full product is computed internally regardless of the
+    output length; truncation only skips the modular reductions.
 
 .. function:: int _mpn_mod_poly_mullow_fft_small(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong len, gr_ctx_t ctx)
+              int _mpn_mod_poly_mulmid_fft_small(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong nlo, slong nhi, gr_ctx_t ctx)
 
     Polynomial multiplication using the small-prime FFT.
     Returns ``GR_UNABLE`` if the small-prime FFT is not available
     or if the coefficients are too large to use this implementation.
 
 .. function:: int _mpn_mod_poly_mullow(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong len, gr_ctx_t ctx)
+              int _mpn_mod_poly_mulmid(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong nlo, slong nhi, gr_ctx_t ctx)
 
     Polynomial multiplication with automatic algorithm selection.
 
+
 Division
 ..............
 

src/fmpz_mod_poly.h

@@ -310,6 +310,9 @@ void fmpz_mod_poly_mul(fmpz_mod_poly_t res, const fmpz_mod_poly_t poly1, const f
 void _fmpz_mod_poly_mullow(fmpz *res, const fmpz *poly1, slong len1, const fmpz *poly2, slong len2, slong n, const fmpz_mod_ctx_t ctx);
 void fmpz_mod_poly_mullow(fmpz_mod_poly_t res, const fmpz_mod_poly_t poly1, const fmpz_mod_poly_t poly2, slong n, const fmpz_mod_ctx_t ctx);
 
+void _fmpz_mod_poly_mulmid(fmpz *res, const fmpz *poly1, slong len1, const fmpz *poly2, slong len2, slong nlo, slong nhi, const fmpz_mod_ctx_t ctx);
+void fmpz_mod_poly_mulmid(fmpz_mod_poly_t res, const fmpz_mod_poly_t poly1, const fmpz_mod_poly_t poly2, slong nlo, slong nhi, const fmpz_mod_ctx_t ctx);
+
 void _fmpz_mod_poly_sqr(fmpz *res, const fmpz *poly, slong len, const fmpz_mod_ctx_t ctx);
 
 void fmpz_mod_poly_mulhigh(fmpz_mod_poly_t res, const fmpz_mod_poly_t poly1, const fmpz_mod_poly_t poly2, slong start, const fmpz_mod_ctx_t ctx);

src/fmpz_mod_poly/mulmid.c

@@ -0,0 +1,64 @@
+/*
+    Copyright (C) 2011 Sebastian Pancratz
+    Copyright (C) 2010 William Hart
+    Copyright (C) 2026 Fredrik Johansson
+
+    This file is part of FLINT.
+
+    FLINT is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 3 of the License, or
+    (at your option) any later version.  See <https://www.gnu.org/licenses/>.
+*/
+
+#include "fmpz_vec.h"
+#include "fmpz_poly.h"
+#include "fmpz_mod.h"
+#include "fmpz_mod_vec.h"
+#include "fmpz_mod_poly.h"
+
+void _fmpz_mod_poly_mulmid(fmpz *res, const fmpz *poly1, slong len1,
+                                      const fmpz *poly2, slong len2,
+                                      slong nlo, slong nhi, const fmpz_mod_ctx_t ctx)
+{
+    _fmpz_poly_mulmid(res, poly1, len1, poly2, len2, nlo, nhi);
+    _fmpz_mod_vec_set_fmpz_vec(res, res, nhi - nlo, ctx);
+}
+
+void fmpz_mod_poly_mulmid(fmpz_mod_poly_t res, const fmpz_mod_poly_t poly1,
+                const fmpz_mod_poly_t poly2, slong nlo, slong nhi, const fmpz_mod_ctx_t ctx)
+{
+    slong len1 = poly1->length;
+    slong len2 = poly2->length;
+    slong len;
+
+    FLINT_ASSERT(nlo >= 0);
+    FLINT_ASSERT(nhi >= 0);
+
+    if (len1 == 0 || len2 == 0 || nlo >= FLINT_MIN(nhi, len1 + len2 - 1))
+    {
+        fmpz_mod_poly_zero(res, ctx);
+        return;
+    }
+
+    nhi = FLINT_MIN(nhi, len1 + len2 - 1);
+    len = nhi - nlo;
+
+    if ((res == poly1) || (res == poly2))
+    {
+        fmpz *t = _fmpz_vec_init(len);
+        _fmpz_mod_poly_mulmid(t, poly1->coeffs, len1, poly2->coeffs, len2, nlo, nhi, ctx);
+        _fmpz_vec_clear(res->coeffs, res->alloc);
+        res->coeffs = t;
+        res->alloc  = len;
+        res->length = len;
+        _fmpz_mod_poly_normalise(res);
+    }
+    else
+    {
+        fmpz_mod_poly_fit_length(res, len, ctx);
+        _fmpz_mod_poly_mulmid(res->coeffs, poly1->coeffs, len1, poly2->coeffs, len2, nlo, nhi, ctx);
+        _fmpz_mod_poly_set_length(res, len);
+        _fmpz_mod_poly_normalise(res);
+    }
+}

src/fmpz_mod_poly/test/main.c

@@ -63,6 +63,7 @@
 #include "t-mul.c"
 #include "t-mulhigh.c"
 #include "t-mullow.c"
+#include "t-mulmid.c"
 #include "t-mulmod.c"
 #include "t-mulmod_preinv.c"
 #include "t-neg.c"
@@ -139,6 +140,7 @@ test_struct tests[] =
     TEST_FUNCTION(fmpz_mod_poly_mul),
     TEST_FUNCTION(fmpz_mod_poly_mulhigh),
     TEST_FUNCTION(fmpz_mod_poly_mullow),
+    TEST_FUNCTION(fmpz_mod_poly_mulmid),
     TEST_FUNCTION(fmpz_mod_poly_mulmod),
     TEST_FUNCTION(fmpz_mod_poly_mulmod_preinv),
     TEST_FUNCTION(fmpz_mod_poly_neg),

src/fmpz_mod_poly/test/t-mulmid.c

@@ -0,0 +1,82 @@
+/*
+    Copyright (C) 2009 William Hart
+    Copyright (C) 2011 Sebastian Pancratz
+
+    This file is part of FLINT.
+
+    FLINT is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 3 of the License, or
+    (at your option) any later version.  See <https://www.gnu.org/licenses/>.
+*/
+
+#include "test_helpers.h"
+#include "fmpz.h"
+#include "fmpz_mod.h"
+#include "fmpz_mod_poly.h"
+
+TEST_FUNCTION_START(fmpz_mod_poly_mulmid, state)
+{
+    int i, result;
+    fmpz_mod_ctx_t ctx;
+
+    fmpz_mod_ctx_init_ui(ctx, 2);
+
+    /* Compare with mullow */
+    for (i = 0; i < 200 * flint_test_multiplier(); i++)
+    {
+        fmpz_t p;
+        fmpz_mod_poly_t a, b, c;
+        slong nlo, nhi;
+
+        fmpz_init(p);
+        fmpz_randtest_unsigned(p, state, 2 * FLINT_BITS);
+        fmpz_add_ui(p, p, 2);
+        fmpz_mod_ctx_set_modulus(ctx, p);
+
+        fmpz_mod_poly_init(a, ctx);
+        fmpz_mod_poly_init(b, ctx);
+        fmpz_mod_poly_init(c, ctx);
+        nlo = n_randint(state, 50);
+        nhi = n_randint(state, 50);
+        fmpz_mod_poly_randtest(b, state, 1 + n_randint(state, 50), ctx);
+        fmpz_mod_poly_randtest(c, state, 1 + n_randint(state, 50), ctx);
+
+        if (n_randint(state, 2))
+        {
+            fmpz_mod_poly_set(a, b, ctx);
+            fmpz_mod_poly_mulmid(a, a, c, nlo, nhi, ctx);
+        }
+        else if (n_randint(state, 2))
+        {
+            fmpz_mod_poly_set(a, c, ctx);
+            fmpz_mod_poly_mulmid(a, b, a, nlo, nhi, ctx);
+        }
+        else
+        {
+            fmpz_mod_poly_mulmid(a, b, c, nlo, nhi, ctx);
+        }
+
+        fmpz_mod_poly_mullow(b, b, c, nhi, ctx);
+        fmpz_mod_poly_shift_right(b, b, nlo, ctx);
+
+        result = (fmpz_mod_poly_equal(a, b, ctx));
+        if (!result)
+        {
+            flint_printf("FAIL:\n");
+            fmpz_mod_poly_print(a, ctx), flint_printf("\n\n");
+            fmpz_mod_poly_print(b, ctx), flint_printf("\n\n");
+            fflush(stdout);
+            flint_abort();
+        }
+
+        fmpz_mod_poly_clear(a, ctx);
+        fmpz_mod_poly_clear(b, ctx);
+        fmpz_mod_poly_clear(c, ctx);
+        fmpz_clear(p);
+    }
+
+    fmpz_mod_ctx_clear(ctx);
+
+    TEST_FUNCTION_END(state);
+}

src/gr/fmpz_mod.c

@@ -554,6 +554,20 @@ _gr_fmpz_mod_poly_mullow(fmpz * res,
     return GR_SUCCESS;
 }
 
+static int
+_gr_fmpz_mod_poly_mulmid(fmpz * res,
+    const fmpz * poly1, slong len1,
+    const fmpz * poly2, slong len2, slong nlo, slong nhi, gr_ctx_t ctx)
+{
+    if (len1 >= len2)
+        _fmpz_mod_poly_mulmid(res, poly1, len1, poly2, len2, nlo, nhi, FMPZ_MOD_CTX(ctx));
+    else
+        _fmpz_mod_poly_mulmid(res, poly2, len2, poly1, len1, nlo, nhi, FMPZ_MOD_CTX(ctx));
+
+    return GR_SUCCESS;
+}
+
+
 /* fixme: duplicate _fmpz_mod_poly methods for error handling */
 static int
 _gr_fmpz_mod_poly_divrem(fmpz * Q, fmpz * R, const fmpz * A, slong lenA,
@@ -852,6 +866,7 @@ gr_method_tab_input _fmpz_mod_methods_input[] =
     {GR_METHOD_VEC_DOT,         (gr_funcptr) _gr_fmpz_mod_vec_dot},
     {GR_METHOD_VEC_DOT_REV,     (gr_funcptr) _gr_fmpz_mod_vec_dot_rev},
     {GR_METHOD_POLY_MULLOW,     (gr_funcptr) _gr_fmpz_mod_poly_mullow},
+    {GR_METHOD_POLY_MULMID,     (gr_funcptr) _gr_fmpz_mod_poly_mulmid},
     {GR_METHOD_POLY_INV_SERIES, (gr_funcptr) _gr_fmpz_mod_poly_inv_series},
     {GR_METHOD_POLY_DIV_SERIES, (gr_funcptr) _gr_fmpz_mod_poly_div_series},
     {GR_METHOD_POLY_DIVREM,     (gr_funcptr) _gr_fmpz_mod_poly_divrem},

src/gr_poly.h

@@ -615,6 +615,7 @@ WARN_UNUSED_RESULT int gr_poly_compose_mod_brent_kung_precomp_preinv(gr_poly_t r
 /* Test functions */
 
 void _gr_poly_test_mullow(gr_method_poly_binary_trunc_op mullow_impl, gr_method_poly_binary_trunc_op mullow_ref, flint_rand_t state, slong iters, slong maxn, gr_ctx_t ctx);
+void _gr_poly_test_mulmid(gr_method_poly_binary_trunc2_op mulmid_impl, gr_method_poly_binary_trunc2_op mulmid_ref, flint_rand_t state, slong iters, slong maxn, gr_ctx_t ctx);
 void _gr_poly_test_divrem(gr_method_poly_binary_binary_op divrem_impl, flint_rand_t state, slong iters, slong maxn, gr_ctx_t ctx);
 void _gr_poly_test_div(gr_method_poly_binary_op div_impl, flint_rand_t state, slong iters, slong maxn, gr_ctx_t ctx);
 void _gr_poly_test_inv_series(gr_method_poly_unary_trunc_op inv_series_impl, flint_rand_t state, slong iters, slong maxn, gr_ctx_t ctx);

src/gr_poly/test_mulmid.c

@@ -0,0 +1,99 @@
+/*
+    Copyright (C) 2024, 2026 Fredrik Johansson
+
+    This file is part of FLINT.
+
+    FLINT is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 3 of the License, or
+    (at your option) any later version.  See <https://www.gnu.org/licenses/>.
+*/
+
+#include "gr_vec.h"
+#include "gr_poly.h"
+
+PUSH_OPTIONS
+OPTIMIZE_OSIZE
+
+void _gr_poly_test_mulmid(gr_method_poly_binary_trunc2_op mulmid_impl, gr_method_poly_binary_trunc2_op mulmid_ref,
+    flint_rand_t state, slong iters, slong maxn, gr_ctx_t ctx)
+{
+    slong iter;
+    gr_ctx_ptr given_ctx = ctx;
+
+    if (mulmid_ref == NULL)
+        mulmid_ref = (gr_method_poly_binary_trunc2_op) _gr_poly_mulmid_generic;
+
+    for (iter = 0; iter < iters; iter++)
+    {
+        gr_ptr A, B, C, D;
+        slong n1, n2, nhi, nlo;
+        int status = GR_SUCCESS;
+        gr_ctx_t my_ctx;
+        gr_ctx_struct * ctx;
+        int squaring;
+
+        if (given_ctx == NULL)
+        {
+            gr_ctx_init_random(my_ctx, state);
+            ctx = my_ctx;
+        }
+        else
+            ctx = given_ctx;
+
+        squaring = n_randint(state, 2);
+
+        n1 = 1 + n_randint(state, 1 + maxn);
+        n2 = squaring ? n1 : 1 + n_randint(state, 1 + maxn);
+        nhi = 1 + n_randint(state, 1 + maxn);
+        nhi = FLINT_MIN(nhi, n1 + n2 - 1);
+        nlo = n_randint(state, nhi);
+
+        A = gr_heap_init_vec(n1, ctx);
+        B = squaring ? A : gr_heap_init_vec(n2, ctx);
+        C = gr_heap_init_vec(nhi - nlo, ctx);
+        D = gr_heap_init_vec(nhi - nlo, ctx);
+
+        GR_MUST_SUCCEED(_gr_vec_randtest(A, state, n1, ctx));
+        if (!squaring)
+            GR_MUST_SUCCEED(_gr_vec_randtest(B, state, n2, ctx));
+        GR_MUST_SUCCEED(_gr_vec_randtest(C, state, nhi - nlo, ctx));
+
+#if 0
+        /* Useful for debugging */
+        slong i;
+        for (i = 0; i < n1; i++)
+            GR_IGNORE(gr_set_ui(GR_ENTRY(A, i, ctx->sizeof_elem), i + 1, ctx));
+        for (i = 0; i < n2; i++)
+            GR_IGNORE(gr_set_ui(GR_ENTRY(B, i, ctx->sizeof_elem), i + 1, ctx));
+#endif
+
+        status |= mulmid_impl(C, A, n1, B, n2, nlo, nhi, ctx);
+
+        if (status == GR_SUCCESS)
+            status |= mulmid_ref(D, A, n1, B, n2, nlo, nhi, ctx);
+
+        if (status == GR_SUCCESS && _gr_vec_equal(C, D, nhi - nlo, ctx) == T_FALSE)
+        {
+            flint_printf("FAIL: mulmid\n");
+            gr_ctx_println(ctx);
+            flint_printf("len1 = %wd, len2 = %wd, nlo = %wd, ni = %wd, squaring = %d\n\n", n1, n2, nlo, nhi, squaring);
+            flint_printf("A:\n"); _gr_vec_print(A, n1, ctx); flint_printf("\n\n");
+            flint_printf("B:\n"); _gr_vec_print(B, n2, ctx); flint_printf("\n\n");
+            flint_printf("C:\n"); _gr_vec_print(C, nhi - nlo, ctx); flint_printf("\n\n");
+            flint_printf("D:\n"); _gr_vec_print(D, nhi - nlo, ctx); flint_printf("\n\n");
+            flint_abort();
+        }
+
+        gr_heap_clear_vec(A, n1, ctx);
+        if (!squaring)
+            gr_heap_clear_vec(B, n2, ctx);
+        gr_heap_clear_vec(C, nhi - nlo, ctx);
+        gr_heap_clear_vec(D, nhi - nlo, ctx);
+
+        if (given_ctx == NULL)
+            gr_ctx_clear(ctx);
+    }
+}
+
+POP_OPTIONS

src/mpn_mod.h

@@ -219,6 +219,12 @@ int mpn_mod_mat_reduce_row(slong * column, gr_mat_t A, slong * P, slong * L, slo
 
 /* Polynomial algorithms */
 
+int _mpn_mod_poly_mulmid_classical(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong nlo, slong nhi, gr_ctx_t ctx);
+int _mpn_mod_poly_mulmid_karatsuba(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong nlo, slong nhi, slong cutoff, gr_ctx_t ctx);
+int _mpn_mod_poly_mulmid_KS(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong nlo, slong nhi, gr_ctx_t ctx);
+int _mpn_mod_poly_mulmid_fft_small(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong nlo, slong nhi, gr_ctx_t ctx);
+int _mpn_mod_poly_mulmid(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong nlo, slong nhi, gr_ctx_t ctx);
+
 int _mpn_mod_poly_mullow_classical(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong len, gr_ctx_t ctx);
 int _mpn_mod_poly_mullow_karatsuba(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong len, slong cutoff, gr_ctx_t ctx);
 int _mpn_mod_poly_mullow_KS(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong len, gr_ctx_t ctx);