Middle product for arb_poly, acb_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.

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,

src/acb_poly/mulmid.c

@@ -0,0 +1,90 @@
+/*
+    Copyright (C) 2012 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 <math.h>
+#include "acb_poly.h"
+
+void
+_acb_poly_mulmid(acb_ptr res,
+    acb_srcptr poly1, slong len1,
+    acb_srcptr poly2, slong len2, slong nlo, slong nhi, slong prec)
+{
+    if (len1 <= 7 || len2 <= 7 || nhi <= 7)
+    {
+        _acb_poly_mulmid_classical(res, poly1, len1, poly2, len2, nlo, nhi, prec);
+    }
+    else
+    {
+        slong cutoff;
+        double p;
+
+        if (prec <= 2 * FLINT_BITS)
+        {
+            cutoff = 110;
+        }
+        else
+        {
+            p = log(prec);
+
+            cutoff = 10000.0 / (p * p * p);
+            cutoff = FLINT_MIN(cutoff, 60);
+            if (poly1 == poly2 && prec >= 256)
+                cutoff *= 1.25;
+            if (poly1 == poly2 && prec >= 4096)
+                cutoff *= 1.25;
+            cutoff = FLINT_MAX(cutoff, 8);
+        }
+
+        /* todo: tuning copied from mullow; needs retuning for small nhi - nlo */
+        if (2 * FLINT_MIN(len1, len2) <= cutoff || nhi <= cutoff)
+            _acb_poly_mulmid_classical(res, poly1, len1, poly2, len2, nlo, nhi, prec);
+        else
+            _acb_poly_mulmid_transpose(res, poly1, len1, poly2, len2, nlo, nhi, 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)
+{
+    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(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(res->coeffs, poly1->coeffs, xlen,
+            poly2->coeffs, ylen, nlo, nhi, prec);
+    }
+
+    _acb_poly_set_length(res, zlen);
+    _acb_poly_normalise(res);
+}