Division with remainder for Ore Polynomials

doc/source/gr_ore_poly.rst

@@ -319,6 +319,11 @@ Arithmetic
     which must be two Ore polynomials in the Ore algebra *ctx*.
     The underscore method assumes *res != poly1*, *res != poly2* (no aliasing) and *len1* != 0, *len2* != 0.
 
+.. function:: int _gr_ore_poly_divrem(gr_ptr Q, gr_ptr R, gr_srcptr U, slong lenU, gr_srcptr V, slong lenV, gr_ore_poly_ctx_t ctx)
+              int gr_ore_poly_divrem(gr_ore_poly_t Q, gr_ore_poly_t R, const gr_ore_poly_t U, gr_ore_poly_t V, gr_ore_poly_ctx_t ctx)
+
+    Sets *(Q, R)* to the unique pair such that `U = QV + R` and `ord(R) < ord(V)`.
+
 .. raw:: latex
 
     \newpage

src/gr_ore_poly.h

@@ -280,6 +280,9 @@ WARN_UNUSED_RESULT int gr_ore_poly_lmul_gen(gr_ore_poly_t res, const gr_ore_poly
 WARN_UNUSED_RESULT int _gr_ore_poly_mul(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, gr_ore_poly_ctx_t ctx);
 WARN_UNUSED_RESULT int gr_ore_poly_mul(gr_ore_poly_t res, const gr_ore_poly_t poly1, const gr_ore_poly_t poly2, gr_ore_poly_ctx_t ctx);
 
+WARN_UNUSED_RESULT int _gr_ore_poly_divrem(gr_ptr Q, gr_ptr R, gr_srcptr U, slong lenU, gr_srcptr V, slong lenV, gr_ore_poly_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_ore_poly_divrem(gr_ore_poly_t Q, gr_ore_poly_t R, const gr_ore_poly_t U, gr_ore_poly_t V, gr_ore_poly_ctx_t ctx);
+
 #ifdef __cplusplus
 }
 #endif

src/gr_ore_poly/divrem.c

@@ -0,0 +1,160 @@
+/*
+    Copyright (C) 2026 Maria Neagoie
+
+    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 "flint.h"
+#include "gr.h"
+#include "gr_ore_poly.h"
+
+// Returns the unique pair (Q, R) such that U = QV + R and ord(R) < ord(V)
+int _gr_ore_poly_divrem(gr_ptr Q, gr_ptr R, gr_srcptr U, slong lenU, gr_srcptr V, slong lenV, gr_ore_poly_ctx_t ctx)
+{
+    gr_ctx_struct *cctx = GR_ORE_POLY_ELEM_CTX(ctx);
+    slong el_size = gr_ctx_sizeof_elem(cctx);
+    slong lenQ, lenR, ordV, ordR;
+    int status = GR_SUCCESS;
+
+    if (GR_ORE_POLY_CTX(ctx)->sigma_delta == NULL)
+        return GR_UNABLE;
+
+    lenQ = lenU - lenV + 1;
+    lenR = lenU;
+    ordV = lenV - 1;
+    ordR = lenR - 1;
+
+    // Set Q to 0
+    for (slong k = 0; k < lenQ; k++)
+        status |= gr_zero(GR_ENTRY(Q, k, el_size), cctx);
+
+    // Set R to U
+    for (slong k = 0; k < lenU; k++)
+        status |= gr_set(GR_ENTRY(R, k, el_size), GR_ENTRY(U, k, el_size), cctx);
+
+    gr_ptr lcR, lcV, denominator, c;
+    GR_TMP_INIT(lcR, cctx);
+    GR_TMP_INIT(lcV, cctx);
+    GR_TMP_INIT(denominator, cctx);
+    GR_TMP_INIT(c, cctx);
+
+    gr_ptr A = flint_malloc(lenQ * el_size);
+    gr_ptr B = flint_malloc(lenU * el_size);
+    _gr_vec_init(A, lenQ, cctx);
+    _gr_vec_init(B, lenU, cctx);
+
+    while (ordR > ordV)
+    {
+        slong k = ordR - ordV;
+
+        status |= gr_set(lcR, GR_ENTRY(R, ordR, el_size), cctx);
+        status |= gr_set(lcV, GR_ENTRY(V, ordV, el_size), cctx);
+
+        // Compute denominator = sigma ^ k (lc(V))
+        status |= gr_set(denominator, lcV, cctx);
+        for (slong i = 0; i < k; i++)
+            status |= gr_ore_poly_sigma(denominator, denominator, ctx);
+
+        // c = lc(R) / denominator
+        status |= gr_div(c, lcR, denominator, cctx);
+
+        // R -= c * x^k * V. We compute A =  c * x^k, then B = A * V
+        // A = c * x^k, so A[k] = c, rest 0
+        slong lenA = k + 1;
+        for (slong i = 0; i < lenA; i++)
+            status |= gr_zero(GR_ENTRY(A, i, el_size), cctx);
+        status |= gr_set(GR_ENTRY(A, k, el_size), c, cctx);
+
+        // B = A * V
+        status |= _gr_ore_poly_mul(B, A, lenA, V, lenV, ctx);
+
+        // R -= B
+        slong lenB = lenA + lenV - 1;
+        for (slong i = 0; i < lenB; i++)
+            status |= gr_sub(GR_ENTRY(R, i, el_size), GR_ENTRY(R, i, el_size), GR_ENTRY(B, i, el_size), cctx);
+
+        // Q += c * x^k, so Q[k] += c
+        status |= gr_add(GR_ENTRY(Q, k, el_size), GR_ENTRY(Q, k, el_size), c, cctx);
+
+        ordR--;
+    }
+
+    gr_clear(lcR, cctx);
+    gr_clear(lcV, cctx);
+    gr_clear(denominator, cctx);
+    gr_clear(c, cctx);
+
+    _gr_vec_clear(A, lenQ, cctx);
+    _gr_vec_clear(B, lenU, cctx);
+    flint_free(A);
+    flint_free(B);
+
+    return status;
+}
+
+int gr_ore_poly_divrem(gr_ore_poly_t Q, gr_ore_poly_t R, const gr_ore_poly_t U, gr_ore_poly_t V, gr_ore_poly_ctx_t ctx)
+{
+    int status = GR_SUCCESS;
+    const slong lenU = U->length;
+    const slong lenV = V->length;
+
+    if (lenV == 0) // division by zero polynomial
+        return GR_DOMAIN;
+
+    if (lenU == 0) // result is 0
+    {
+        status |= gr_ore_poly_zero(Q, ctx);
+        status |= gr_ore_poly_zero(R, ctx);
+        return status;
+    }
+
+    if (lenU < lenV)
+    {
+        gr_ore_poly_t tU;
+        // take care of aliasing case with temp
+        gr_ore_poly_init(tU, ctx);
+        status |= gr_ore_poly_set(tU, U, ctx);
+        status |= gr_ore_poly_zero(Q, ctx);
+        status |= gr_ore_poly_set(R, tU, ctx);
+        gr_ore_poly_clear(tU, ctx);
+        return status;
+    }
+
+    slong lenQ = lenU - lenV + 1;
+
+    if (Q == U || Q == V || R == U || R == V) // treat aliasing case separately
+    {
+        gr_ore_poly_t tQ, tR;
+        gr_ore_poly_init(tQ, ctx);
+        gr_ore_poly_init(tR, ctx);
+
+        gr_ore_poly_fit_length(tQ, lenQ, ctx);
+        gr_ore_poly_fit_length(tR, lenU, ctx);
+        status = _gr_ore_poly_divrem(tQ->coeffs, tR->coeffs, U->coeffs, lenU, V->coeffs, lenV, ctx);
+        _gr_ore_poly_set_length(tQ, lenQ, ctx);
+        _gr_ore_poly_set_length(tR, lenU, ctx);
+        _gr_ore_poly_normalise(tQ, ctx);
+        _gr_ore_poly_normalise(tR, ctx);
+
+        gr_ore_poly_swap(Q, tQ, ctx);
+        gr_ore_poly_swap(R, tR, ctx);
+        gr_ore_poly_clear(tQ, ctx);
+        gr_ore_poly_clear(tR, ctx);
+    }
+    else
+    {
+        gr_ore_poly_fit_length(Q, lenQ, ctx);
+        gr_ore_poly_fit_length(R, lenU, ctx);
+        status = _gr_ore_poly_divrem(Q->coeffs, R->coeffs, U->coeffs, lenU, V->coeffs, lenV, ctx);
+        _gr_ore_poly_set_length(Q, lenQ, ctx);
+        _gr_ore_poly_set_length(R, lenU, ctx);
+        _gr_ore_poly_normalise(Q, ctx);
+        _gr_ore_poly_normalise(R, ctx);
+    }
+    return status;
+}

src/gr_ore_poly/test/main.c

@@ -15,6 +15,7 @@
 #include "t-set_str.c"
 #include "t-sigma_delta.c"
 #include "t-mul.c"
+#include "t-divrem.c"
 
 /* Array of test functions ***************************************************/
 
@@ -24,6 +25,7 @@ test_struct tests[] =
     TEST_FUNCTION(gr_ore_poly_set_str),
     TEST_FUNCTION(gr_ore_poly_sigma_delta),
     TEST_FUNCTION(gr_ore_poly_mul),
+    TEST_FUNCTION(gr_ore_poly_divrem)
 };
 
 /* main function *************************************************************/

src/gr_ore_poly/test/t-divrem.c

@@ -0,0 +1,114 @@
+/*
+    Copyright (C) 2026 Maria Neagoie
+
+    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 "ulong_extras.h"
+#include "gr_ore_poly.h"
+
+TEST_FUNCTION_START(gr_ore_poly_divrem, state)
+{
+    slong iter;
+    // slong success = 0, domain = 0, unable = 0, combined = 0;
+    for (iter = 0; iter < 3000 * flint_test_multiplier(); iter++)
+    {
+        int status = GR_SUCCESS;
+        gr_ctx_t ctx;
+        gr_ore_poly_ctx_t ore_ctx;
+        gr_ore_poly_t A, B, Q, R, QBR;
+
+        gr_ore_poly_ctx_init_randtest2(ctx, ore_ctx, state);
+
+        gr_ore_poly_init(A, ore_ctx);
+        gr_ore_poly_init(B, ore_ctx);
+        gr_ore_poly_init(Q, ore_ctx);
+        gr_ore_poly_init(R, ore_ctx);
+        gr_ore_poly_init(QBR, ore_ctx);
+
+        status |= gr_ore_poly_randtest(A, state, 1 + n_randint(state, 6), ore_ctx);
+        status |= gr_ore_poly_randtest(B, state, 1 + n_randint(state, 6), ore_ctx);
+        status |= gr_ore_poly_randtest(Q, state, 1 + n_randint(state, 6), ore_ctx);
+        status |= gr_ore_poly_randtest(R, state, 1 + n_randint(state, 6), ore_ctx);
+
+        if (n_randint(state, 3) == 0)
+        {
+            status |= gr_ore_poly_mul(A, A, B, ore_ctx);
+            status |= gr_ore_poly_add(A, A, R, ore_ctx);
+        }
+
+        if (status == GR_SUCCESS)
+        {
+            /* test aliasing */
+            switch (n_randint(state, 5))
+            {
+                case 0:
+                    status |= gr_ore_poly_set(Q, A, ore_ctx);
+                    status |= gr_ore_poly_divrem(Q, R, Q, B, ore_ctx);
+                    break;
+                case 1:
+                    status |= gr_ore_poly_set(R, A, ore_ctx);
+                    status |= gr_ore_poly_divrem(Q, R, R, B, ore_ctx);
+                    break;
+                case 2:
+                    status |= gr_ore_poly_set(Q, B, ore_ctx);
+                    status |= gr_ore_poly_divrem(Q, R, A, Q, ore_ctx);
+                    break;
+                case 3:
+                    status |= gr_ore_poly_set(R, B, ore_ctx);
+                    status |= gr_ore_poly_divrem(Q, R, A, R, ore_ctx);
+                    break;
+                default:
+                    status |= gr_ore_poly_divrem(Q, R, A, B, ore_ctx);
+                    break;
+            }
+
+            if (status == GR_SUCCESS)
+            {
+                status |= gr_ore_poly_mul(QBR, Q, B, ore_ctx);
+                status |= gr_ore_poly_add(QBR, QBR, R, ore_ctx);
+
+                if (status == GR_SUCCESS && gr_ore_poly_equal(QBR, A, ore_ctx) == T_FALSE)
+                {
+                    flint_printf("FAIL\n\n");
+                    flint_printf("A = "); gr_ore_poly_print(A, ore_ctx); flint_printf("\n");
+                    flint_printf("B = "); gr_ore_poly_print(B, ore_ctx); flint_printf("\n");
+                    flint_printf("Q = "); gr_ore_poly_print(Q, ore_ctx); flint_printf("\n");
+                    flint_printf("R = "); gr_ore_poly_print(R, ore_ctx); flint_printf("\n");
+                    flint_printf("Q*B + R = "); gr_ore_poly_print(QBR, ore_ctx); flint_printf("\n");
+                    flint_abort();
+                }
+                // success++;
+            }
+        }
+
+        /*if (status == GR_DOMAIN)
+            domain++;
+
+        if (status == GR_UNABLE)
+            unable++;
+
+        if (status == 3)
+            combined++;*/
+
+        gr_ore_poly_clear(A, ore_ctx);
+        gr_ore_poly_clear(B, ore_ctx);
+        gr_ore_poly_clear(Q, ore_ctx);
+        gr_ore_poly_clear(R, ore_ctx);
+        gr_ore_poly_clear(QBR, ore_ctx);
+
+        gr_ore_poly_ctx_clear(ore_ctx);
+        gr_ctx_clear(ctx);
+    }
+    /*flint_printf("GR_SUCCESS = %d\n", success);
+    flint_printf("GR_DOMAIN = %d\n", domain);
+    flint_printf("GR_UNABLE = %d\n", unable);
+    flint_printf("Combined (3) = %d\n", combined);*/
+    TEST_FUNCTION_END(state);
+}