diff --git a/doc/source/gr_ore_poly.rst b/doc/source/gr_ore_poly.rst index f294571c1c..9b281fa218 100644 --- a/doc/source/gr_ore_poly.rst +++ b/doc/source/gr_ore_poly.rst @@ -319,6 +319,19 @@ 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)`. + +.. function:: int _gr_ore_poly_div(gr_ptr Q, gr_srcptr U, slong lenU, gr_srcptr V, slong lenV, gr_ore_poly_ctx_t ctx) + + Version of the divrem function which outputs only the quotient. + +.. function:: int _gr_ore_poly_rem(gr_ptr R, gr_srcptr U, slong lenU, gr_srcptr V, slong lenV, gr_ore_poly_ctx_t ctx) + + Version of the divrem function which outputs only the remainder. + .. raw:: latex \newpage diff --git a/src/gr_ore_poly.h b/src/gr_ore_poly.h index 4a2df9cb3c..96e4e223e3 100644 --- a/src/gr_ore_poly.h +++ b/src/gr_ore_poly.h @@ -280,6 +280,11 @@ 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); +WARN_UNUSED_RESULT int gr_ore_poly_div(gr_ore_poly_t Q, const gr_ore_poly_t U, gr_ore_poly_t V, gr_ore_poly_ctx_t ctx); +WARN_UNUSED_RESULT int gr_ore_poly_rem(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 diff --git a/src/gr_ore_poly/ctx.c b/src/gr_ore_poly/ctx.c index 627000c3c9..0bdd9d80d8 100644 --- a/src/gr_ore_poly/ctx.c +++ b/src/gr_ore_poly/ctx.c @@ -296,11 +296,11 @@ gr_method_tab_input _gr_ore_poly_methods_input[] = {GR_METHOD_POW_FMPZ, (gr_funcptr) gr_ore_poly_pow_fmpz}, {GR_METHOD_DIV, (gr_funcptr) gr_ore_poly_div}, {GR_METHOD_INV, (gr_funcptr) gr_ore_poly_inv}, - - {GR_METHOD_EUCLIDEAN_DIV, (gr_funcptr) gr_ore_poly_euclidean_div}, - {GR_METHOD_EUCLIDEAN_REM, (gr_funcptr) gr_ore_poly_euclidean_rem}, - {GR_METHOD_EUCLIDEAN_DIVREM, (gr_funcptr) gr_ore_poly_euclidean_divrem}, - +*/ + {GR_METHOD_EUCLIDEAN_DIV, (gr_funcptr) gr_ore_poly_div}, + {GR_METHOD_EUCLIDEAN_REM, (gr_funcptr) gr_ore_poly_rem}, + {GR_METHOD_EUCLIDEAN_DIVREM, (gr_funcptr) gr_ore_poly_divrem}, +/* {GR_METHOD_GCD, (gr_funcptr) gr_ore_poly_gcd}, {GR_METHOD_FACTOR, (gr_funcptr) gr_ore_poly_factor}, diff --git a/src/gr_ore_poly/divrem.c b/src/gr_ore_poly/divrem.c new file mode 100644 index 0000000000..b144e290e4 --- /dev/null +++ b/src/gr_ore_poly/divrem.c @@ -0,0 +1,159 @@ +/* + 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 . +*/ + +#include "flint.h" +#include "gr.h" +#include "gr_ore_poly.h" +#include "gr_vec.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; + + lenQ = lenU - lenV + 1; + lenR = lenU; + ordV = lenV - 1; + ordR = lenR - 1; + + // Set Q to 0 + status |= _gr_vec_zero(Q, lenQ, cctx); + + // Set R to U + status |= _gr_vec_set(R, U, lenU, cctx); + + gr_ptr c; + GR_TMP_INIT(c, cctx); + + gr_ptr A, B; + GR_TMP_INIT_VEC(A, lenQ, cctx); + GR_TMP_INIT_VEC(B, lenU, cctx); + + gr_ptr sigma_pows; + GR_TMP_INIT_VEC(sigma_pows, lenQ, cctx); + status |= gr_set(GR_ENTRY(sigma_pows, 0, el_size), GR_ENTRY(V, ordV, el_size), cctx); + for (slong i = 1; i < lenQ; i++) + status |= gr_ore_poly_sigma(GR_ENTRY(sigma_pows, i, el_size), GR_ENTRY(sigma_pows, i - 1, el_size), ctx); + + while (ordR > ordV) + { + slong k = ordR - ordV; + + // c = lc(R) / sigma ^ k (lc(V)) using sigma_pows computed above + status |= gr_div(c, GR_ENTRY(R, ordR, el_size), GR_ENTRY(sigma_pows, k, el_size), 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; + status |= _gr_vec_zero(A, lenA, 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; + status |= _gr_vec_sub(R, R, B, lenB, 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_TMP_CLEAR_VEC(sigma_pows, lenQ, cctx); + GR_TMP_CLEAR(c, cctx); + + GR_TMP_CLEAR_VEC(A, lenQ, cctx); + GR_TMP_CLEAR_VEC(B, lenU, cctx); + + 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) + { + status |= gr_ore_poly_set(R, U, ctx); + status |= gr_ore_poly_zero(Q, ctx); + return status; + } + + slong lenQ = lenU - lenV + 1; + + if (Q == U || Q == V || R == U || R == V || Q == R) // 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; +} + +int gr_ore_poly_div(gr_ore_poly_t Q, const gr_ore_poly_t U, gr_ore_poly_t V, gr_ore_poly_ctx_t ctx) +{ + gr_ore_poly_t R; + int status; + gr_ore_poly_init(R, ctx); + status = gr_ore_poly_divrem(Q, R, U, V, ctx); + gr_ore_poly_clear(R, ctx); + return status; +} + +int gr_ore_poly_rem(gr_ore_poly_t R, const gr_ore_poly_t U, gr_ore_poly_t V, gr_ore_poly_ctx_t ctx) +{ + gr_ore_poly_t Q; + int status; + gr_ore_poly_init(Q, ctx); + status = gr_ore_poly_divrem(Q, R, U, V, ctx); + gr_ore_poly_clear(Q, ctx); + return status; +} diff --git a/src/gr_ore_poly/test/main.c b/src/gr_ore_poly/test/main.c index fbe884c835..79e983f851 100644 --- a/src/gr_ore_poly/test/main.c +++ b/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 *************************************************************/ diff --git a/src/gr_ore_poly/test/t-divrem.c b/src/gr_ore_poly/test/t-divrem.c new file mode 100644 index 0000000000..82ca288701 --- /dev/null +++ b/src/gr_ore_poly/test/t-divrem.c @@ -0,0 +1,103 @@ +/* + 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 . +*/ + +#include "test_helpers.h" +#include "ulong_extras.h" +#include "gr_ore_poly.h" + +TEST_GR_FUNCTION_START(gr_ore_poly_divrem, state, count_success, count_domain, count_unable) +{ + slong iter; + 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(); + } + } + } + + count_success += (status == GR_SUCCESS); + count_domain += ((status & GR_DOMAIN) != 0); + count_unable += ((status & GR_UNABLE) != 0); + + 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); + } + TEST_GR_FUNCTION_END(state, count_success, count_domain, count_unable); +}