more gr_poly division: basecase divexact, series divconquer

doc/source/gr_poly.rst

@@ -277,6 +277,8 @@ Division uses the Karp-Markstein algorithm.
 
 .. function:: int _gr_poly_div_series_newton(gr_ptr res, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, slong len, slong cutoff, gr_ctx_t ctx)
               int gr_poly_div_series_newton(gr_poly_t res, const gr_poly_t A, const gr_poly_t B, slong len, slong cutoff, gr_ctx_t ctx)
+              int _gr_poly_div_series_divconquer(gr_ptr res, gr_srcptr B, slong Blen, gr_srcptr A, slong Alen, slong len, slong cutoff, gr_ctx_t ctx);
+              int gr_poly_div_series_divconquer(gr_poly_t Q, const gr_poly_t A, const gr_poly_t B, slong len, slong cutoff, gr_ctx_t ctx);
               int _gr_poly_div_series_invmul(gr_ptr res, gr_srcptr B, slong Blen, gr_srcptr A, slong Alen, slong len, gr_ctx_t ctx)
               int gr_poly_div_series_invmul(gr_poly_t res, const gr_poly_t A, const gr_poly_t B, slong len, gr_ctx_t ctx)
               int _gr_poly_div_series_basecase_preinv1(gr_ptr Q, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, gr_srcptr Binv, slong len, gr_ctx_t ctx)
@@ -286,6 +288,25 @@ Division uses the Karp-Markstein algorithm.
               int _gr_poly_div_series(gr_ptr res, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, slong len, gr_ctx_t ctx)
               int gr_poly_div_series(gr_poly_t res, const gr_poly_t A, const gr_poly_t B, slong len, gr_ctx_t ctx)
 
+Exact division
+--------------------------------------------------------------------------------
+
+These functions compute a quotient `A / B` which is known to be exact
+(without remainder) in `R[x]` (or in `R[[x]] / x^n` in the case of series
+division). Given a nonexact division, they may return gibberish instead of
+returning an error flag.
+There are no ``preinv1`` versions because those are identical to the
+``div`` counterparts.
+
+.. function:: int _gr_poly_divexact_basecase_noinv(gr_ptr Q, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, gr_ctx_t ctx)
+              int _gr_poly_divexact_basecase(gr_ptr Q, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, gr_ctx_t ctx)
+              int gr_poly_divexact_basecase(gr_poly_t Q, const gr_poly_t A, const gr_poly_t B, gr_ctx_t ctx)
+
+.. function:: int _gr_poly_divexact_series_basecase_noinv(gr_ptr Q, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, slong len, gr_ctx_t ctx)
+              int _gr_poly_divexact_series_basecase(gr_ptr Q, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, slong len, gr_ctx_t ctx)
+              int gr_poly_divexact_series_basecase(gr_poly_t Q, const gr_poly_t A, const gr_poly_t B, slong len, gr_ctx_t ctx)
+
+
 Square roots
 --------------------------------------------------------------------------------
 

src/gr_poly.h

@@ -212,6 +212,8 @@ WARN_UNUSED_RESULT int gr_poly_inv_series(gr_poly_t Qinv, const gr_poly_t Q, slo
 
 WARN_UNUSED_RESULT int _gr_poly_div_series_newton(gr_ptr res, gr_srcptr B, slong Blen, gr_srcptr A, slong Alen, slong len, slong cutoff, gr_ctx_t ctx);
 WARN_UNUSED_RESULT int gr_poly_div_series_newton(gr_poly_t Q, const gr_poly_t A, const gr_poly_t B, slong len, slong cutoff, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int _gr_poly_div_series_divconquer(gr_ptr res, gr_srcptr B, slong Blen, gr_srcptr A, slong Alen, slong len, slong cutoff, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_poly_div_series_divconquer(gr_poly_t Q, const gr_poly_t A, const gr_poly_t B, slong len, slong cutoff, gr_ctx_t ctx);
 WARN_UNUSED_RESULT int _gr_poly_div_series_basecase_generic(gr_ptr Q, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, slong len, gr_ctx_t ctx);
 WARN_UNUSED_RESULT int gr_poly_div_series_basecase(gr_poly_t Q, const gr_poly_t A, const gr_poly_t B, slong len, gr_ctx_t ctx);
 WARN_UNUSED_RESULT int _gr_poly_div_series_basecase_preinv1(gr_ptr Q, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, gr_srcptr Binv, slong len, gr_ctx_t ctx);

src/gr_poly/div_basecase.c

@@ -85,7 +85,7 @@ _gr_poly_div_basecase(gr_ptr Q,
 
     /* constant is a unit; can multiply by inverse */
     if (status == GR_SUCCESS)
-        status |= _gr_poly_div_basecase_preinv1(Q, A, Alen, B, Blen, invB, ctx);
+        status = _gr_poly_div_basecase_preinv1(Q, A, Alen, B, Blen, invB, ctx);
     else
         status = _gr_poly_div_basecase_noinv(Q, A, Alen, B, Blen, ctx);
 

src/gr_poly/div_series.c

@@ -19,15 +19,24 @@ _gr_poly_div_series_generic(gr_ptr Q,
     gr_srcptr A, slong Alen,
     gr_srcptr B, slong Blen, slong len, gr_ctx_t ctx)
 {
+    int status;
+
     /* todo */
     if (FLINT_MIN(Blen, len) <= DEFAULT_CUTOFF || ctx->methods[GR_METHOD_POLY_MULLOW] == (gr_funcptr) _gr_poly_mullow_generic)
     {
-        return _gr_poly_div_series_basecase(Q, A, Alen, B, Blen, len, ctx);
+        status = _gr_poly_div_series_basecase(Q, A, Alen, B, Blen, len, ctx);
     }
     else
     {
-        return _gr_poly_div_series_newton(Q, A, Alen, B, Blen, len, DEFAULT_CUTOFF, ctx);
+        status = _gr_poly_div_series_newton(Q, A, Alen, B, Blen, len, DEFAULT_CUTOFF, ctx);
+
+        /* Newton requires invertible constant term of B; basecase and divide-and-conquer
+           may yet succeed without it. */
+        if (status == GR_DOMAIN)
+            status = _gr_poly_div_series_divconquer(Q, A, Alen, B, Blen, len, DEFAULT_CUTOFF, ctx);
     }
+
+    return status;
 }
 
 int

src/gr_poly/div_series_basecase.c

@@ -211,9 +211,9 @@ _gr_poly_div_series_basecase_generic(gr_ptr Q,
         status = gr_inv(q, B, ctx);
 
         if (status == GR_SUCCESS)
-            status |= _gr_poly_div_series_basecase_preinv1(Q, A, Alen, B, Blen, q, len, ctx);
+            status = _gr_poly_div_series_basecase_preinv1(Q, A, Alen, B, Blen, q, len, ctx);
         else
-            status |= _gr_poly_div_series_basecase_noinv(Q, A, Alen, B, Blen, len, ctx);
+            status = _gr_poly_div_series_basecase_noinv(Q, A, Alen, B, Blen, len, ctx);
 
         GR_TMP_CLEAR(q, ctx);
     }

src/gr_poly/div_series_divconquer.c

@@ -0,0 +1,81 @@
+/*
+    Copyright (C) 2010 Sebastian Pancratz
+    Copyright (C) 2019 William Hart
+    Copyright (C) 2014, 2021, 2023 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 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "gr_vec.h"
+#include "gr_poly.h"
+
+int
+_gr_poly_div_series_divconquer(gr_ptr res, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, slong len, slong cutoff, gr_ctx_t ctx)
+{
+    gr_ptr Arev, Brev;
+    slong Arevlen;
+    int status = GR_SUCCESS;
+
+    Alen = FLINT_MIN(Alen, len);
+    Blen = FLINT_MIN(Blen, len);
+
+    Arevlen = Blen + len - 1;
+
+    /* todo: algorithm without zero padding */
+    /* todo: shallow reversals */
+    GR_TMP_INIT_VEC(Arev, Arevlen, ctx);
+    GR_TMP_INIT_VEC(Brev, Blen, ctx);
+
+    status |= _gr_poly_reverse(Arev, A, Alen, Arevlen, ctx);
+    status |= _gr_poly_reverse(Brev, B, Blen, Blen, ctx);
+    status |= _gr_poly_div_divconquer(res, Arev, Arevlen, Brev, Blen, cutoff, ctx);
+    status |= _gr_poly_reverse(res, res, len, len, ctx);
+
+    GR_TMP_CLEAR_VEC(Arev, Arevlen, ctx);
+    GR_TMP_CLEAR_VEC(Brev, Blen, ctx);
+
+    return status;
+}
+
+int
+gr_poly_div_series_divconquer(gr_poly_t Q, const gr_poly_t A, const gr_poly_t B, slong len, slong cutoff, gr_ctx_t ctx)
+{
+    int status = GR_SUCCESS;
+
+    if (len == 0)
+        return gr_poly_zero(Q, ctx);
+
+    if (B->length == 0)
+        return GR_DOMAIN;
+
+    if (A->length == 0)
+    {
+        truth_t is_zero = gr_poly_is_zero(B, ctx);
+
+        if (is_zero == T_FALSE)
+            return gr_poly_zero(Q, ctx);
+
+        return GR_UNABLE;
+    }
+
+    if (Q == A || Q == B)
+    {
+        gr_poly_t t;
+        gr_poly_init(t, ctx);
+        status = gr_poly_div_series_divconquer(t, A, B, len, cutoff, ctx);
+        gr_poly_swap(Q, t, ctx);
+        gr_poly_clear(t, ctx);
+        return status;
+    }
+
+    gr_poly_fit_length(Q, len, ctx);
+    status = _gr_poly_div_series_divconquer(Q->coeffs, A->coeffs, A->length, B->coeffs, B->length, len, cutoff, ctx);
+    _gr_poly_set_length(Q, len, ctx);
+    _gr_poly_normalise(Q, ctx);
+    return status;
+}

src/gr_poly/divexact_basecase.c

@@ -0,0 +1,105 @@
+/*
+    Copyright (C) 2023 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 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "gr_vec.h"
+#include "gr_poly.h"
+
+int
+_gr_poly_divexact_basecase_noinv(gr_ptr Q,
+    gr_srcptr A, slong Alen,
+    gr_srcptr B, slong Blen, gr_ctx_t ctx)
+{
+    int status = GR_SUCCESS;
+    slong sz = ctx->sizeof_elem;
+    slong Qlen;
+    slong i, l;
+
+    if (Blen == 1)
+        return _gr_vec_div_scalar(Q, A, Alen, B, ctx);
+
+    Qlen = Alen - Blen + 1;
+    status = gr_divexact(GR_ENTRY(Q, Qlen - 1, sz), GR_ENTRY(A, Alen - 1, sz), GR_ENTRY(B, Blen - 1, sz), ctx);
+
+    for (i = 1; status == GR_SUCCESS && i < Qlen; i++)
+    {
+        l = FLINT_MIN(i, Blen - 1);
+        status |= _gr_vec_dot_rev(GR_ENTRY(Q, Qlen - 1 - i, sz), GR_ENTRY(A, Alen - 1 - i, sz), 1,
+            GR_ENTRY(B, Blen - 1 - l, sz), GR_ENTRY(Q, Qlen - i, sz), l, ctx);
+        status |= gr_divexact(GR_ENTRY(Q, Qlen - 1 - i, sz), GR_ENTRY(Q, Qlen - 1 - i, sz), GR_ENTRY(B, Blen - 1, sz), ctx);
+    }
+
+    return status;
+}
+
+int
+_gr_poly_divexact_basecase(gr_ptr Q,
+    gr_srcptr A, slong Alen,
+    gr_srcptr B, slong Blen, gr_ctx_t ctx)
+{
+    int status = GR_SUCCESS;
+    slong sz = ctx->sizeof_elem;
+    gr_ptr invB;
+
+    GR_TMP_INIT(invB, ctx);
+
+    /* todo: we sometimes want to keep dividing, e.g. over RR with small coefficient */
+    status = gr_inv(invB, GR_ENTRY(B, Blen - 1, sz), ctx);
+
+    /* constant is a unit; can multiply by inverse */
+    if (status == GR_SUCCESS)
+        status = _gr_poly_div_basecase_preinv1(Q, A, Alen, B, Blen, invB, ctx);
+    else
+        status = _gr_poly_divexact_basecase_noinv(Q, A, Alen, B, Blen, ctx);
+
+    GR_TMP_CLEAR(invB, ctx);
+
+    return status;
+}
+
+int
+gr_poly_divexact_basecase(gr_poly_t Q, const gr_poly_t A, const gr_poly_t B, gr_ctx_t ctx)
+{
+    slong Alen, Blen, Qlen;
+    int status = GR_SUCCESS;
+    slong sz = ctx->sizeof_elem;
+
+    Alen = A->length;
+    Blen = B->length;
+
+    if (Blen == 0)
+        return GR_DOMAIN;
+
+    if (gr_is_zero(GR_ENTRY(B->coeffs, Blen - 1, sz), ctx) != T_FALSE)
+        return GR_UNABLE;
+
+    if (Alen < Blen)
+        return gr_poly_zero(Q, ctx);
+
+    Qlen = Alen - Blen + 1;
+
+    if (Q == A || Q == B)
+    {
+        gr_poly_t t;
+        gr_poly_init2(t, Qlen, ctx);
+        status = _gr_poly_divexact_basecase(t->coeffs, A->coeffs, A->length, B->coeffs, B->length, ctx);
+        gr_poly_swap(Q, t, ctx);
+        gr_poly_clear(t, ctx);
+    }
+    else
+    {
+        gr_poly_fit_length(Q, Qlen, ctx);
+        status = _gr_poly_divexact_basecase(Q->coeffs, A->coeffs, A->length, B->coeffs, B->length, ctx);
+    }
+
+    _gr_poly_set_length(Q, Qlen, ctx);
+    _gr_poly_normalise(Q, ctx);
+    return status;
+}