Merge pull request #2613 from fredrik-johansson/mulmid2

Middle products for `arb_poly`, `acb_poly`, `nmod_poly`, `gr_poly`

Author: fredrik-johansson

doc/source/acb_poly.rst

@@ -347,6 +347,16 @@ Arithmetic
     If the same variable is passed for *A* and *B*, sets *C* to
     the square of *A*.
 
+.. function:: void _acb_poly_mulmid_transpose(acb_ptr z, acb_srcptr x, slong xlen, acb_srcptr y, slong ylen, slong nlo, slong nhi, slong prec)
+              void acb_poly_mulmid_transpose(acb_poly_t res, const acb_poly_t poly1, const acb_poly_t poly2, slong nlo, slong nhi, slong prec)
+              void _acb_poly_mulmid_classical(acb_ptr z, acb_srcptr x, slong xlen, acb_srcptr y, slong ylen, slong nlo, slong nhi, slong prec)
+              void acb_poly_mulmid_classical(acb_poly_t res, const acb_poly_t poly1, const acb_poly_t poly2, slong nlo, slong nhi, slong prec)
+              void _acb_poly_mulmid(acb_ptr z, acb_srcptr x, slong xlen, acb_srcptr y, slong ylen, slong nlo, slong nhi, slong prec)
+              void acb_poly_mulmid(acb_poly_t res, const acb_poly_t poly1, const acb_poly_t poly2, slong nlo, slong nhi, slong prec)
+
+    Analogous to *mullow* functions, but compute the product truncated
+    at length *nhi* and right-shifted by *nlo*.
+
 .. function:: void _acb_poly_inv_series(acb_ptr Qinv, acb_srcptr Q, slong Qlen, slong len, slong prec)
 
     Sets *{Qinv, len}* to the power series inverse of *{Q, Qlen}*. Uses Newton iteration.

doc/source/arb_poly.rst

@@ -343,6 +343,16 @@ Arithmetic
     If the same variable is passed for *A* and *B*, sets *C* to the
     square of *A*.
 
+.. function:: void _arb_poly_mulmid_block(arb_ptr z, arb_srcptr x, slong xlen, arb_srcptr y, slong ylen, slong nlo, slong nhi, slong prec)
+              void arb_poly_mulmid_block(arb_poly_t res, const arb_poly_t poly1, const arb_poly_t poly2, slong nlo, slong nhi, slong prec)
+              void _arb_poly_mulmid_classical(arb_ptr z, arb_srcptr x, slong xlen, arb_srcptr y, slong ylen, slong nlo, slong nhi, slong prec)
+              void arb_poly_mulmid_classical(arb_poly_t res, const arb_poly_t poly1, const arb_poly_t poly2, slong nlo, slong nhi, slong prec)
+              void _arb_poly_mulmid(arb_ptr z, arb_srcptr x, slong xlen, arb_srcptr y, slong ylen, slong nlo, slong nhi, slong prec)
+              void arb_poly_mulmid(arb_poly_t res, const arb_poly_t poly1, const arb_poly_t poly2, slong nlo, slong nhi, slong prec)
+
+    Analogous to *mullow* functions, but compute the product truncated
+    at length *nhi* and right-shifted by *nlo*.
+
 .. function:: void _arb_poly_inv_series(arb_ptr Q, arb_srcptr A, slong Alen, slong len, slong prec)
 
     Sets *{Q, len}* to the power series inverse of *{A, Alen}*. Uses Newton iteration.

doc/source/gr_poly.rst

@@ -183,6 +183,10 @@ Arithmetic
     by default, which currently always delegates to :func:`_gr_poly_mullow_classical`.
     This can be overridden by specific rings.
 
+.. function:: int _gr_poly_mulmid(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, slong nlo, slong nhi, gr_ctx_t ctx)
+              int _gr_poly_mulmid_generic(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, slong nlo, slong nhi, gr_ctx_t ctx)
+              int gr_poly_mulmid(gr_poly_t res, const gr_poly_t poly1, const gr_poly_t poly2, slong nlo, slong nhi, gr_ctx_t ctx)
+
 Multiplication algorithms
 -------------------------------------------------------------------------------
 

doc/source/nmod_poly.rst

@@ -554,12 +554,12 @@ Bit packing and unpacking
     coefficient of ``poly`` is bigger than ``bits/2`` bits. We
     also assume ``bits < 3 * FLINT_BITS``.
 
-.. function:: void _nmod_poly_bit_unpack(nn_ptr res, slong len, nn_srcptr mpn, ulong bits, nmod_t mod)
+.. function:: void _nmod_poly_bit_unpack(nn_ptr res, slong nlo, slong nhi, nn_srcptr mpn, ulong bits, nmod_t mod)
 
-    Unpacks ``len`` coefficients stored in the big integer ``mpn``
-    in bit fields of the given number of bits, reduces them modulo the
-    given modulus, then stores them in the polynomial ``res``.
-    We assume ``len > 0`` and ``3 * FLINT_BITS > bits > 0``.
+    Unpacks ``nhi - nlo`` coefficients stored in the big integer ``mpn``
+    in bit fields of the given number of bits, starting at offset ``nlo``,
+    reduces them modulo the given modulus, then stores them in the polynomial ``res``.
+    We assume ``3 * FLINT_BITS > bits > 0``.
     There are no restrictions on the size of the actual coefficients as
     stored within the bitfields.
 
@@ -679,19 +679,14 @@ Multiplication
     coefficients from ``start`` onwards into the high coefficients of
     ``res``, the remaining coefficients being arbitrary but reduced.
 
-.. function:: void _nmod_poly_mul_KS(nn_ptr out, nn_srcptr in1, slong len1, nn_srcptr in2, slong len2, flint_bitcnt_t bits, nmod_t mod)
+.. function:: void _nmod_poly_mul_KS(nn_ptr out, nn_srcptr in1, slong len1, nn_srcptr in2, slong len2, nmod_t mod)
 
-    Sets ``res`` to the product of ``in1`` and ``in2``
-    assuming the output coefficients are at most the given number of
-    bits wide. If ``bits`` is set to `0` an appropriate value is
-    computed automatically.  Assumes that ``len1 >= len2 > 0``.
+    Sets ``res`` to the product of ``in1`` and ``in2``.
+    Assumes that ``len1 >= len2 > 0``.
 
-.. function:: void nmod_poly_mul_KS(nmod_poly_t res, const nmod_poly_t poly1, const nmod_poly_t poly2, flint_bitcnt_t bits)
+.. function:: void nmod_poly_mul_KS(nmod_poly_t res, const nmod_poly_t poly1, const nmod_poly_t poly2)
 
-    Sets ``res`` to the product of ``poly1`` and ``poly2``
-    assuming the output coefficients are at most the given number of
-    bits wide. If ``bits`` is set to `0` an appropriate value
-    is computed automatically.
+    Sets ``res`` to the product of ``poly1`` and ``poly2``.
 
 .. function:: void _nmod_poly_mul_KS2(nn_ptr res, nn_srcptr op1, slong n1, nn_srcptr op2, slong n2, nmod_t mod)
 
@@ -711,14 +706,14 @@ Multiplication
 
     Sets ``res`` to the product of ``poly1`` and ``poly2``.
 
-.. function:: void _nmod_poly_mullow_KS(nn_ptr out, nn_srcptr in1, slong len1, nn_srcptr in2, slong len2, flint_bitcnt_t bits, slong n, nmod_t mod)
+.. function:: void _nmod_poly_mullow_KS(nn_ptr out, nn_srcptr in1, slong len1, nn_srcptr in2, slong len2, slong n, nmod_t mod)
 
     Sets ``out`` to the low `n` coefficients of ``in1`` of length
     ``len1`` times ``in2`` of length ``len2``. The output must have
     space for ``n`` coefficients. We assume that ``len1 >= len2 > 0``
     and that ``0 < n <= len1 + len2 - 1``.
 
-.. function:: void nmod_poly_mullow_KS(nmod_poly_t res, const nmod_poly_t poly1, const nmod_poly_t poly2, flint_bitcnt_t bits, slong n)
+.. function:: void nmod_poly_mullow_KS(nmod_poly_t res, const nmod_poly_t poly1, const nmod_poly_t poly2, slong n)
 
     Set ``res`` to the low `n` coefficients of ``in1`` of length
     ``len1`` times ``in2`` of length ``len2``.
@@ -746,6 +741,20 @@ Multiplication
     Sets ``res`` to the first ``trunc`` coefficients of the
     product of ``poly1`` and ``poly2``.
 
+.. function:: void _nmod_poly_mulmid(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong nlo, slong nhi, nmod_t mod)
+              void nmod_poly_mulmid(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong nlo, slong nhi, nmod_t mod)
+              void _nmod_poly_mulmid_classical(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong nlo, slong nhi, nmod_t mod)
+              void nmod_poly_mulmid_classical(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong nlo, slong nhi, nmod_t mod)
+              void _nmod_poly_mulmid_KS(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong nlo, slong nhi, nmod_t mod)
+              void nmod_poly_mulmid_KS(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong nlo, slong nhi, nmod_t mod)
+
+    Sets ``res`` to the first ``nhi - nlo`` middle coefficients of the
+    product of ``poly1`` of length ``len1`` and ``poly2`` of
+    length ``len2`` starting at offset ``nlo``.
+    It is assumed that ``0 <= nlo < nhi <= len1 + len2 - 1``.
+    The function :func:`_nmod_poly_mulmid_classical` does not support
+    aliasing between inputs and outputs, but all others do.
+
 .. function:: void _nmod_poly_mulhigh(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong n, nmod_t mod)
 
     Sets all but the low `n` coefficients of ``res`` to the
@@ -767,16 +776,17 @@ Multiplication
     other multiplication algorithms, given inputs of length *len1* and *len2*
     and output truncation to length *n*.
 
-.. function:: int _nmod_poly_mullow_fft_small_repack(nn_ptr z, nn_srcptr a, slong an, nn_srcptr b, slong bn, slong zn, nmod_t mod)
+.. function:: int _nmod_poly_mulmid_fft_small_repack(nn_ptr z, nn_srcptr a, slong an, nn_srcptr b, slong bn, slong znlo, slong zn, nmod_t mod)
 
     Internal helper function for :func:`_nmod_poly_mullow_fft_small`: if the
     inputs are small enough to perform a repacked convolution of half the
     length, multiply and return 1, otherwise do nothing and return 0.
-    The conditions on the arguments are the same as for :func:`_nmod_poly_mullow`.
+    The conditions on the arguments are the same as for :func:`_nmod_poly_mulmid_fft_small`.
 
-.. function:: void _nmod_poly_mullow_fft_small(nn_ptr z, nn_srcptr a, slong an, nn_srcptr b, slong bn, slong zn, nmod_t mod)
+.. function:: void _nmod_poly_mulmid_fft_small(nn_ptr z, nn_srcptr a, slong an, nn_srcptr b, slong bn, slong znlo, slong zn, nmod_t mod)
+              void _nmod_poly_mullow_fft_small(nn_ptr z, nn_srcptr a, slong an, nn_srcptr b, slong bn, slong zn, nmod_t mod)
 
-    Low multiplication via the *fft_small* module. Throws an error
+    Multiplication via the *fft_small* module. Throws an error
     if *fft_small* is not available. The conditions on the arguments
     are the same as for :func:`_nmod_poly_mullow`.
 

src/acb_poly.h

@@ -333,6 +333,13 @@ void _acb_poly_mul(acb_ptr C,
 void acb_poly_mul(acb_poly_t res, const acb_poly_t poly1,
               const acb_poly_t poly2, slong prec);
 
+void _acb_poly_mulmid_transpose(acb_ptr z, acb_srcptr x, slong xlen, acb_srcptr y, slong ylen, slong nlo, slong nhi, slong prec);
+void acb_poly_mulmid_transpose(acb_poly_t res, const acb_poly_t poly1, const acb_poly_t poly2, slong nlo, slong nhi, slong prec);
+void _acb_poly_mulmid_classical(acb_ptr z, acb_srcptr x, slong xlen, acb_srcptr y, slong ylen, slong nlo, slong nhi, slong prec);
+void acb_poly_mulmid_classical(acb_poly_t res, const acb_poly_t poly1, const acb_poly_t poly2, slong nlo, slong nhi, slong prec);
+void _acb_poly_mulmid(acb_ptr z, acb_srcptr x, slong xlen, acb_srcptr y, slong ylen, slong nlo, slong nhi, slong prec);
+void acb_poly_mulmid(acb_poly_t res, const acb_poly_t poly1, const acb_poly_t poly2, slong nlo, slong nhi, slong prec);
+
 ACB_POLY_INLINE void
 _acb_poly_mul_monic(acb_ptr res, acb_srcptr poly1, slong len1,
     acb_srcptr poly2, slong len2, slong prec)

src/acb_poly/exp_series.c

@@ -49,16 +49,16 @@ _acb_poly_exp_series_newton(acb_ptr f, acb_ptr g,
     slong l = m - 1; /* shifted for derivative */
 
     /* g := exp(-h) + O(x^m) */
-    _acb_poly_mullow(T, f, m, g, m2, m, prec);
-    _acb_poly_mullow(g + m2, g, m2, T + m2, m - m2, m - m2, prec);
+    _acb_poly_mulmid(T, f, m, g, m2, m2, m, prec);
+    _acb_poly_mullow(g + m2, g, m2, T, m - m2, m - m2, prec);
     _acb_vec_neg(g + m2, g + m2, m - m2);
 
     /* U := h' + g (f' - f h') + O(x^(n-1))
         Note: should replace h' by h' mod x^(m-1) */
     _acb_vec_zero(f + m, n - m);
-    _acb_poly_mullow(T, f, n, hprime, n, n, prec); /* should be mulmid */
+    _acb_poly_mulmid(T, f, n, hprime, n, l, n, prec);
     _acb_poly_derivative(U, f, n, prec); acb_zero(U + n - 1); /* should skip low terms */
-    _acb_vec_sub(U + l, U + l, T + l, n - l, prec);
+    _acb_vec_sub(U + l, U + l, T, n - l, prec);
     _acb_poly_mullow(T + l, g, n - m, U + l, n - m, n - m, prec);
     _acb_vec_add(U + l, hprime + l, T + l, n - m, prec);
 
@@ -71,8 +71,8 @@ _acb_poly_exp_series_newton(acb_ptr f, acb_ptr g,
     /* not needed if we only want exp(x) */
     if (n == len && inverse)
     {
-        _acb_poly_mullow(T, f, n, g, m, n, prec);
-        _acb_poly_mullow(g + m, g, m, T + m, n - m, n - m, prec);
+        _acb_poly_mulmid(T, f, n, g, m, m, n, prec);
+        _acb_poly_mullow(g + m, g, m, T, n - m, n - m, prec);
         _acb_vec_neg(g + m, g + m, n - m);
     }
 

src/acb_poly/inv_series.c

@@ -12,12 +12,6 @@
 #include "arb_poly.h"
 #include "acb_poly.h"
 
-#define MULLOW(z, x, xn, y, yn, nn, prec) \
-    if ((xn) >= (yn)) \
-        _acb_poly_mullow(z, x, xn, y, yn, nn, prec); \
-    else \
-        _acb_poly_mullow(z, y, yn, x, xn, nn, prec); \
-
 void
 _acb_poly_inv_series(acb_ptr Qinv,
     acb_srcptr Q, slong Qlen, slong len, slong prec)
@@ -60,22 +54,22 @@ _acb_poly_inv_series(acb_ptr Qinv,
             slong Qnlen, Wlen, W2len;
             acb_ptr W;
 
-            W = _acb_vec_init(len);
+            W = _acb_vec_init(len / 2);
 
             NEWTON_INIT(blen, len)
             NEWTON_LOOP(m, n)
 
             Qnlen = FLINT_MIN(Qlen, n);
             Wlen = FLINT_MIN(Qnlen + m - 1, n);
             W2len = Wlen - m;
-            MULLOW(W, Q, Qnlen, Qinv, m, Wlen, prec);
-            MULLOW(Qinv + m, Qinv, m, W + m, W2len, n - m, prec);
+            _acb_poly_mulmid(W, Q, Qnlen, Qinv, m, m, Wlen, prec);
+            _acb_poly_mullow(Qinv + m, Qinv, m, W, W2len, n - m, prec);
             _acb_vec_neg(Qinv + m, Qinv + m, n - m);
 
             NEWTON_END_LOOP
             NEWTON_END
 
-            _acb_vec_clear(W, len);
+            _acb_vec_clear(W, len / 2);
         }
     }
 }

src/acb_poly/mullow_transpose.c

@@ -13,27 +13,47 @@
 #include "acb_poly.h"
 
 void
-_acb_poly_mullow_transpose(acb_ptr res,
+_acb_poly_mulmid_transpose(acb_ptr res,
     acb_srcptr poly1, slong len1,
-    acb_srcptr poly2, slong len2, slong n, slong prec)
+    acb_srcptr poly2, slong len2, slong nlo, slong nhi, slong prec)
 {
     arb_ptr a, b, c, d, e, f, w;
     arb_ptr t;
     slong i;
 
-    len1 = FLINT_MIN(len1, n);
-    len2 = FLINT_MIN(len2, n);
+    len1 = FLINT_MIN(len1, nhi);
+    len2 = FLINT_MIN(len2, nhi);
+
+    slong nlo2 = (len1 + len2 - 1) - nlo;
 
-    w = flint_malloc(sizeof(arb_struct) * (2 * (len1 + len2 + n)));
+    if (len1 > nlo2)
+    {
+        slong trunc = len1 - nlo2;
+        poly1 += trunc;
+        len1 -= trunc;
+        nlo -= trunc;
+        nhi -= trunc;
+    }
+
+    if (len2 > nlo2)
+    {
+        slong trunc = len2 - nlo2;
+        poly2 += trunc;
+        len2 -= trunc;
+        nlo -= trunc;
+        nhi -= trunc;
+    }
+
+    w = flint_malloc(sizeof(arb_struct) * (2 * (len1 + len2 + (nhi - nlo))));
     a = w;
     b = a + len1;
     c = b + len1;
     d = c + len2;
     e = d + len2;
-    f = e + n;
+    f = e + (nhi - nlo);
 
     /* (e+fi) = (a+bi)(c+di) = (ac - bd) + (ad + bc)i */
-    t = _arb_vec_init(n);
+    t = _arb_vec_init(nhi - nlo);
 
     for (i = 0; i < len1; i++)
     {
@@ -47,38 +67,84 @@ _acb_poly_mullow_transpose(acb_ptr res,
         d[i] = *acb_imagref(poly2 + i);
     }
 
-    for (i = 0; i < n; i++)
+    for (i = 0; i < nhi - nlo; i++)
     {
         e[i] = *acb_realref(res + i);
         f[i] = *acb_imagref(res + i);
     }
 
-    _arb_poly_mullow(e, a, len1, c, len2, n, prec);
-    _arb_poly_mullow(t, b, len1, d, len2, n, prec);
-    _arb_vec_sub(e, e, t, n, prec);
+    _arb_poly_mulmid(e, a, len1, c, len2, nlo, nhi, prec);
+    _arb_poly_mulmid(t, b, len1, d, len2, nlo, nhi, prec);
+    _arb_vec_sub(e, e, t, nhi - nlo, prec);
 
-    _arb_poly_mullow(f, a, len1, d, len2, n, prec);
+    _arb_poly_mulmid(f, a, len1, d, len2, nlo, nhi, prec);
     /* squaring */
     if (poly1 == poly2 && len1 == len2)
     {
-        _arb_vec_scalar_mul_2exp_si(f, f, n, 1);
+        _arb_vec_scalar_mul_2exp_si(f, f, nhi - nlo, 1);
     }
     else
     {
-        _arb_poly_mullow(t, b, len1, c, len2, n, prec);
-        _arb_vec_add(f, f, t, n, prec);
+        _arb_poly_mulmid(t, b, len1, c, len2, nlo, nhi, prec);
+        _arb_vec_add(f, f, t, nhi - nlo, prec);
     }
 
-    for (i = 0; i < n; i++)
+    for (i = 0; i < nhi - nlo; i++)
     {
         *acb_realref(res + i) = e[i];
         *acb_imagref(res + i) = f[i];
     }
 
-    _arb_vec_clear(t, n);
+    _arb_vec_clear(t, nhi - nlo);
     flint_free(w);
 }
 
+void
+_acb_poly_mullow_transpose(acb_ptr res,
+    acb_srcptr poly1, slong len1,
+    acb_srcptr poly2, slong len2, slong n, slong prec)
+{
+    _acb_poly_mulmid_transpose(res, poly1, len1, poly2, len2, 0, n, prec);
+}
+
+void
+acb_poly_mulmid_transpose(acb_poly_t res, const acb_poly_t poly1,
+              const acb_poly_t poly2, slong nlo, slong nhi, slong prec)
+{
+    slong xlen, ylen, zlen;
+
+    xlen = poly1->length;
+    ylen = poly2->length;
+
+    if (xlen == 0 || ylen == 0 || nlo >= FLINT_MIN(nhi, xlen + ylen - 1))
+    {
+        acb_poly_zero(res);
+        return;
+    }
+
+    nhi = FLINT_MIN(nhi, xlen + ylen - 1);
+    zlen = nhi - nlo;
+
+    if (res == poly1 || res == poly2)
+    {
+        acb_poly_t tmp;
+        acb_poly_init2(tmp, zlen);
+        _acb_poly_mulmid_transpose(tmp->coeffs, poly1->coeffs, xlen,
+            poly2->coeffs, ylen, nlo, nhi, prec);
+        acb_poly_swap(res, tmp);
+        acb_poly_clear(tmp);
+    }
+    else
+    {
+        acb_poly_fit_length(res, zlen);
+        _acb_poly_mulmid_transpose(res->coeffs, poly1->coeffs, xlen,
+            poly2->coeffs, ylen, nlo, nhi, prec);
+    }
+
+    _acb_poly_set_length(res, zlen);
+    _acb_poly_normalise(res);
+}
+
 void
 acb_poly_mullow_transpose(acb_poly_t res, const acb_poly_t poly1,
                                             const acb_poly_t poly2,