From ed877e7165fafc9a36c02601ed9e408aed659e22 Mon Sep 17 00:00:00 2001 From: "Marc Mezzarobba (AI)" Date: Tue, 30 Jun 2026 17:08:14 +0200 Subject: [PATCH 1/7] gr_ore_poly: apply an Ore polynomial to a base ring element Co-Authored-By: Claude Opus 4.8 --- doc/source/gr_ore_poly.rst | 14 ++ src/gr_ore_poly.h | 3 + src/gr_ore_poly/apply.c | 103 ++++++++++++ src/gr_ore_poly/test/main.c | 4 +- src/gr_ore_poly/test/t-apply.c | 285 +++++++++++++++++++++++++++++++++ 5 files changed, 408 insertions(+), 1 deletion(-) create mode 100644 src/gr_ore_poly/apply.c create mode 100644 src/gr_ore_poly/test/t-apply.c diff --git a/doc/source/gr_ore_poly.rst b/doc/source/gr_ore_poly.rst index 47c4c2fb09..6eed64ff3a 100644 --- a/doc/source/gr_ore_poly.rst +++ b/doc/source/gr_ore_poly.rst @@ -267,6 +267,20 @@ Action A pointer to a function with the same specification as :func:`gr_ore_poly_sigma_delta`. +.. function:: int gr_ore_poly_apply_custom(gr_ptr res, const gr_ore_poly_t P, gr_srcptr f, gr_srcptr d1, gr_ore_poly_ctx_t ctx) + + Sets *res* to the result of applying *P* to the base ring element *f*, where + the generator `D` acts by `g \mapsto \sigma(g) \cdot d1 + \delta(g)` for the + given value ``d1`` of `D(1)`. Any ``d1`` defines a valid action. Aliasing of + *res* with *f* or *d1* is allowed. + +.. function:: int gr_ore_poly_apply(gr_ptr res, const gr_ore_poly_t P, gr_srcptr f, gr_ore_poly_ctx_t ctx) + + As :func:`gr_ore_poly_apply_custom`, but with the standard `D(1)` for the + algebra type (so derivative operators differentiate, shift operators shift, + etc.). Returns ``GR_UNABLE`` when no standard action exists + (`ORE_ALGEBRA_COMMUTATIVE`, `ORE_ALGEBRA_CUSTOM`) or it cannot be computed. + Arithmetic ------------------------------------------------------------------------------- diff --git a/src/gr_ore_poly.h b/src/gr_ore_poly.h index cf189b74df..93d690a6f1 100644 --- a/src/gr_ore_poly.h +++ b/src/gr_ore_poly.h @@ -247,6 +247,9 @@ gr_ore_poly_delta(gr_ptr res, gr_srcptr a, gr_ore_poly_ctx_t ctx) extern const gr_ore_poly_sigma_delta_t _gr_ore_poly_default_sigma_delta[]; +WARN_UNUSED_RESULT int gr_ore_poly_apply_custom(gr_ptr res, const gr_ore_poly_t P, gr_srcptr f, gr_srcptr d1, gr_ore_poly_ctx_t ctx); +WARN_UNUSED_RESULT int gr_ore_poly_apply(gr_ptr res, const gr_ore_poly_t P, gr_srcptr f, gr_ore_poly_ctx_t ctx); + /* Arithmetic */ WARN_UNUSED_RESULT int gr_ore_poly_neg(gr_ore_poly_t res, const gr_ore_poly_t src, gr_ore_poly_ctx_t ctx); diff --git a/src/gr_ore_poly/apply.c b/src/gr_ore_poly/apply.c new file mode 100644 index 0000000000..1a98d1ac9e --- /dev/null +++ b/src/gr_ore_poly/apply.c @@ -0,0 +1,103 @@ +/* + 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 . +*/ + +/* generated using Claude Opus 4.8 */ + +#include "flint.h" +#include "gr.h" +#include "gr_ore_poly.h" + +int +gr_ore_poly_apply_custom(gr_ptr res, const gr_ore_poly_t P, gr_srcptr f, gr_srcptr d1, gr_ore_poly_ctx_t ctx) +{ + gr_ctx_struct * base = GR_ORE_POLY_ELEM_CTX(ctx); + slong el = base->sizeof_elem; + slong len = P->length; + int status = GR_SUCCESS; + truth_t d1_is_zero, d1_is_one; + gr_ptr g, acc, sig, del, term; + + if (len == 0) + return gr_zero(res, base); + + d1_is_zero = gr_is_zero(d1, base); + d1_is_one = gr_is_one(d1, base); + + GR_TMP_INIT5(g, acc, sig, del, term, base); + + status |= gr_set(g, f, base); + status |= gr_zero(acc, base); + + for (slong i = 0; i < len; i++) + { + /* acc += p_i * (D^i f) */ + status |= gr_mul(term, GR_ENTRY(P->coeffs, i, el), g, base); + status |= gr_add(acc, acc, term, base); + + /* g <- D(g) = sigma(g)*d1 + delta(g) for the next iteration */ + if (i + 1 < len) + { + if (d1_is_zero == T_TRUE) + { + status |= gr_ore_poly_delta(del, g, ctx); + status |= gr_set(g, del, base); + } + else + { + status |= gr_ore_poly_sigma_delta(sig, del, g, ctx); + if (d1_is_one != T_TRUE) + status |= gr_mul(sig, sig, d1, base); + status |= gr_add(g, sig, del, base); + } + } + } + + status |= gr_set(res, acc, base); + + GR_TMP_CLEAR5(g, acc, sig, del, term, base); + + return status; +} + +int +gr_ore_poly_apply(gr_ptr res, const gr_ore_poly_t P, gr_srcptr f, gr_ore_poly_ctx_t ctx) +{ + gr_ctx_struct * base = GR_ORE_POLY_ELEM_CTX(ctx); + int status = GR_SUCCESS; + gr_ptr d1; + + GR_TMP_INIT(d1, base); + + /* Standard value of D(1) for the algebra type. */ + switch (GR_ORE_POLY_CTX(ctx)->which_algebra) + { + case ORE_ALGEBRA_DERIVATIVE: + case ORE_ALGEBRA_EULER_DERIVATIVE: + case ORE_ALGEBRA_FORWARD_DIFFERENCE: + case ORE_ALGEBRA_BACKWARD_DIFFERENCE: + status |= gr_zero(d1, base); + break; + case ORE_ALGEBRA_FORWARD_SHIFT: + case ORE_ALGEBRA_BACKWARD_SHIFT: + case ORE_ALGEBRA_Q_SHIFT: + case ORE_ALGEBRA_MAHLER: + case ORE_ALGEBRA_FROBENIUS: + status |= gr_one(d1, base); + break; + default: /* ORE_ALGEBRA_COMMUTATIVE, ORE_ALGEBRA_CUSTOM: no standard action */ + status = GR_UNABLE; + } + + if (status == GR_SUCCESS) + status |= gr_ore_poly_apply_custom(res, P, f, d1, ctx); + + GR_TMP_CLEAR(d1, base); + + return status; +} diff --git a/src/gr_ore_poly/test/main.c b/src/gr_ore_poly/test/main.c index 79e983f851..838e5889c1 100644 --- a/src/gr_ore_poly/test/main.c +++ b/src/gr_ore_poly/test/main.c @@ -16,6 +16,7 @@ #include "t-sigma_delta.c" #include "t-mul.c" #include "t-divrem.c" +#include "t-apply.c" /* Array of test functions ***************************************************/ @@ -25,7 +26,8 @@ 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) + TEST_FUNCTION(gr_ore_poly_divrem), + TEST_FUNCTION(gr_ore_poly_apply) }; /* main function *************************************************************/ diff --git a/src/gr_ore_poly/test/t-apply.c b/src/gr_ore_poly/test/t-apply.c new file mode 100644 index 0000000000..f498e55964 --- /dev/null +++ b/src/gr_ore_poly/test/t-apply.c @@ -0,0 +1,285 @@ +/* + 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 "fmpz.h" +#include "test_helpers.h" +#include "gr_ore_poly.h" + +static void +check_gen_action(gr_srcptr x, gr_srcptr expected, gr_ore_poly_ctx_t octx, gr_ctx_t cctx) +{ + int status = GR_SUCCESS; + gr_ore_poly_t D; + gr_ptr got; + + gr_ore_poly_init(D, octx); + got = gr_heap_init(cctx); + + status |= gr_ore_poly_gen(D, octx); + status |= gr_ore_poly_apply(got, D, x, octx); + + if (status != GR_SUCCESS || gr_equal(got, expected, cctx) != T_TRUE) + { + flint_printf("FAIL: got %{gr}, expected %{gr}, in %{gr_ctx}\n", + got, cctx, expected, cctx, cctx); + flint_abort(); + } + + gr_heap_clear(got, cctx); + gr_ore_poly_clear(D, octx); +} + +TEST_GR_FUNCTION_START(gr_ore_poly_apply, state, count_success, count_domain, count_unable) +{ + { + int status = GR_SUCCESS; + gr_ctx_t zctx, cctx; + gr_ore_poly_ctx_t octx; + gr_ptr x, q, expected; + + gr_ctx_init_fmpz(zctx); + gr_ctx_init_gr_poly(cctx, zctx); + x = gr_heap_init(cctx); + q = gr_heap_init(cctx); + expected = gr_heap_init(cctx); + status |= gr_gen(x, cctx); + + gr_ore_poly_ctx_init(octx, cctx, 0, ORE_ALGEBRA_DERIVATIVE); + status |= gr_one(expected, cctx); + check_gen_action(x, expected, octx, cctx); + gr_ore_poly_ctx_clear(octx); + + gr_ore_poly_ctx_init(octx, cctx, 0, ORE_ALGEBRA_EULER_DERIVATIVE); + status |= gr_set(expected, x, cctx); + check_gen_action(x, expected, octx, cctx); + gr_ore_poly_ctx_clear(octx); + + gr_ore_poly_ctx_init(octx, cctx, 0, ORE_ALGEBRA_FORWARD_SHIFT); + status |= gr_add_si(expected, x, 1, cctx); + check_gen_action(x, expected, octx, cctx); + gr_ore_poly_ctx_clear(octx); + + gr_ore_poly_ctx_init(octx, cctx, 0, ORE_ALGEBRA_FORWARD_DIFFERENCE); + status |= gr_one(expected, cctx); + check_gen_action(x, expected, octx, cctx); + gr_ore_poly_ctx_clear(octx); + + gr_ore_poly_ctx_init(octx, cctx, 0, ORE_ALGEBRA_BACKWARD_SHIFT); + status |= gr_add_si(expected, x, -1, cctx); + check_gen_action(x, expected, octx, cctx); + gr_ore_poly_ctx_clear(octx); + + gr_ore_poly_ctx_init(octx, cctx, 0, ORE_ALGEBRA_BACKWARD_DIFFERENCE); + status |= gr_one(expected, cctx); + check_gen_action(x, expected, octx, cctx); + gr_ore_poly_ctx_clear(octx); + + status |= gr_set_si(q, 3, cctx); + status |= gr_ore_poly_ctx_init_q_shift(octx, cctx, 0, q); + status |= gr_mul(expected, q, x, cctx); + check_gen_action(x, expected, octx, cctx); + gr_ore_poly_ctx_clear(octx); + + status |= gr_ore_poly_ctx_init_mahler(octx, cctx, 0, 3); + status |= gr_pow_ui(expected, x, 3, cctx); + check_gen_action(x, expected, octx, cctx); + gr_ore_poly_ctx_clear(octx); + + if (status != GR_SUCCESS) + { + flint_printf("FAIL: unexpected failure in %{gr_ctx}\n", cctx, cctx); + flint_abort(); + } + + gr_heap_clear(x, cctx); + gr_heap_clear(q, cctx); + gr_heap_clear(expected, cctx); + gr_ctx_clear(cctx); + gr_ctx_clear(zctx); + + { + fmpz_t p; + gr_ctx_t fctx; + gr_ptr fx, fexpected; + + fmpz_init_set_ui(p, 17); + gr_ctx_init_fq(fctx, p, 3, NULL); + fmpz_clear(p); + + fx = gr_heap_init(fctx); + fexpected = gr_heap_init(fctx); + + status = GR_SUCCESS; + status |= gr_gen(fx, fctx); + status |= gr_pow_ui(fexpected, fx, 17, fctx); + + gr_ore_poly_ctx_init(octx, fctx, 0, ORE_ALGEBRA_FROBENIUS); + check_gen_action(fx, fexpected, octx, fctx); + gr_ore_poly_ctx_clear(octx); + + if (status != GR_SUCCESS) + { + flint_printf("FAIL: unexpected failure in %{gr_ctx}\n", fctx, fctx); + flint_abort(); + } + + gr_heap_clear(fx, fctx); + gr_heap_clear(fexpected, fctx); + gr_ctx_clear(fctx); + } + } + + + for (slong iter = 0; iter < 1000 * flint_test_multiplier(); iter++) + { + gr_ctx_t cctx, ctx; + slong maxlen; + int status = GR_SUCCESS; + + gr_ore_poly_ctx_init_randtest2(cctx, ctx, state); + + /* Repeated application of D blows up the degree for these, so keep the + operators short. */ + if (GR_ORE_POLY_CTX(ctx)->which_algebra == ORE_ALGEBRA_MAHLER + || GR_ORE_POLY_CTX(ctx)->which_algebra == ORE_ALGEBRA_Q_SHIFT) + maxlen = 2; + else + maxlen = 5; + + gr_ore_poly_t P, Q, PQ, PpQ, D; + gr_ore_poly_init(P, ctx); + gr_ore_poly_init(Q, ctx); + gr_ore_poly_init(PQ, ctx); + gr_ore_poly_init(PpQ, ctx); + gr_ore_poly_init(D, ctx); + + gr_ptr f = gr_heap_init(cctx); + gr_ptr g = gr_heap_init(cctx); + gr_ptr d1 = gr_heap_init(cctx); + gr_ptr u = gr_heap_init(cctx); + gr_ptr v = gr_heap_init(cctx); + gr_ptr w = gr_heap_init(cctx); + gr_ptr lhs = gr_heap_init(cctx); + gr_ptr rhs = gr_heap_init(cctx); + gr_ptr sf = gr_heap_init(cctx); + gr_ptr df = gr_heap_init(cctx); + gr_ptr one = gr_heap_init(cctx); + gr_ptr c = gr_heap_init(cctx); + + status |= gr_ore_poly_randtest(P, state, 1 + n_randint(state, maxlen), ctx); + status |= gr_ore_poly_randtest(Q, state, 1 + n_randint(state, maxlen), ctx); + status |= gr_ore_poly_gen(D, ctx); + status |= gr_randtest(f, state, cctx); + status |= gr_randtest(g, state, cctx); + status |= gr_randtest(d1, state, cctx); + + /* D(f) = sigma(f)*d1 + delta(f) */ + status |= gr_ore_poly_apply_custom(u, D, f, d1, ctx); + status |= gr_ore_poly_sigma_delta(sf, df, f, ctx); + status |= gr_mul(sf, sf, d1, cctx); + status |= gr_add(v, sf, df, cctx); + if (status == GR_SUCCESS && gr_equal(u, v, cctx) == T_FALSE) + { + flint_printf("FAIL: D(f) = sigma(f)*d1 + delta(f)\n"); + flint_abort(); + } + + /* (P*Q)(f) = P(Q(f)) */ + status |= gr_ore_poly_mul(PQ, P, Q, ctx); + status |= gr_ore_poly_apply_custom(lhs, PQ, f, d1, ctx); + status |= gr_ore_poly_apply_custom(u, Q, f, d1, ctx); + status |= gr_ore_poly_apply_custom(rhs, P, u, d1, ctx); + if (status == GR_SUCCESS && gr_equal(lhs, rhs, cctx) == T_FALSE) + { + flint_printf("FAIL: (P*Q)(f) = P(Q(f))\n"); + flint_abort(); + } + + /* (P+Q)(f) = P(f) + Q(f) */ + status |= gr_ore_poly_add(PpQ, P, Q, ctx); + status |= gr_ore_poly_apply_custom(lhs, PpQ, f, d1, ctx); + status |= gr_ore_poly_apply_custom(u, P, f, d1, ctx); + status |= gr_ore_poly_apply_custom(v, Q, f, d1, ctx); + status |= gr_add(rhs, u, v, cctx); + if (status == GR_SUCCESS && gr_equal(lhs, rhs, cctx) == T_FALSE) + { + flint_printf("FAIL: (P+Q)(f) = P(f) + Q(f)\n"); + flint_abort(); + } + + /* P(f+g) = P(f) + P(g) */ + status |= gr_add(w, f, g, cctx); + status |= gr_ore_poly_apply_custom(lhs, P, w, d1, ctx); + status |= gr_ore_poly_apply_custom(u, P, f, d1, ctx); + status |= gr_ore_poly_apply_custom(v, P, g, d1, ctx); + status |= gr_add(rhs, u, v, cctx); + if (status == GR_SUCCESS && gr_equal(lhs, rhs, cctx) == T_FALSE) + { + flint_printf("FAIL: P(f+g) = P(f) + P(g)\n"); + flint_abort(); + } + + /* aliasing res == f */ + status |= gr_ore_poly_apply_custom(u, P, f, d1, ctx); + status |= gr_set(v, f, cctx); + status |= gr_ore_poly_apply_custom(v, P, v, d1, ctx); + if (status == GR_SUCCESS && gr_equal(u, v, cctx) == T_FALSE) + { + flint_printf("FAIL: aliasing res == f\n"); + flint_abort(); + } + + /* aliasing res == d1 */ + status |= gr_ore_poly_apply_custom(u, P, f, d1, ctx); + status |= gr_set(v, d1, cctx); + status |= gr_ore_poly_apply_custom(v, P, f, v, ctx); + if (status == GR_SUCCESS && gr_equal(u, v, cctx) == T_FALSE) + { + flint_printf("FAIL: aliasing res == d1\n"); + flint_abort(); + } + + /* standard apply matches custom with c = D(1) = apply(D, 1) */ + status |= gr_one(one, cctx); + status |= gr_ore_poly_apply(c, D, one, ctx); + status |= gr_ore_poly_apply(lhs, P, f, ctx); + status |= gr_ore_poly_apply_custom(rhs, P, f, c, ctx); + if (status == GR_SUCCESS && gr_equal(lhs, rhs, cctx) == T_FALSE) + { + flint_printf("FAIL: apply(P, f) = apply_custom(P, f, D(1))\n"); + flint_abort(); + } + + count_success += (status == GR_SUCCESS); + count_domain += ((status & GR_DOMAIN) != 0); + count_unable += ((status & GR_UNABLE) != 0); + + gr_heap_clear(f, cctx); + gr_heap_clear(g, cctx); + gr_heap_clear(d1, cctx); + gr_heap_clear(u, cctx); + gr_heap_clear(v, cctx); + gr_heap_clear(w, cctx); + gr_heap_clear(lhs, cctx); + gr_heap_clear(rhs, cctx); + gr_heap_clear(sf, cctx); + gr_heap_clear(df, cctx); + gr_heap_clear(one, cctx); + gr_heap_clear(c, cctx); + gr_ore_poly_clear(P, ctx); + gr_ore_poly_clear(Q, ctx); + gr_ore_poly_clear(PQ, ctx); + gr_ore_poly_clear(PpQ, ctx); + gr_ore_poly_clear(D, ctx); + gr_ore_poly_ctx_clear(ctx); + gr_ctx_clear(cctx); + } + + TEST_GR_FUNCTION_END(state, count_success, count_domain, count_unable); +} From 28228cf7b36cc339b45fb6e4838218585b79bd91 Mon Sep 17 00:00:00 2001 From: "Marc Mezzarobba (AI)" Date: Tue, 30 Jun 2026 17:08:14 +0200 Subject: [PATCH 2/7] =?UTF-8?q?gr=5Fore=5Fpoly:=20conversions=20between=20?= =?UTF-8?q?d/dx,=20=CE=B8,=20shift,=20=CE=94,=20etc.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Opus 4.8 --- doc/source/gr_ore_poly.rst | 93 +++++++- src/gr_ore_poly.h | 46 ++++ src/gr_ore_poly/backshift_to_euler_univar.c | 69 ++++++ src/gr_ore_poly/convert.c | 82 +++++++ src/gr_ore_poly/ctx.c | 11 +- src/gr_ore_poly/ddx_to_euler.c | 69 ++++++ src/gr_ore_poly/differential_to_shift.c | 93 ++++++++ src/gr_ore_poly/euler_to_backshift_univar.c | 62 +++++ src/gr_ore_poly/euler_to_ddx.c | 67 ++++++ src/gr_ore_poly/shift_convert.c | 180 +++++++++++++++ src/gr_ore_poly/shift_to_differential.c | 85 +++++++ src/gr_ore_poly/test/main.c | 12 +- src/gr_ore_poly/test/t-convert.c | 168 ++++++++++++++ src/gr_ore_poly/test/t-ddx_to_euler.c | 145 ++++++++++++ src/gr_ore_poly/test/t-differential_shift.c | 242 ++++++++++++++++++++ src/gr_ore_poly/test/t-euler_to_ddx.c | 142 ++++++++++++ src/gr_ore_poly/test/t-shift_convert.c | 154 +++++++++++++ 17 files changed, 1715 insertions(+), 5 deletions(-) create mode 100644 src/gr_ore_poly/backshift_to_euler_univar.c create mode 100644 src/gr_ore_poly/convert.c create mode 100644 src/gr_ore_poly/ddx_to_euler.c create mode 100644 src/gr_ore_poly/differential_to_shift.c create mode 100644 src/gr_ore_poly/euler_to_backshift_univar.c create mode 100644 src/gr_ore_poly/euler_to_ddx.c create mode 100644 src/gr_ore_poly/shift_convert.c create mode 100644 src/gr_ore_poly/shift_to_differential.c create mode 100644 src/gr_ore_poly/test/t-convert.c create mode 100644 src/gr_ore_poly/test/t-ddx_to_euler.c create mode 100644 src/gr_ore_poly/test/t-differential_shift.c create mode 100644 src/gr_ore_poly/test/t-euler_to_ddx.c create mode 100644 src/gr_ore_poly/test/t-shift_convert.c diff --git a/doc/source/gr_ore_poly.rst b/doc/source/gr_ore_poly.rst index 6eed64ff3a..a15a0a59ba 100644 --- a/doc/source/gr_ore_poly.rst +++ b/doc/source/gr_ore_poly.rst @@ -5,8 +5,7 @@ .. note:: - This module is under construction. Functionality is currently limited to - memory management, additive arithmetic, and multiplication. + This module is under construction. A :type:`gr_ore_poly_t` represents a univariate Ore polynomial `L \in R[D]` implemented as a dense array of coefficients in a generic ring *R*. @@ -281,6 +280,96 @@ Action etc.). Returns ``GR_UNABLE`` when no standard action exists (`ORE_ALGEBRA_COMMUTATIVE`, `ORE_ALGEBRA_CUSTOM`) or it cannot be computed. +Conversions +------------------------------------------------------------------------------- + +These functions rewrite an operator from one Ore algebra into another. The +underscore versions act on raw coefficient arrays over the *coefficient* ring +*ctx* (the base ring of the Ore ring, not the Ore ring itself); *x* denotes its +generator of index *var*, that is, entry *var* of :func:`gr_gens`. Unless stated +otherwise, *res* has the same length *len* as *op* and must not alias it. + +A conversion is an exact rewriting up to a unit factor. Since *x* and the +forward shift `S` need not act invertibly in the chosen representation, the +result may be the input premultiplied **on the left** by a power of *x* or `S`, +returned in an output parameter; left multiplication does not change the +solution space of the operator. Conversions that substitute into or +Taylor-shift the variable require a :macro:`GR_CTX_GR_POLY` base ring and return +``GR_UNABLE`` otherwise. + +The first pair rewrites between the derivative `D = d/dx` and the Euler +derivative `\theta = x D`. + +.. function:: int _gr_ore_poly_ddx_to_euler(gr_ptr res, gr_srcptr op, slong len, slong var, gr_ctx_t ctx) + int _gr_ore_poly_euler_to_ddx(gr_ptr res, gr_srcptr op, slong len, slong var, gr_ctx_t ctx) + + Rewrite *op* between `D` and `\theta`. As `x` may not be invertible, + :func:`_gr_ore_poly_ddx_to_euler` returns `x^{len-1} \cdot \text{op}`; the + reverse rewriting needs no factor of `x`. + +The next pair converts among the four shift/difference algebras: the forward and +backward shifts `S`, `S^{-1}` and the forward and backward differences `S-1`, +`1-S^{-1}`. + +.. function:: int _gr_ore_poly_shift_convert(gr_ptr res, slong * epow, gr_srcptr op, slong len, ore_algebra_t src_alg, ore_algebra_t dst_alg, slong var, gr_ctx_t ctx) + + Convert *op* from *src_alg* to *dst_alg*. Each generator is a degree-one + Laurent polynomial in `S`, so on output ``*epow`` `= p` and *res* satisfy + `\text{res}_{dst} = S^{p} \cdot \text{op}_{src}`. A nonzero `p` (hence a + polynomial base ring) is needed only when crossing between the forward side + `S`, `S-1` and the backward side `S^{-1}`, `1-S^{-1}`, where `p = \mp(len-1)`. + Converting between the two difference types uses a single binomial transform. + The variable is the base ring generator of index *var*; only ``var == 0`` + over a univariate polynomial ring is implemented, and other indices return + ``GR_UNABLE``. Returns ``GR_DOMAIN`` if either algebra is not a + shift/difference type. + +The remaining functions bridge differential and shift operators through the +generating-series isomorphism `\theta \mapsto n`, `x \mapsto S^{-1}`: for +`f = \sum_n a_n x^n`, the Euler derivative `\theta = x d/dx` acts on `(a_n)` as +multiplication by `n` and `x` as the backward shift `S^{-1}`. All require a +`GR_CTX_GR_POLY` base ring and exchange the operator order with the coefficient +degree. + +.. function:: int _gr_ore_poly_euler_to_backshift_univar(gr_ptr res, slong reslen, gr_srcptr op, slong len, gr_ctx_t ctx) + int _gr_ore_poly_backshift_to_euler_univar(gr_ptr res, slong reslen, gr_srcptr op, slong len, gr_ctx_t ctx) + + The two inverse rewritings of the isomorphism between the Euler operator and + the backward shift `S^{-1}`. The caller allocates *res* to *reslen*, one more + than the largest coefficient degree of *op*. + +.. function:: int gr_ore_poly_differential_to_shift(gr_ore_poly_t res, slong * power, const gr_ore_poly_t op, gr_ore_poly_ctx_t res_ctx, gr_ore_poly_ctx_t op_ctx) + int gr_ore_poly_shift_to_differential(gr_ore_poly_t res, slong * power, const gr_ore_poly_t op, gr_ore_poly_ctx_t res_ctx, gr_ore_poly_ctx_t op_ctx) + + Convert a whole operator between a differential algebra + (`ORE_ALGEBRA_DERIVATIVE` or `ORE_ALGEBRA_EULER_DERIVATIVE`) and a + shift/difference algebra, routing through the Euler and backward-shift pivots + above. On output ``*power`` `= p` is the left power absorbed (`S^p` for + :func:`gr_ore_poly_differential_to_shift`, `x^p` for + :func:`gr_ore_poly_shift_to_differential`), so `\text{op}(f) = 0` and + `\text{res} = 0` have the same solutions up to the index shift by `p`. Both + contexts must share a `GR_CTX_GR_POLY` base ring (else ``GR_UNABLE``); returns + ``GR_DOMAIN`` if a context algebra is outside its family. + +.. function:: int gr_ore_poly_convert(gr_ore_poly_t res, slong * power, const gr_ore_poly_t op, gr_ore_poly_ctx_t res_ctx, gr_ore_poly_ctx_t op_ctx) + + Convert *op* from the algebra of *op_ctx* to that of *res_ctx*, provided both + belong to the same family: the differential family + (`ORE_ALGEBRA_DERIVATIVE`, `ORE_ALGEBRA_EULER_DERIVATIVE`) or the + shift/difference family (the four shift and difference types). This + dispatches to the appropriate conversion above and handles memory + management. On output ``*power`` `= p` is the left power absorbed, **but its + meaning depends on the pair of algebras**: for a conversion within the + differential family `\text{res} = x^{p} \cdot \text{op}` (a power of the base + generator `x`), whereas for one within the shift/difference family + `\text{res} = S^{p} \cdot \text{op}` (a power of the forward shift `S`). A + conversion to the same algebra is the identity, with `p = 0`. Returns + ``GR_UNABLE`` for an unsupported pair: one that crosses the + differential/shift boundary (use :func:`gr_ore_poly_differential_to_shift` or + :func:`gr_ore_poly_shift_to_differential` for those) or that involves an + algebra outside the two families, such as `ORE_ALGEBRA_Q_SHIFT`, + `ORE_ALGEBRA_MAHLER`, `ORE_ALGEBRA_FROBENIUS` or `ORE_ALGEBRA_COMMUTATIVE`. + Arithmetic ------------------------------------------------------------------------------- diff --git a/src/gr_ore_poly.h b/src/gr_ore_poly.h index 93d690a6f1..efeecc84ab 100644 --- a/src/gr_ore_poly.h +++ b/src/gr_ore_poly.h @@ -250,6 +250,52 @@ extern const gr_ore_poly_sigma_delta_t _gr_ore_poly_default_sigma_delta[]; WARN_UNUSED_RESULT int gr_ore_poly_apply_custom(gr_ptr res, const gr_ore_poly_t P, gr_srcptr f, gr_srcptr d1, gr_ore_poly_ctx_t ctx); WARN_UNUSED_RESULT int gr_ore_poly_apply(gr_ptr res, const gr_ore_poly_t P, gr_srcptr f, gr_ore_poly_ctx_t ctx); +/* Conversions between operator types. + + These act on raw coefficient arrays and take the coefficient ring context + (not the Ore polynomial ring context). The variable x is the generator of the + coefficient ring of index var, i.e. entry var of gr_gens(ctx). Since x need + not be invertible, ddx_to_euler multiplies the operator by x^(len-1); + euler_to_ddx leaves it unchanged. The result res has the same length len as + op. */ +WARN_UNUSED_RESULT int _gr_ore_poly_ddx_to_euler(gr_ptr res, gr_srcptr op, slong len, slong var, gr_ctx_t ctx); +WARN_UNUSED_RESULT int _gr_ore_poly_euler_to_ddx(gr_ptr res, gr_srcptr op, slong len, slong var, gr_ctx_t ctx); + +/* Convert a length-len operator between the shift/difference algebra types + src_alg and dst_alg. Each generator is a degree-one Laurent polynomial in the + forward shift S, so on output *p and res satisfy + S^p * res_in_dst = op_in_src (a power of S appears only when crossing between + the forward side S, S-1 and the backward side S^-1, 1-S^-1; it then needs a + GR_CTX_GR_POLY base ring). Returns GR_DOMAIN if either algebra is not a + shift/difference type. The variable is the base ring generator of index var; + only var == 0 over a univariate polynomial ring is implemented (else + GR_UNABLE). */ +WARN_UNUSED_RESULT int _gr_ore_poly_shift_convert(gr_ptr res, slong * p, gr_srcptr op, slong len, ore_algebra_t src_alg, ore_algebra_t dst_alg, slong var, gr_ctx_t ctx); + +/* Exact bridge of the isomorphism theta <-> n, x <-> backward shift S^{-1}, on + raw coefficient vectors over a GR_CTX_GR_POLY base ring (else GR_UNABLE). The + caller allocates res to reslen = 1 + the maximum coefficient degree of op. */ +WARN_UNUSED_RESULT int _gr_ore_poly_euler_to_backshift_univar(gr_ptr res, slong reslen, gr_srcptr op, slong len, gr_ctx_t ctx); +WARN_UNUSED_RESULT int _gr_ore_poly_backshift_to_euler_univar(gr_ptr res, slong reslen, gr_srcptr op, slong len, gr_ctx_t ctx); + +/* Convert between differential operators (op_ctx/res_ctx in {DERIVATIVE, + EULER_DERIVATIVE}) and shift/difference operators (in the four shift types), + via the generating-series isomorphism theta <-> n, x <-> backward shift. + Both contexts must share a GR_CTX_GR_POLY base ring. On output *p is such that + C^p * res corresponds to op under the isomorphism, where C is the companion + generator on the res side (S for *_to_shift, x for *_to_differential). + Returns GR_DOMAIN if a context algebra is outside its family. */ +WARN_UNUSED_RESULT int gr_ore_poly_differential_to_shift(gr_ore_poly_t res, slong * p, const gr_ore_poly_t op, gr_ore_poly_ctx_t res_ctx, gr_ore_poly_ctx_t op_ctx); +WARN_UNUSED_RESULT int gr_ore_poly_shift_to_differential(gr_ore_poly_t res, slong * p, const gr_ore_poly_t op, gr_ore_poly_ctx_t res_ctx, gr_ore_poly_ctx_t op_ctx); + +/* Convert an operator between two algebras of the same family (both differential + or both shift/difference), dispatching to the conversions above. On output *p + satisfies C^p * res = op, where the companion generator C depends on the pair + (x for differential conversions, the forward shift S for shift ones). + Returns GR_UNABLE for an unsupported pair (crossing the differential/shift + boundary, or involving an algebra outside those two families). */ +WARN_UNUSED_RESULT int gr_ore_poly_convert(gr_ore_poly_t res, slong * p, const gr_ore_poly_t op, gr_ore_poly_ctx_t res_ctx, gr_ore_poly_ctx_t op_ctx); + /* Arithmetic */ WARN_UNUSED_RESULT int gr_ore_poly_neg(gr_ore_poly_t res, const gr_ore_poly_t src, gr_ore_poly_ctx_t ctx); diff --git a/src/gr_ore_poly/backshift_to_euler_univar.c b/src/gr_ore_poly/backshift_to_euler_univar.c new file mode 100644 index 0000000000..8f6ff89cc9 --- /dev/null +++ b/src/gr_ore_poly/backshift_to_euler_univar.c @@ -0,0 +1,69 @@ +/* + 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 . +*/ + +/* generated using Claude Opus 4.8 */ + +#include "gr.h" +#include "gr_poly.h" +#include "gr_ore_poly.h" + +/* Bridge from the backward shift operator (K[n]) to the Euler operator + (K[x]), the inverse of _gr_ore_poly_euler_to_backshift: Taylor-shift + column k by +k, then transpose the (S^{-1}-degree, n-degree) coefficient + matrix. The caller must allocate res to reslen, the number of + theta-coefficients (= 1 + max n-degree of op). */ +int +_gr_ore_poly_backshift_to_euler_univar(gr_ptr res, slong reslen, gr_srcptr op, slong len, gr_ctx_t ctx) +{ + int status = GR_SUCCESS; + gr_ctx_struct * sctx; + slong bsz = ctx->sizeof_elem, ssz; + slong i, k; + gr_ptr c, cols; + + if (ctx->which_ring != GR_CTX_GR_POLY) + return GR_UNABLE; + + sctx = POLYNOMIAL_ELEM_CTX(ctx); + ssz = sctx->sizeof_elem; + + GR_TMP_INIT(c, sctx); + GR_TMP_INIT_VEC(cols, len, ctx); + for (k = 0; k < len; k++) + { + const gr_poly_struct * ok = (const gr_poly_struct *) GR_ENTRY(op, k, bsz); + gr_poly_struct * colk = (gr_poly_struct *) GR_ENTRY(cols, k, bsz); + status |= gr_poly_set(colk, ok, sctx); + status |= gr_set_si(c, k, sctx); + status |= gr_poly_taylor_shift(colk, colk, c, sctx); + } + + for (i = 0; i < reslen; i++) + { + gr_poly_struct * ri = (gr_poly_struct *) GR_ENTRY(res, i, bsz); + + gr_poly_fit_length(ri, len, sctx); + for (k = 0; k < len; k++) + { + gr_poly_struct * colk = (gr_poly_struct *) GR_ENTRY(cols, k, bsz); + gr_ptr dst = GR_ENTRY(ri->coeffs, k, ssz); + if (i < colk->length) + status |= gr_set(dst, GR_ENTRY(colk->coeffs, i, ssz), sctx); + else + status |= gr_zero(dst, sctx); + } + _gr_poly_set_length(ri, len, sctx); + _gr_poly_normalise(ri, sctx); + } + + GR_TMP_CLEAR_VEC(cols, len, ctx); + GR_TMP_CLEAR(c, sctx); + + return status; +} diff --git a/src/gr_ore_poly/convert.c b/src/gr_ore_poly/convert.c new file mode 100644 index 0000000000..03d96a1fcd --- /dev/null +++ b/src/gr_ore_poly/convert.c @@ -0,0 +1,82 @@ +/* + 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 . +*/ + +/* generated using Claude Opus 4.8 */ + +#include "gr.h" +#include "gr_vec.h" +#include "gr_ore_poly.h" + +int +gr_ore_poly_convert(gr_ore_poly_t res, slong * p, const gr_ore_poly_t op, + gr_ore_poly_ctx_t res_ctx, gr_ore_poly_ctx_t op_ctx) +{ + gr_ctx_struct * base = GR_ORE_POLY_ELEM_CTX(op_ctx); + ore_algebra_t sa = GR_ORE_POLY_CTX(op_ctx)->which_algebra; + ore_algebra_t da = GR_ORE_POLY_CTX(res_ctx)->which_algebra; + int status = GR_SUCCESS; + slong len, sp = 0; + + *p = 0; + + if (GR_ORE_POLY_ELEM_CTX(res_ctx) != base) + return GR_UNABLE; + + len = op->length; + if (len == 0) + return gr_ore_poly_zero(res, res_ctx); + + gr_ore_poly_fit_length(res, len, res_ctx); + + if ((sa == ORE_ALGEBRA_DERIVATIVE || sa == ORE_ALGEBRA_EULER_DERIVATIVE) + && (da == ORE_ALGEBRA_DERIVATIVE || da == ORE_ALGEBRA_EULER_DERIVATIVE)) + { + /* differential family {d/dx, theta} */ + slong var = GR_ORE_POLY_ORE_DATA(op_ctx)->base_var; + + if (sa == da) + status |= _gr_vec_set(res->coeffs, op->coeffs, len, base); + else if (sa == ORE_ALGEBRA_DERIVATIVE) /* d/dx -> theta */ + { + status |= _gr_ore_poly_ddx_to_euler(res->coeffs, op->coeffs, len, var, base); + /* res = x^(len-1) * op, i.e. x^p * res = op with p = -(len-1) */ + *p = -(len - 1); + } + else /* theta -> d/dx */ + status |= _gr_ore_poly_euler_to_ddx(res->coeffs, op->coeffs, len, var, base); + } + else + { + /* _gr_ore_poly_shift_convert handles the shift/difference family (the + same-algebra case included) and returns GR_DOMAIN for anything that is + not a shift/difference type, which is exactly the remaining unsupported + case: a boundary-crossing pair, or an algebra such as Q_SHIFT, MAHLER + or COMMUTATIVE (whose identity we deliberately do not claim, since its + parameters could differ between the two contexts). */ + /* a custom algebra has no ore_data; shift_convert then returns GR_DOMAIN + regardless of var, so the placeholder 0 is harmless */ + gr_ore_poly_ore_data_t * od = GR_ORE_POLY_ORE_DATA(op_ctx); + slong var = (od != NULL) ? od->base_var : 0; + + status = _gr_ore_poly_shift_convert(res->coeffs, &sp, op->coeffs, len, + sa, da, var, base); + if (status == GR_DOMAIN) + status = GR_UNABLE; + else + *p = sp; + } + + if (status == GR_SUCCESS) + { + _gr_ore_poly_set_length(res, len, res_ctx); + _gr_ore_poly_normalise(res, res_ctx); + } + + return status; +} diff --git a/src/gr_ore_poly/ctx.c b/src/gr_ore_poly/ctx.c index fcdf4b40a9..3dc997aa9c 100644 --- a/src/gr_ore_poly/ctx.c +++ b/src/gr_ore_poly/ctx.c @@ -409,11 +409,18 @@ gr_ore_poly_ctx_init_randtest(gr_ore_poly_ctx_t ctx, flint_rand_t state, gr_ctx_t base_ring) { ore_algebra_t alg = ore_algebra_randtest(state); - /* todo: base_var = n_randint(gr_ctx_ngens(base_ring)) */ - slong base_var = 0; + slong ngens = 0; + slong base_var; int status = GR_SUCCESS; + if (gr_ctx_ngens(&ngens, base_ring) != GR_SUCCESS) + ngens = 0; + /* an invalid index when the base ring has no generator, so that a spurious + reliance on a default generator surfaces as GR_UNABLE rather than silently + using generator 0 */ + base_var = (ngens > 0) ? (slong) n_randint(state, ngens) : -1; + if (alg == ORE_ALGEBRA_CUSTOM) { /* sigma_delta_commutative does not use ore_data */ diff --git a/src/gr_ore_poly/ddx_to_euler.c b/src/gr_ore_poly/ddx_to_euler.c new file mode 100644 index 0000000000..f6f8cf61b3 --- /dev/null +++ b/src/gr_ore_poly/ddx_to_euler.c @@ -0,0 +1,69 @@ +/* + 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 . +*/ + +/* generated using Claude Opus 4.8 */ + +#include "fmpz.h" +#include "gr.h" +#include "gr_mat.h" +#include "gr_vec.h" +#include "gr_ore_poly.h" + +/* x^(len-1) * sum_i op[i] D^i, expanded via x^i D^i = sum_k s(i,k) theta^k with + s(i,k) the signed Stirling numbers of the first kind. Here x is the generator + of index var. */ +int +_gr_ore_poly_ddx_to_euler(gr_ptr res, gr_srcptr op, slong len, slong var, gr_ctx_t ctx) +{ + int status = GR_SUCCESS; + slong sz = ctx->sizeof_elem; + slong i, k; + gr_ctx_t zz; + gr_mat_t S; + gr_vec_t gens; + gr_ptr x, xp, c; + + if (len <= 0) + return GR_SUCCESS; + + gr_ctx_init_fmpz(zz); + gr_mat_init(S, len, len, zz); + /* s(i,k) = 0 for k > i, so only the lower triangle of S is used below */ + status |= gr_mat_stirling(S, 1, zz); + + GR_TMP_INIT3(x, xp, c, ctx); + + status |= _gr_vec_zero(res, len, ctx); + + gr_vec_init(gens, 0, ctx); + status |= gr_gens(gens, ctx); + if (var >= 0 && var < gens->length) + status |= gr_set(x, gr_vec_entry_srcptr(gens, var, ctx), ctx); + else + status |= GR_UNABLE; + gr_vec_clear(gens, ctx); + + status |= gr_one(xp, ctx); + + for (i = len - 1; i >= 0; i--) + { + status |= gr_mul(c, xp, GR_ENTRY(op, i, sz), ctx); + for (k = 0; k <= i; k++) + status |= gr_addmul_fmpz(GR_ENTRY(res, k, sz), c, + gr_mat_entry_srcptr(S, i, k, zz), ctx); + if (i > 0) + status |= gr_mul(xp, xp, x, ctx); + } + + GR_TMP_CLEAR3(x, xp, c, ctx); + gr_mat_clear(S, zz); + gr_ctx_clear(zz); + + return status; +} diff --git a/src/gr_ore_poly/differential_to_shift.c b/src/gr_ore_poly/differential_to_shift.c new file mode 100644 index 0000000000..ff736a776e --- /dev/null +++ b/src/gr_ore_poly/differential_to_shift.c @@ -0,0 +1,93 @@ +/* + 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 . +*/ + +/* generated using Claude Opus 4.8 */ + +#include "gr.h" +#include "gr_poly.h" +#include "gr_vec.h" +#include "gr_ore_poly.h" + +int +gr_ore_poly_differential_to_shift(gr_ore_poly_t res, slong * p, const gr_ore_poly_t op, + gr_ore_poly_ctx_t res_ctx, gr_ore_poly_ctx_t op_ctx) +{ + gr_ctx_struct * base = GR_ORE_POLY_ELEM_CTX(op_ctx); + ore_algebra_t sa = GR_ORE_POLY_CTX(op_ctx)->which_algebra; + ore_algebra_t da = GR_ORE_POLY_CTX(res_ctx)->which_algebra; + slong var = GR_ORE_POLY_ORE_DATA(op_ctx)->base_var; + slong bsz = base->sizeof_elem; + int status = GR_SUCCESS; + slong len, rlen = 0, a = 0, sp = 0, i; + gr_srcptr eul; + gr_ptr cur = NULL, rec = NULL; + + *p = 0; + + if (GR_ORE_POLY_ELEM_CTX(res_ctx) != base || base->which_ring != GR_CTX_GR_POLY) + return GR_UNABLE; + + len = op->length; + if (len == 0) + return gr_ore_poly_zero(res, res_ctx); + + /* differential -> Euler (the Euler operator is op itself when sa is already + the Euler derivative) */ + if (sa == ORE_ALGEBRA_DERIVATIVE) + { + GR_TMP_INIT_VEC(cur, len, base); + status |= _gr_ore_poly_ddx_to_euler(cur, op->coeffs, len, var, base); + eul = cur; + a = len - 1; + } + else if (sa == ORE_ALGEBRA_EULER_DERIVATIVE) + { + eul = op->coeffs; + a = 0; + } + else + { + status = GR_DOMAIN; + goto cleanup; + } + + for (i = 0; i < len; i++) + { + const gr_poly_struct * q = (const gr_poly_struct *) GR_ENTRY(eul, i, bsz); + if (q->length > rlen) + rlen = q->length; + } + if (rlen == 0) + { + status |= gr_ore_poly_zero(res, res_ctx); + goto cleanup; + } + + /* Euler -> backward shift, then backward shift -> destination shift type */ + GR_TMP_INIT_VEC(rec, rlen, base); + status |= _gr_ore_poly_euler_to_backshift_univar(rec, rlen, eul, len, base); + + gr_ore_poly_fit_length(res, rlen, res_ctx); + status |= _gr_ore_poly_shift_convert(res->coeffs, &sp, rec, rlen, + ORE_ALGEBRA_BACKWARD_SHIFT, da, var, base); + if (status == GR_SUCCESS) + { + _gr_ore_poly_set_length(res, rlen, res_ctx); + _gr_ore_poly_normalise(res, res_ctx); + *p = sp + a; + } + +cleanup: + if (cur != NULL) + GR_TMP_CLEAR_VEC(cur, len, base); + if (rec != NULL) + GR_TMP_CLEAR_VEC(rec, rlen, base); + + return status; +} diff --git a/src/gr_ore_poly/euler_to_backshift_univar.c b/src/gr_ore_poly/euler_to_backshift_univar.c new file mode 100644 index 0000000000..555b8c04aa --- /dev/null +++ b/src/gr_ore_poly/euler_to_backshift_univar.c @@ -0,0 +1,62 @@ +/* + 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 . +*/ + +/* generated using Claude Opus 4.8 */ + +#include "gr.h" +#include "gr_poly.h" +#include "gr_ore_poly.h" + +/* Bridge from the Euler operator (K[x], theta = x d/dx) to the backward + shift operator (K[n], S^{-1} a(n) = a(n-1)) under the exact + isomorphism theta <-> n, x <-> S^{-1}. Writing op = sum_k x^k Q_k(theta), the + image is sum_k Q_k(n-k) S^{-k}: transpose the (theta-degree, x-degree) + coefficient matrix and Taylor-shift column k by -k. The caller must allocate + res to reslen, the number of S^{-1}-coefficients (= 1 + max x-degree of + op). */ +int +_gr_ore_poly_euler_to_backshift_univar(gr_ptr res, slong reslen, gr_srcptr op, slong len, gr_ctx_t ctx) +{ + int status = GR_SUCCESS; + gr_ctx_struct * sctx; + slong bsz = ctx->sizeof_elem, ssz; + slong i, k; + gr_ptr c; + + if (ctx->which_ring != GR_CTX_GR_POLY) + return GR_UNABLE; + + sctx = POLYNOMIAL_ELEM_CTX(ctx); + ssz = sctx->sizeof_elem; + + GR_TMP_INIT(c, sctx); + for (k = 0; k < reslen; k++) + { + gr_poly_struct * rk = (gr_poly_struct *) GR_ENTRY(res, k, bsz); + + gr_poly_fit_length(rk, len, sctx); + for (i = 0; i < len; i++) + { + const gr_poly_struct * oi = (const gr_poly_struct *) GR_ENTRY(op, i, bsz); + gr_ptr dst = GR_ENTRY(rk->coeffs, i, ssz); + if (k < oi->length) + status |= gr_set(dst, GR_ENTRY(oi->coeffs, k, ssz), sctx); + else + status |= gr_zero(dst, sctx); + } + _gr_poly_set_length(rk, len, sctx); + _gr_poly_normalise(rk, sctx); + + status |= gr_set_si(c, -k, sctx); + status |= gr_poly_taylor_shift(rk, rk, c, sctx); + } + GR_TMP_CLEAR(c, sctx); + + return status; +} diff --git a/src/gr_ore_poly/euler_to_ddx.c b/src/gr_ore_poly/euler_to_ddx.c new file mode 100644 index 0000000000..f7e21581a5 --- /dev/null +++ b/src/gr_ore_poly/euler_to_ddx.c @@ -0,0 +1,67 @@ +/* + 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 . +*/ + +/* generated using Claude Opus 4.8 */ + +#include "fmpz.h" +#include "gr.h" +#include "gr_mat.h" +#include "gr_vec.h" +#include "gr_ore_poly.h" + +/* sum_k op[k] theta^k expanded via theta^k = sum_i S(k,i) x^i D^i with S(k,i) + the Stirling numbers of the second kind. Here x is the generator of index + var. */ +int +_gr_ore_poly_euler_to_ddx(gr_ptr res, gr_srcptr op, slong len, slong var, gr_ctx_t ctx) +{ + int status = GR_SUCCESS; + slong sz = ctx->sizeof_elem; + slong i, k; + gr_ctx_t zz; + gr_mat_t S; + gr_vec_t gens; + gr_ptr x, xp, t; + + if (len <= 0) + return GR_SUCCESS; + + gr_ctx_init_fmpz(zz); + gr_mat_init(S, len, len, zz); + /* S(k,i) = 0 for i > k, so only the lower triangle of S is used below */ + status |= gr_mat_stirling(S, 2, zz); + + GR_TMP_INIT3(x, xp, t, ctx); + + gr_vec_init(gens, 0, ctx); + status |= gr_gens(gens, ctx); + if (var >= 0 && var < gens->length) + status |= gr_set(x, gr_vec_entry_srcptr(gens, var, ctx), ctx); + else + status |= GR_UNABLE; + gr_vec_clear(gens, ctx); + + status |= gr_one(xp, ctx); + + for (i = 0; i < len; i++) + { + status |= gr_zero(t, ctx); + for (k = i; k < len; k++) + status |= gr_addmul_fmpz(t, GR_ENTRY(op, k, sz), + gr_mat_entry_srcptr(S, k, i, zz), ctx); + status |= gr_mul(GR_ENTRY(res, i, sz), xp, t, ctx); + status |= gr_mul(xp, xp, x, ctx); + } + + GR_TMP_CLEAR3(x, xp, t, ctx); + gr_mat_clear(S, zz); + gr_ctx_clear(zz); + + return status; +} diff --git a/src/gr_ore_poly/shift_convert.c b/src/gr_ore_poly/shift_convert.c new file mode 100644 index 0000000000..ed2f3032c3 --- /dev/null +++ b/src/gr_ore_poly/shift_convert.c @@ -0,0 +1,180 @@ +/* + 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 . +*/ + +/* generated using Claude Opus 4.8 */ + +#include "fmpz.h" +#include "gr.h" +#include "gr_poly.h" +#include "gr_vec.h" +#include "gr_ore_poly.h" + +static int +shift_alg_is_shift(ore_algebra_t a) +{ + return a == ORE_ALGEBRA_FORWARD_SHIFT || a == ORE_ALGEBRA_BACKWARD_SHIFT + || a == ORE_ALGEBRA_FORWARD_DIFFERENCE || a == ORE_ALGEBRA_BACKWARD_DIFFERENCE; +} + +static int +shift_alg_is_backward(ore_algebra_t a) +{ + return a == ORE_ALGEBRA_BACKWARD_SHIFT || a == ORE_ALGEBRA_BACKWARD_DIFFERENCE; +} + +static int +shift_alg_is_difference(ore_algebra_t a) +{ + return a == ORE_ALGEBRA_FORWARD_DIFFERENCE || a == ORE_ALGEBRA_BACKWARD_DIFFERENCE; +} + +/* In-place change of basis between a shift generator and the associated + finite-difference generator, computing + vec[m] <- sum_{k >= m} w(k, m) * binomial(k, m) * vec[k], + where the sign w is 1 for sign == 0, (-1)^(k-m) for sign == +1, and (-1)^m + for sign == -1. The binomial coefficients along each row are built up + incrementally, as in _arb_poly_binomial_transform_basecase. The update reads + only entries of index >= m, so writing vec[m] in place is safe. */ +static int +binomial_transform(gr_ptr vec, slong len, int sign, gr_ctx_t ctx) +{ + int status = GR_SUCCESS; + slong sz = ctx->sizeof_elem, m, k; + gr_ptr acc; + fmpz_t t; + + GR_TMP_INIT(acc, ctx); + fmpz_init(t); + + for (m = 0; m < len; m++) + { + status |= gr_zero(acc, ctx); + for (k = m; k < len; k++) + { + if (k == m) + fmpz_set_si(t, (sign == -1 && (m & 1)) ? -1 : 1); + else + { + fmpz_mul_si(t, t, (sign == 1) ? -k : k); + fmpz_divexact_ui(t, t, (ulong) (k - m)); + } + status |= gr_addmul_fmpz(acc, GR_ENTRY(vec, k, sz), t, ctx); + } + status |= gr_set(GR_ENTRY(vec, m, sz), acc, ctx); + } + + fmpz_clear(t); + GR_TMP_CLEAR(acc, ctx); + + return status; +} + +static void +reverse(gr_ptr vec, slong len, gr_ctx_t ctx) +{ + slong sz = ctx->sizeof_elem, k; + + for (k = 0; k < len / 2; k++) + gr_swap(GR_ENTRY(vec, k, sz), GR_ENTRY(vec, len - 1 - k, sz), ctx); +} + +/* Taylor-shift every coefficient by p; requires a polynomial base ring. */ +static int +taylor_shift_all(gr_ptr vec, slong len, slong p, gr_ctx_t ctx) +{ + int status = GR_SUCCESS; + slong sz = ctx->sizeof_elem, k; + gr_ctx_struct * sctx = POLYNOMIAL_ELEM_CTX(ctx); + gr_ptr c; + + GR_TMP_INIT(c, sctx); + status |= gr_set_si(c, p, sctx); + for (k = 0; k < len; k++) + { + gr_poly_struct * vk = (gr_poly_struct *) GR_ENTRY(vec, k, sz); + status |= gr_poly_taylor_shift(vk, vk, c, sctx); + } + GR_TMP_CLEAR(c, sctx); + + return status; +} + +int +_gr_ore_poly_shift_convert(gr_ptr res, slong * p, gr_srcptr op, slong len, + ore_algebra_t src_alg, ore_algebra_t dst_alg, slong var, gr_ctx_t ctx) +{ + int status = GR_SUCCESS; + int src_back, dst_back, crossing; + slong s; + + *p = 0; + + if (!shift_alg_is_shift(src_alg) || !shift_alg_is_shift(dst_alg)) + return GR_DOMAIN; + + if (len <= 0) + return GR_SUCCESS; + + /* The implementation is limited to a univariate polynomial base ring, whose + only generator has index 0. */ + if (var != 0) + return GR_UNABLE; + + if (src_alg == dst_alg) + return _gr_vec_set(res, op, len, ctx); + + src_back = shift_alg_is_backward(src_alg); + dst_back = shift_alg_is_backward(dst_alg); + crossing = (src_back != dst_back); + /* res is computed as S^s * op, so the reported p (satisfying S^p * res = op) + is -s. */ + s = dst_back ? -(len - 1) : (len - 1); + + /* Crossing between the forward side S, S - 1 and the backward side S^{-1}, + 1 - S^{-1} conjugates by a power of S, realized as a Taylor shift of the + coefficients; this needs a polynomial base ring once len >= 2. */ + if (crossing && len >= 2 && ctx->which_ring != GR_CTX_GR_POLY) + return GR_UNABLE; + + status |= _gr_vec_set(res, op, len, ctx); + + /* The two difference generators sit on opposite sides, so converting between + them always crosses. Because S - 1 = S (1 - S^{-1}), a single binomial + transform suffices: reversing the coefficients turns the combination into + the standard one handled by binomial_transform. */ + if (shift_alg_is_difference(src_alg) && shift_alg_is_difference(dst_alg)) + { + reverse(res, len, ctx); + status |= binomial_transform(res, len, dst_back ? 1 : 0, ctx); + reverse(res, len, ctx); + if (s != 0) + status |= taylor_shift_all(res, len, s, ctx); + *p = -s; + return status; + } + + /* General pipeline: source difference -> its shift, cross sides if needed, + then shift -> destination difference. At most one binomial transform runs + (the difference <-> difference case above is the only one needing two). */ + if (shift_alg_is_difference(src_alg)) + status |= binomial_transform(res, len, src_back ? -1 : 1, ctx); + + if (crossing) + { + if (s != 0) + status |= taylor_shift_all(res, len, s, ctx); + reverse(res, len, ctx); + *p = -s; + } + + if (shift_alg_is_difference(dst_alg)) + status |= binomial_transform(res, len, dst_back ? -1 : 0, ctx); + + return status; +} diff --git a/src/gr_ore_poly/shift_to_differential.c b/src/gr_ore_poly/shift_to_differential.c new file mode 100644 index 0000000000..35413ff6cc --- /dev/null +++ b/src/gr_ore_poly/shift_to_differential.c @@ -0,0 +1,85 @@ +/* + 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 . +*/ + +/* generated using Claude Opus 4.8 */ + +#include "gr.h" +#include "gr_poly.h" +#include "gr_vec.h" +#include "gr_ore_poly.h" + +int +gr_ore_poly_shift_to_differential(gr_ore_poly_t res, slong * p, const gr_ore_poly_t op, + gr_ore_poly_ctx_t res_ctx, gr_ore_poly_ctx_t op_ctx) +{ + gr_ctx_struct * base = GR_ORE_POLY_ELEM_CTX(op_ctx); + ore_algebra_t sa = GR_ORE_POLY_CTX(op_ctx)->which_algebra; + ore_algebra_t da = GR_ORE_POLY_CTX(res_ctx)->which_algebra; + slong var = GR_ORE_POLY_ORE_DATA(res_ctx)->base_var; + slong bsz = base->sizeof_elem; + int status = GR_SUCCESS; + slong len, elen = 0, sp = 0, k; + gr_ptr rec = NULL, eul = NULL; + + *p = 0; + + if (GR_ORE_POLY_ELEM_CTX(res_ctx) != base || base->which_ring != GR_CTX_GR_POLY) + return GR_UNABLE; + + if (da != ORE_ALGEBRA_DERIVATIVE && da != ORE_ALGEBRA_EULER_DERIVATIVE) + return GR_DOMAIN; + + len = op->length; + if (len == 0) + return gr_ore_poly_zero(res, res_ctx); + + /* source shift type -> backward shift (the algebra is validated here) */ + GR_TMP_INIT_VEC(rec, len, base); + status |= _gr_ore_poly_shift_convert(rec, &sp, op->coeffs, len, + sa, ORE_ALGEBRA_BACKWARD_SHIFT, var, base); + if (status != GR_SUCCESS) + goto cleanup; + + for (k = 0; k < len; k++) + { + const gr_poly_struct * q = (const gr_poly_struct *) GR_ENTRY(rec, k, bsz); + if (q->length > elen) + elen = q->length; + } + if (elen == 0) + { + status |= gr_ore_poly_zero(res, res_ctx); + goto cleanup; + } + + /* backward shift -> Euler -> destination differential type. For the Euler + destination, backshift_to_euler_univar writes the result in place. */ + gr_ore_poly_fit_length(res, elen, res_ctx); + if (da == ORE_ALGEBRA_EULER_DERIVATIVE) + { + status |= _gr_ore_poly_backshift_to_euler_univar(res->coeffs, elen, rec, len, base); + } + else + { + GR_TMP_INIT_VEC(eul, elen, base); + status |= _gr_ore_poly_backshift_to_euler_univar(eul, elen, rec, len, base); + status |= _gr_ore_poly_euler_to_ddx(res->coeffs, eul, elen, var, base); + } + _gr_ore_poly_set_length(res, elen, res_ctx); + _gr_ore_poly_normalise(res, res_ctx); + *p = -sp; + +cleanup: + if (rec != NULL) + GR_TMP_CLEAR_VEC(rec, len, base); + if (eul != NULL) + GR_TMP_CLEAR_VEC(eul, elen, base); + + return status; +} diff --git a/src/gr_ore_poly/test/main.c b/src/gr_ore_poly/test/main.c index 838e5889c1..edbf9648a3 100644 --- a/src/gr_ore_poly/test/main.c +++ b/src/gr_ore_poly/test/main.c @@ -17,6 +17,11 @@ #include "t-mul.c" #include "t-divrem.c" #include "t-apply.c" +#include "t-ddx_to_euler.c" +#include "t-euler_to_ddx.c" +#include "t-shift_convert.c" +#include "t-differential_shift.c" +#include "t-convert.c" /* Array of test functions ***************************************************/ @@ -27,7 +32,12 @@ test_struct tests[] = TEST_FUNCTION(gr_ore_poly_sigma_delta), TEST_FUNCTION(gr_ore_poly_mul), TEST_FUNCTION(gr_ore_poly_divrem), - TEST_FUNCTION(gr_ore_poly_apply) + TEST_FUNCTION(gr_ore_poly_apply), + TEST_FUNCTION(gr_ore_poly_ddx_to_euler), + TEST_FUNCTION(gr_ore_poly_euler_to_ddx), + TEST_FUNCTION(gr_ore_poly_shift_convert), + TEST_FUNCTION(gr_ore_poly_differential_to_shift), + TEST_FUNCTION(gr_ore_poly_convert) }; /* main function *************************************************************/ diff --git a/src/gr_ore_poly/test/t-convert.c b/src/gr_ore_poly/test/t-convert.c new file mode 100644 index 0000000000..282fe765bb --- /dev/null +++ b/src/gr_ore_poly/test/t-convert.c @@ -0,0 +1,168 @@ +/* + 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 . +*/ + +/* generated using Claude Opus 4.8 */ + +#include "test_helpers.h" +#include "gr.h" +#include "gr_poly.h" +#include "gr_vec.h" +#include "gr_ore_poly.h" + +TEST_GR_FUNCTION_START(gr_ore_poly_convert, state, count_success, count_domain, count_unable) +{ + for (slong iter = 0; iter < 1000 * flint_test_multiplier(); iter++) + { + gr_ctx_t cctx, cctx2; + gr_ore_poly_ctx_t ctx_src, ctx_dst; + ore_algebra_t sa, da; + slong power = 0, j; + int independent, cstatus, status = GR_SUCCESS; + + /* Most of the time convert between two algebras over a shared base ring; + occasionally use unrelated base rings to check that the function copes + sensibly (it returns GR_UNABLE) rather than misbehaving. */ + independent = (n_randint(state, 4) == 0); + + gr_ore_poly_ctx_init_randtest2(cctx, ctx_src, state); + if (independent) + gr_ore_poly_ctx_init_randtest2(cctx2, ctx_dst, state); + else + gr_ore_poly_ctx_init_randtest(ctx_dst, state, cctx); + + sa = GR_ORE_POLY_CTX(ctx_src)->which_algebra; + da = GR_ORE_POLY_CTX(ctx_dst)->which_algebra; + + gr_ore_poly_t op, res; + gr_ore_poly_init(op, ctx_src); + gr_ore_poly_init(res, ctx_dst); + + status |= gr_ore_poly_randtest(op, state, 1 + n_randint(state, 5), ctx_src); + + cstatus = gr_ore_poly_convert(res, &power, op, ctx_dst, ctx_src); + status |= cstatus; + + /* Over these representative base rings the conversion never fails, as + long as both algebras stay within one family (differential or + shift/difference) and the two contexts share the base ring. */ + int sa_diff = (sa == ORE_ALGEBRA_DERIVATIVE || sa == ORE_ALGEBRA_EULER_DERIVATIVE); + int da_diff = (da == ORE_ALGEBRA_DERIVATIVE || da == ORE_ALGEBRA_EULER_DERIVATIVE); + int sa_shift = (sa == ORE_ALGEBRA_FORWARD_SHIFT || sa == ORE_ALGEBRA_BACKWARD_SHIFT + || sa == ORE_ALGEBRA_FORWARD_DIFFERENCE || sa == ORE_ALGEBRA_BACKWARD_DIFFERENCE); + int da_shift = (da == ORE_ALGEBRA_FORWARD_SHIFT || da == ORE_ALGEBRA_BACKWARD_SHIFT + || da == ORE_ALGEBRA_FORWARD_DIFFERENCE || da == ORE_ALGEBRA_BACKWARD_DIFFERENCE); + int expect_success = (!independent && ((sa_diff && da_diff) || (sa_shift && da_shift)) + && cctx->which_ring == GR_CTX_GR_POLY + && (POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_FMPZ + || POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_CC_ACB + || POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_NMOD)); + if (expect_success && cstatus != GR_SUCCESS) + { + flint_printf("FAIL: convert unexpected failure\n"); + flint_printf("sa = %d, da = %d\n", sa, da); + flint_abort(); + } + + /* When the conversion succeeds, the converted operator must act on the + base ring consistently with the original, up to the reported left + power: x^power for the differential family, the forward shift S^power + (a Taylor shift) for the shift/difference family. */ + if (cstatus == GR_SUCCESS) + { + int differential = (sa == ORE_ALGEBRA_DERIVATIVE || sa == ORE_ALGEBRA_EULER_DERIVATIVE) + && (da == ORE_ALGEBRA_DERIVATIVE || da == ORE_ALGEBRA_EULER_DERIVATIVE); + + gr_ptr f = gr_heap_init(cctx); + gr_ptr g_src = gr_heap_init(cctx); + gr_ptr g_dst = gr_heap_init(cctx); + gr_ptr corrected = gr_heap_init(cctx); + + for (j = 0; j < 2; j++) + { + int s1, s2; + + status |= gr_randtest(f, state, cctx); + s1 = gr_ore_poly_apply(g_src, op, f, ctx_src); + s2 = gr_ore_poly_apply(g_dst, res, f, ctx_dst); + status |= s1 | s2; + + if (s1 != GR_SUCCESS || s2 != GR_SUCCESS) + continue; + + /* power == 0 is the identity correction and is the only case + reachable over a non-polynomial base ring (apply succeeds there + only for a constant operator, whose conversion has power 0). A + nonzero power implies a GR_CTX_GR_POLY base ring. */ + if (power == 0) + { + status |= gr_set(corrected, g_src, cctx); + } + else if (differential) + { + /* x^power * res = op, so corrected = x^{-power} * g_src, x the + generator of index var */ + slong var = GR_ORE_POLY_ORE_DATA(ctx_src)->base_var; + gr_vec_t gens; + gr_ptr xp = gr_heap_init(cctx); + + gr_vec_init(gens, 0, cctx); + status |= gr_gens(gens, cctx); + if (var < gens->length) + status |= gr_set(xp, gr_vec_entry_srcptr(gens, var, cctx), cctx); + else + status |= GR_UNABLE; + gr_vec_clear(gens, cctx); + + status |= gr_pow_ui(xp, xp, (ulong) (-power), cctx); + status |= gr_mul(corrected, xp, g_src, cctx); + + gr_heap_clear(xp, cctx); + } + else + { + /* S^power * res = op, so corrected = S^{-power}(g_src), a + Taylor shift of the base ring element (a polynomial) */ + gr_ctx_struct * sctx = POLYNOMIAL_ELEM_CTX(cctx); + gr_ptr c = gr_heap_init(sctx); + + status |= gr_set_si(c, -power, sctx); + status |= gr_poly_taylor_shift((gr_poly_struct *) corrected, + (gr_poly_struct *) g_src, c, sctx); + gr_heap_clear(c, sctx); + } + + if (status == GR_SUCCESS && gr_equal(g_dst, corrected, cctx) == T_FALSE) + { + flint_printf("FAIL: convert application inconsistent\n"); + flint_printf("sa = %d, da = %d, power = %wd\n", sa, da, power); + flint_abort(); + } + } + + gr_heap_clear(f, cctx); + gr_heap_clear(g_src, cctx); + gr_heap_clear(g_dst, cctx); + gr_heap_clear(corrected, cctx); + } + + count_success += (status == GR_SUCCESS); + count_domain += ((status & GR_DOMAIN) != 0); + count_unable += ((status & GR_UNABLE) != 0); + + gr_ore_poly_clear(op, ctx_src); + gr_ore_poly_clear(res, ctx_dst); + gr_ore_poly_ctx_clear(ctx_src); + gr_ore_poly_ctx_clear(ctx_dst); + gr_ctx_clear(cctx); + if (independent) + gr_ctx_clear(cctx2); + } + + TEST_GR_FUNCTION_END(state, count_success, count_domain, count_unable); +} diff --git a/src/gr_ore_poly/test/t-ddx_to_euler.c b/src/gr_ore_poly/test/t-ddx_to_euler.c new file mode 100644 index 0000000000..04b581f377 --- /dev/null +++ b/src/gr_ore_poly/test/t-ddx_to_euler.c @@ -0,0 +1,145 @@ +/* + 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 . +*/ + +/* generated using Claude Opus 4.8 */ + +#include "test_helpers.h" +#include "gr.h" +#include "gr_vec.h" +#include "gr_ore_poly.h" + +TEST_GR_FUNCTION_START(gr_ore_poly_ddx_to_euler, state, count_success, count_domain, count_unable) +{ + for (slong iter = 0; iter < 1000 * flint_test_multiplier(); iter++) + { + gr_ctx_t cctx; + gr_ore_poly_ctx_t ctx_d, ctx_e; + gr_vec_t gens; + slong len, i, j, sz, ngens, var; + int status = GR_SUCCESS; + + /* occasionally use a multivariate base ring so that the conversions are + exercised with a generator other than the default one */ + switch (n_randint(state, 8)) + { + case 0: + gr_ctx_init_random(cctx, state); + break; + case 1: + case 2: + gr_ctx_init_random_mpoly(cctx, state); + break; + default: + gr_ctx_init_random_poly(cctx, state); + break; + } + + sz = cctx->sizeof_elem; + + ngens = 0; + status |= gr_ctx_ngens(&ngens, cctx); + var = (ngens >= 1) ? (slong) n_randint(state, ngens) : 0; + + gr_ore_poly_ctx_init(ctx_d, cctx, var, ORE_ALGEBRA_DERIVATIVE); + gr_ore_poly_ctx_init(ctx_e, cctx, var, ORE_ALGEBRA_EULER_DERIVATIVE); + + gr_ore_poly_t P_d, P_e, P_d2, scaled_P; + gr_ore_poly_init(P_d, ctx_d); + gr_ore_poly_init(P_e, ctx_e); + gr_ore_poly_init(P_d2, ctx_d); + gr_ore_poly_init(scaled_P, ctx_d); + + status |= gr_ore_poly_randtest(P_d, state, 1 + n_randint(state, 5), ctx_d); + len = P_d->length; + + gr_ore_poly_fit_length(P_e, len, ctx_e); + status |= _gr_ore_poly_ddx_to_euler(P_e->coeffs, P_d->coeffs, len, var, cctx); + _gr_ore_poly_set_length(P_e, len, ctx_e); + _gr_ore_poly_normalise(P_e, ctx_e); + + gr_ore_poly_fit_length(P_d2, len, ctx_d); + status |= _gr_ore_poly_euler_to_ddx(P_d2->coeffs, P_e->coeffs, len, var, cctx); + _gr_ore_poly_set_length(P_d2, len, ctx_d); + _gr_ore_poly_normalise(P_d2, ctx_d); + + gr_ptr xpow = gr_heap_init(cctx); + gr_ptr f = gr_heap_init(cctx); + gr_ptr g_d = gr_heap_init(cctx); + gr_ptr g_e = gr_heap_init(cctx); + gr_ptr scaled = gr_heap_init(cctx); + + /* x is the generator of index var */ + gr_vec_init(gens, 0, cctx); + status |= gr_gens(gens, cctx); + if (var < gens->length) + status |= gr_set(xpow, gr_vec_entry_srcptr(gens, var, cctx), cctx); + else + status |= GR_UNABLE; + gr_vec_clear(gens, cctx); + status |= gr_pow_ui(xpow, xpow, (len >= 1) ? (ulong) (len - 1) : 0, cctx); + + /* round trip euler_to_ddx(ddx_to_euler(P_d)) = x^(len-1) P_d */ + gr_ore_poly_fit_length(scaled_P, len, ctx_d); + for (i = 0; i < len; i++) + status |= gr_mul(GR_ENTRY(scaled_P->coeffs, i, sz), xpow, GR_ENTRY(P_d->coeffs, i, sz), cctx); + _gr_ore_poly_set_length(scaled_P, len, ctx_d); + _gr_ore_poly_normalise(scaled_P, ctx_d); + + if (status == GR_SUCCESS && gr_ore_poly_equal(P_d2, scaled_P, ctx_d) == T_FALSE) + { + flint_printf("FAIL: ddx_to_euler round trip\n"); + flint_abort(); + } + + for (j = 0; j < 4; j++) + { + status |= gr_randtest(f, state, cctx); + status |= gr_ore_poly_apply(g_d, P_d, f, ctx_d); + status |= gr_ore_poly_apply(g_e, P_e, f, ctx_e); + + /* P_e(f) = x^(len-1) P_d(f) */ + status |= gr_mul(scaled, xpow, g_d, cctx); + if (status == GR_SUCCESS && gr_equal(g_e, scaled, cctx) == T_FALSE) + { + flint_printf("FAIL: ddx_to_euler application identity\n"); + flint_abort(); + } + } + + /* these representative base rings never fail for the implemented var 0 */ + int expect_success = (var == 0 && cctx->which_ring == GR_CTX_GR_POLY + && (POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_FMPZ + || POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_CC_ACB + || POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_NMOD)); + if (expect_success && status != GR_SUCCESS) + { + flint_printf("FAIL: ddx_to_euler unexpected failure\n"); + flint_abort(); + } + + count_success += (status == GR_SUCCESS); + count_domain += ((status & GR_DOMAIN) != 0); + count_unable += ((status & GR_UNABLE) != 0); + + gr_heap_clear(xpow, cctx); + gr_heap_clear(f, cctx); + gr_heap_clear(g_d, cctx); + gr_heap_clear(g_e, cctx); + gr_heap_clear(scaled, cctx); + gr_ore_poly_clear(P_d, ctx_d); + gr_ore_poly_clear(P_e, ctx_e); + gr_ore_poly_clear(P_d2, ctx_d); + gr_ore_poly_clear(scaled_P, ctx_d); + gr_ore_poly_ctx_clear(ctx_d); + gr_ore_poly_ctx_clear(ctx_e); + gr_ctx_clear(cctx); + } + + TEST_GR_FUNCTION_END(state, count_success, count_domain, count_unable); +} diff --git a/src/gr_ore_poly/test/t-differential_shift.c b/src/gr_ore_poly/test/t-differential_shift.c new file mode 100644 index 0000000000..68ac9411c1 --- /dev/null +++ b/src/gr_ore_poly/test/t-differential_shift.c @@ -0,0 +1,242 @@ +/* + 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 . +*/ + +/* generated using Claude Opus 4.8 */ + +#include "test_helpers.h" +#include "gr.h" +#include "gr_poly.h" +#include "gr_ore_poly.h" + +static const ore_algebra_t ds_diff_algs[2] = { + ORE_ALGEBRA_DERIVATIVE, ORE_ALGEBRA_EULER_DERIVATIVE +}; + +static const ore_algebra_t ds_shift_algs[4] = { + ORE_ALGEBRA_FORWARD_SHIFT, ORE_ALGEBRA_BACKWARD_SHIFT, + ORE_ALGEBRA_FORWARD_DIFFERENCE, ORE_ALGEBRA_BACKWARD_DIFFERENCE +}; + +/* Returns 0 only if a == x^net * b is provably false (x = gen of cctx). */ +static int +check_scaled(const gr_ore_poly_t a, const gr_ore_poly_t b, slong net, + gr_ore_poly_ctx_t octx, gr_ctx_t cctx) +{ + int ok = 1, status = GR_SUCCESS; + slong sz = cctx->sizeof_elem, i; + const gr_ore_poly_struct * lo = (net >= 0) ? b : a; + const gr_ore_poly_struct * hi = (net >= 0) ? a : b; + gr_ptr xp; + gr_ore_poly_t scaled; + + GR_TMP_INIT(xp, cctx); + gr_ore_poly_init(scaled, octx); + + status |= gr_gen(xp, cctx); + status |= gr_pow_ui(xp, xp, (ulong) FLINT_ABS(net), cctx); + gr_ore_poly_fit_length(scaled, lo->length, octx); + for (i = 0; i < lo->length; i++) + status |= gr_mul(GR_ENTRY(scaled->coeffs, i, sz), xp, GR_ENTRY(lo->coeffs, i, sz), cctx); + _gr_ore_poly_set_length(scaled, lo->length, octx); + _gr_ore_poly_normalise(scaled, octx); + + if (status == GR_SUCCESS && gr_ore_poly_equal(scaled, hi, octx) == T_FALSE) + ok = 0; + + gr_ore_poly_clear(scaled, octx); + GR_TMP_CLEAR(xp, cctx); + return ok; +} + +/* Apply a forward-shift recurrence R = sum_k q_k(nu) S^k to the coefficients + (a_m = coeff_m(f)) of a polynomial: s_m = sum_k q_k(m) a_{m+k}. */ +static int +apply_forward_recurrence(gr_poly_t s, const gr_ore_poly_t R, const gr_poly_t f, gr_ctx_t cctx) +{ + gr_ctx_struct * sctx = POLYNOMIAL_ELEM_CTX(cctx); + slong bsz = cctx->sizeof_elem, ssz = sctx->sizeof_elem; + slong df = f->length, m, k; + int status = GR_SUCCESS; + gr_ptr mval, qval, term; + + GR_TMP_INIT3(mval, qval, term, sctx); + gr_poly_fit_length(s, FLINT_MAX(df, 1), sctx); + for (m = 0; m < df; m++) + { + gr_ptr sm = GR_ENTRY(s->coeffs, m, ssz); + status |= gr_zero(sm, sctx); + status |= gr_set_si(mval, m, sctx); + for (k = 0; k < R->length && m + k < df; k++) + { + const gr_poly_struct * qk = (const gr_poly_struct *) GR_ENTRY(R->coeffs, k, bsz); + status |= gr_poly_evaluate(qval, qk, mval, sctx); + status |= gr_mul(term, qval, GR_ENTRY(f->coeffs, m + k, ssz), sctx); + status |= gr_add(sm, sm, term, sctx); + } + } + _gr_poly_set_length(s, df, sctx); + _gr_poly_normalise(s, sctx); + GR_TMP_CLEAR3(mval, qval, term, sctx); + return status; +} + +TEST_GR_FUNCTION_START(gr_ore_poly_differential_to_shift, state, count_success, count_domain, count_unable) +{ + for (slong iter = 0; iter < 1000 * flint_test_multiplier(); iter++) + { + gr_ctx_t cctx; + gr_ore_poly_ctx_t ctx_d, ctx_s, ctx_fwd, ctx_eul, ctx_back; + ore_algebra_t diff_alg, shift_alg; + slong p1, p2, pf, pL, pM, pLM, ngens, var; + int status = GR_SUCCESS; + + /* occasionally use a multivariate base ring with a non-default + generator; the bridges are only implemented over a univariate + polynomial base and return GR_UNABLE otherwise */ + switch (n_randint(state, 8)) + { + case 0: + gr_ctx_init_random(cctx, state); + break; + case 1: + case 2: + gr_ctx_init_random_mpoly(cctx, state); + break; + default: + gr_ctx_init_random_poly(cctx, state); + break; + } + + ngens = 0; + status |= gr_ctx_ngens(&ngens, cctx); + var = (ngens >= 1) ? (slong) n_randint(state, ngens) : 0; + + diff_alg = ds_diff_algs[n_randint(state, 2)]; + shift_alg = ds_shift_algs[n_randint(state, 4)]; + + gr_ore_poly_ctx_init(ctx_d, cctx, var, diff_alg); + gr_ore_poly_ctx_init(ctx_s, cctx, var, shift_alg); + gr_ore_poly_ctx_init(ctx_fwd, cctx, var, ORE_ALGEBRA_FORWARD_SHIFT); + gr_ore_poly_ctx_init(ctx_eul, cctx, var, ORE_ALGEBRA_EULER_DERIVATIVE); + gr_ore_poly_ctx_init(ctx_back, cctx, var, ORE_ALGEBRA_BACKWARD_SHIFT); + + gr_ore_poly_t op, res_s, op2; + gr_ore_poly_init(op, ctx_d); + gr_ore_poly_init(res_s, ctx_s); + gr_ore_poly_init(op2, ctx_d); + + status |= gr_ore_poly_randtest(op, state, 1 + n_randint(state, 5), ctx_d); + + /* round trip recovers x^(p1-p2) * op (the corrections are left powers) */ + status |= gr_ore_poly_differential_to_shift(res_s, &p1, op, ctx_s, ctx_d); + status |= gr_ore_poly_shift_to_differential(op2, &p2, res_s, ctx_d, ctx_s); + if (status == GR_SUCCESS && !check_scaled(op2, op, p1 - p2, ctx_d, cctx)) + { + flint_printf("FAIL: differential<->shift round trip\n"); + flint_abort(); + } + + /* series action: a forward-shift recurrence R from op satisfies + (R a)_m = coeff_{m-pf}(op . f) for the coefficient sequence of f */ + if (cctx->which_ring == GR_CTX_GR_POLY) + { + gr_ctx_struct * sctx = POLYNOMIAL_ELEM_CTX(cctx); + gr_ore_poly_t Rf; + gr_ptr f = gr_heap_init(cctx); + gr_ptr g = gr_heap_init(cctx); + gr_poly_t s, eg; + + gr_ore_poly_init(Rf, ctx_fwd); + gr_poly_init(s, sctx); + gr_poly_init(eg, sctx); + + status |= gr_ore_poly_differential_to_shift(Rf, &pf, op, ctx_fwd, ctx_d); + status |= gr_randtest(f, state, cctx); + status |= gr_ore_poly_apply(g, op, f, ctx_d); + status |= apply_forward_recurrence(s, Rf, (gr_poly_struct *) f, cctx); + if (pf >= 0) + status |= gr_poly_shift_left(eg, (gr_poly_struct *) g, pf, sctx); + else + status |= gr_poly_shift_right(eg, (gr_poly_struct *) g, -pf, sctx); + if (status == GR_SUCCESS && gr_poly_equal(s, eg, sctx) == T_FALSE) + { + flint_printf("FAIL: differential->shift series action\n"); + flint_abort(); + } + + gr_ore_poly_clear(Rf, ctx_fwd); + gr_poly_clear(s, sctx); + gr_poly_clear(eg, sctx); + gr_heap_clear(f, cctx); + gr_heap_clear(g, cctx); + } + + /* the Euler -> backward shift bridge is an exact algebra homomorphism */ + gr_ore_poly_t L, M, LM, Lb, Mb, LMb, LbMb; + gr_ore_poly_init(L, ctx_eul); + gr_ore_poly_init(M, ctx_eul); + gr_ore_poly_init(LM, ctx_eul); + gr_ore_poly_init(Lb, ctx_back); + gr_ore_poly_init(Mb, ctx_back); + gr_ore_poly_init(LMb, ctx_back); + gr_ore_poly_init(LbMb, ctx_back); + + status |= gr_ore_poly_randtest(L, state, 1 + n_randint(state, 3), ctx_eul); + status |= gr_ore_poly_randtest(M, state, 1 + n_randint(state, 3), ctx_eul); + status |= gr_ore_poly_mul(LM, L, M, ctx_eul); + status |= gr_ore_poly_differential_to_shift(Lb, &pL, L, ctx_back, ctx_eul); + status |= gr_ore_poly_differential_to_shift(Mb, &pM, M, ctx_back, ctx_eul); + status |= gr_ore_poly_differential_to_shift(LMb, &pLM, LM, ctx_back, ctx_eul); + status |= gr_ore_poly_mul(LbMb, Lb, Mb, ctx_back); + if (status == GR_SUCCESS && (pL != 0 || pM != 0 || pLM != 0)) + { + flint_printf("FAIL: Euler -> backward shift not exact\n"); + flint_abort(); + } + if (status == GR_SUCCESS && gr_ore_poly_equal(LMb, LbMb, ctx_back) == T_FALSE) + { + flint_printf("FAIL: bridge not multiplicative\n"); + flint_abort(); + } + + /* these representative base rings never fail for the implemented var 0 */ + int expect_success = (var == 0 && cctx->which_ring == GR_CTX_GR_POLY + && (POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_FMPZ + || POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_CC_ACB + || POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_NMOD)); + if (expect_success && status != GR_SUCCESS) + { + flint_printf("FAIL: differential<->shift unexpected failure\n"); + flint_abort(); + } + + count_success += (status == GR_SUCCESS); + count_domain += ((status & GR_DOMAIN) != 0); + count_unable += ((status & GR_UNABLE) != 0); + + gr_ore_poly_clear(op, ctx_d); + gr_ore_poly_clear(res_s, ctx_s); + gr_ore_poly_clear(op2, ctx_d); + gr_ore_poly_clear(L, ctx_eul); + gr_ore_poly_clear(M, ctx_eul); + gr_ore_poly_clear(LM, ctx_eul); + gr_ore_poly_clear(Lb, ctx_back); + gr_ore_poly_clear(Mb, ctx_back); + gr_ore_poly_clear(LMb, ctx_back); + gr_ore_poly_clear(LbMb, ctx_back); + gr_ore_poly_ctx_clear(ctx_d); + gr_ore_poly_ctx_clear(ctx_s); + gr_ore_poly_ctx_clear(ctx_fwd); + gr_ore_poly_ctx_clear(ctx_eul); + gr_ore_poly_ctx_clear(ctx_back); + gr_ctx_clear(cctx); + } + + TEST_GR_FUNCTION_END(state, count_success, count_domain, count_unable); +} diff --git a/src/gr_ore_poly/test/t-euler_to_ddx.c b/src/gr_ore_poly/test/t-euler_to_ddx.c new file mode 100644 index 0000000000..c5693537fe --- /dev/null +++ b/src/gr_ore_poly/test/t-euler_to_ddx.c @@ -0,0 +1,142 @@ +/* + 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 . +*/ + +/* generated using Claude Opus 4.8 */ + +#include "test_helpers.h" +#include "gr.h" +#include "gr_vec.h" +#include "gr_ore_poly.h" + +TEST_GR_FUNCTION_START(gr_ore_poly_euler_to_ddx, state, count_success, count_domain, count_unable) +{ + for (slong iter = 0; iter < 1000 * flint_test_multiplier(); iter++) + { + gr_ctx_t cctx; + gr_ore_poly_ctx_t ctx_d, ctx_e; + gr_vec_t gens; + slong len, i, j, sz, ngens, var; + int status = GR_SUCCESS; + + /* occasionally use a multivariate base ring so that the conversions are + exercised with a generator other than the default one */ + switch (n_randint(state, 8)) + { + case 0: + gr_ctx_init_random(cctx, state); + break; + case 1: + case 2: + gr_ctx_init_random_mpoly(cctx, state); + break; + default: + gr_ctx_init_random_poly(cctx, state); + break; + } + + sz = cctx->sizeof_elem; + + ngens = 0; + status |= gr_ctx_ngens(&ngens, cctx); + var = (ngens >= 1) ? (slong) n_randint(state, ngens) : 0; + + gr_ore_poly_ctx_init(ctx_d, cctx, var, ORE_ALGEBRA_DERIVATIVE); + gr_ore_poly_ctx_init(ctx_e, cctx, var, ORE_ALGEBRA_EULER_DERIVATIVE); + + gr_ore_poly_t P_e, P_d, P_e2, scaled_P; + gr_ore_poly_init(P_e, ctx_e); + gr_ore_poly_init(P_d, ctx_d); + gr_ore_poly_init(P_e2, ctx_e); + gr_ore_poly_init(scaled_P, ctx_e); + + status |= gr_ore_poly_randtest(P_e, state, 1 + n_randint(state, 5), ctx_e); + len = P_e->length; + + gr_ore_poly_fit_length(P_d, len, ctx_d); + status |= _gr_ore_poly_euler_to_ddx(P_d->coeffs, P_e->coeffs, len, var, cctx); + _gr_ore_poly_set_length(P_d, len, ctx_d); + _gr_ore_poly_normalise(P_d, ctx_d); + + gr_ore_poly_fit_length(P_e2, len, ctx_e); + status |= _gr_ore_poly_ddx_to_euler(P_e2->coeffs, P_d->coeffs, len, var, cctx); + _gr_ore_poly_set_length(P_e2, len, ctx_e); + _gr_ore_poly_normalise(P_e2, ctx_e); + + gr_ptr xpow = gr_heap_init(cctx); + gr_ptr f = gr_heap_init(cctx); + gr_ptr g_e = gr_heap_init(cctx); + gr_ptr g_d = gr_heap_init(cctx); + + /* x is the generator of index var */ + gr_vec_init(gens, 0, cctx); + status |= gr_gens(gens, cctx); + if (var < gens->length) + status |= gr_set(xpow, gr_vec_entry_srcptr(gens, var, cctx), cctx); + else + status |= GR_UNABLE; + gr_vec_clear(gens, cctx); + status |= gr_pow_ui(xpow, xpow, (len >= 1) ? (ulong) (len - 1) : 0, cctx); + + /* round trip ddx_to_euler(euler_to_ddx(P_e)) = x^(len-1) P_e */ + gr_ore_poly_fit_length(scaled_P, len, ctx_e); + for (i = 0; i < len; i++) + status |= gr_mul(GR_ENTRY(scaled_P->coeffs, i, sz), xpow, GR_ENTRY(P_e->coeffs, i, sz), cctx); + _gr_ore_poly_set_length(scaled_P, len, ctx_e); + _gr_ore_poly_normalise(scaled_P, ctx_e); + + if (status == GR_SUCCESS && gr_ore_poly_equal(P_e2, scaled_P, ctx_e) == T_FALSE) + { + flint_printf("FAIL: euler_to_ddx round trip\n"); + flint_abort(); + } + + for (j = 0; j < 4; j++) + { + status |= gr_randtest(f, state, cctx); + status |= gr_ore_poly_apply(g_e, P_e, f, ctx_e); + status |= gr_ore_poly_apply(g_d, P_d, f, ctx_d); + + /* euler_to_ddx leaves the operator unchanged: P_d(f) = P_e(f) */ + if (status == GR_SUCCESS && gr_equal(g_d, g_e, cctx) == T_FALSE) + { + flint_printf("FAIL: euler_to_ddx application identity\n"); + flint_abort(); + } + } + + /* these representative base rings never fail for the implemented var 0 */ + int expect_success = (var == 0 && cctx->which_ring == GR_CTX_GR_POLY + && (POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_FMPZ + || POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_CC_ACB + || POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_NMOD)); + if (expect_success && status != GR_SUCCESS) + { + flint_printf("FAIL: euler_to_ddx unexpected failure\n"); + flint_abort(); + } + + count_success += (status == GR_SUCCESS); + count_domain += ((status & GR_DOMAIN) != 0); + count_unable += ((status & GR_UNABLE) != 0); + + gr_heap_clear(xpow, cctx); + gr_heap_clear(f, cctx); + gr_heap_clear(g_e, cctx); + gr_heap_clear(g_d, cctx); + gr_ore_poly_clear(P_e, ctx_e); + gr_ore_poly_clear(P_d, ctx_d); + gr_ore_poly_clear(P_e2, ctx_e); + gr_ore_poly_clear(scaled_P, ctx_e); + gr_ore_poly_ctx_clear(ctx_d); + gr_ore_poly_ctx_clear(ctx_e); + gr_ctx_clear(cctx); + } + + TEST_GR_FUNCTION_END(state, count_success, count_domain, count_unable); +} diff --git a/src/gr_ore_poly/test/t-shift_convert.c b/src/gr_ore_poly/test/t-shift_convert.c new file mode 100644 index 0000000000..b30b7979dc --- /dev/null +++ b/src/gr_ore_poly/test/t-shift_convert.c @@ -0,0 +1,154 @@ +/* + 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 . +*/ + +/* generated using Claude Opus 4.8 */ + +#include "test_helpers.h" +#include "gr.h" +#include "gr_poly.h" +#include "gr_ore_poly.h" + +static const ore_algebra_t shift_algs[4] = { + ORE_ALGEBRA_FORWARD_SHIFT, ORE_ALGEBRA_BACKWARD_SHIFT, + ORE_ALGEBRA_FORWARD_DIFFERENCE, ORE_ALGEBRA_BACKWARD_DIFFERENCE +}; + +TEST_GR_FUNCTION_START(gr_ore_poly_shift_convert, state, count_success, count_domain, count_unable) +{ + for (slong iter = 0; iter < 1000 * flint_test_multiplier(); iter++) + { + gr_ctx_t cctx; + gr_ore_poly_ctx_t ctx_src, ctx_dst; + ore_algebra_t src_alg, dst_alg; + slong len, j, epow, epow2, ngens, var; + int status = GR_SUCCESS; + + /* occasionally use a multivariate base ring so that the conversion is + exercised with a generator other than the default one (which is not + yet implemented and returns GR_UNABLE) */ + switch (n_randint(state, 8)) + { + case 0: + gr_ctx_init_random(cctx, state); + break; + case 1: + case 2: + gr_ctx_init_random_mpoly(cctx, state); + break; + default: + gr_ctx_init_random_poly(cctx, state); + break; + } + + ngens = 0; + status |= gr_ctx_ngens(&ngens, cctx); + var = (ngens >= 1) ? (slong) n_randint(state, ngens) : 0; + + src_alg = shift_algs[n_randint(state, 4)]; + dst_alg = shift_algs[n_randint(state, 4)]; + + gr_ore_poly_ctx_init(ctx_src, cctx, var, src_alg); + gr_ore_poly_ctx_init(ctx_dst, cctx, var, dst_alg); + + gr_ore_poly_t P, Pd, P2; + gr_ore_poly_init(P, ctx_src); + gr_ore_poly_init(Pd, ctx_dst); + gr_ore_poly_init(P2, ctx_src); + + status |= gr_ore_poly_randtest(P, state, 1 + n_randint(state, 5), ctx_src); + len = P->length; + + gr_ore_poly_fit_length(Pd, len, ctx_dst); + status |= _gr_ore_poly_shift_convert(Pd->coeffs, &epow, P->coeffs, len, src_alg, dst_alg, var, cctx); + _gr_ore_poly_set_length(Pd, len, ctx_dst); + _gr_ore_poly_normalise(Pd, ctx_dst); + + gr_ore_poly_fit_length(P2, len, ctx_src); + status |= _gr_ore_poly_shift_convert(P2->coeffs, &epow2, Pd->coeffs, len, dst_alg, src_alg, var, cctx); + _gr_ore_poly_set_length(P2, len, ctx_src); + _gr_ore_poly_normalise(P2, ctx_src); + + /* these representative base rings never fail for the implemented var 0 */ + int expect_success = (var == 0 && cctx->which_ring == GR_CTX_GR_POLY + && (POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_FMPZ + || POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_CC_ACB + || POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_NMOD)); + if (expect_success && status != GR_SUCCESS) + { + flint_printf("FAIL: shift_convert unexpected failure\n"); + flint_abort(); + } + + /* the round trip recovers the operator: the powers of S cancel */ + if (status == GR_SUCCESS && epow + epow2 != 0) + { + flint_printf("FAIL: shift_convert round-trip power\n"); + flint_abort(); + } + if (status == GR_SUCCESS && gr_ore_poly_equal(P2, P, ctx_src) == T_FALSE) + { + flint_printf("FAIL: shift_convert round trip\n"); + flint_abort(); + } + + /* application identity dst(f) = S^{-epow}(src(f)), over polynomial base + rings (S^epow * dst = src as operators, so the correction is S^{-epow} + applied on the left of the value) */ + if (cctx->which_ring == GR_CTX_GR_POLY) + { + gr_ctx_struct * sctx = POLYNOMIAL_ELEM_CTX(cctx); + gr_ptr f = gr_heap_init(cctx); + gr_ptr shifted = gr_heap_init(cctx); + gr_ptr g1 = gr_heap_init(cctx); + gr_ptr g2 = gr_heap_init(cctx); + gr_ptr c = gr_heap_init(sctx); + + status |= gr_set_si(c, -epow, sctx); + + for (j = 0; j < 4; j++) + { + status |= gr_randtest(f, state, cctx); + status |= gr_ore_poly_apply(g1, Pd, f, ctx_dst); + status |= gr_ore_poly_apply(g2, P, f, ctx_src); + status |= gr_poly_taylor_shift((gr_poly_struct *) shifted, (gr_poly_struct *) g2, c, sctx); + if (status == GR_SUCCESS && gr_equal(g1, shifted, cctx) == T_FALSE) + { + flint_printf("FAIL: shift_convert application identity\n"); + flint_abort(); + } + } + + gr_heap_clear(f, cctx); + gr_heap_clear(shifted, cctx); + gr_heap_clear(g1, cctx); + gr_heap_clear(g2, cctx); + gr_heap_clear(c, sctx); + } + + /* a non-shift algebra type is rejected cleanly */ + if (_gr_ore_poly_shift_convert(P2->coeffs, &epow2, P->coeffs, len, ORE_ALGEBRA_DERIVATIVE, dst_alg, var, cctx) != GR_DOMAIN) + { + flint_printf("FAIL: shift_convert expected GR_DOMAIN\n"); + flint_abort(); + } + + count_success += (status == GR_SUCCESS); + count_domain += ((status & GR_DOMAIN) != 0); + count_unable += ((status & GR_UNABLE) != 0); + + gr_ore_poly_clear(P, ctx_src); + gr_ore_poly_clear(Pd, ctx_dst); + gr_ore_poly_clear(P2, ctx_src); + gr_ore_poly_ctx_clear(ctx_src); + gr_ore_poly_ctx_clear(ctx_dst); + gr_ctx_clear(cctx); + } + + TEST_GR_FUNCTION_END(state, count_success, count_domain, count_unable); +} From 070023b5b65be024a7f9e0acc903a6e81ec79f9e Mon Sep 17 00:00:00 2001 From: Marc Mezzarobba Date: Wed, 1 Jul 2026 13:14:26 +0200 Subject: [PATCH 3/7] gr_ore_poly: revise AI-generated code --- doc/source/gr_ore_poly.rst | 201 ++++++++-------- src/gr_ore_poly.h | 38 +-- src/gr_ore_poly/apply.c | 10 +- src/gr_ore_poly/backshift_to_euler_univar.c | 31 ++- src/gr_ore_poly/convert.c | 68 +++--- src/gr_ore_poly/ddx_to_euler.c | 29 ++- src/gr_ore_poly/differential_to_shift.c | 59 ++--- src/gr_ore_poly/euler_to_backshift_univar.c | 17 +- src/gr_ore_poly/euler_to_ddx.c | 26 +-- src/gr_ore_poly/shift_convert.c | 126 +++++----- src/gr_ore_poly/shift_to_differential.c | 75 +++--- src/gr_ore_poly/test/main.c | 2 +- src/gr_ore_poly/test/t-apply.c | 217 +++++++++--------- src/gr_ore_poly/test/t-ddx_to_euler.c | 67 +++--- ...tial_shift.c => t-differential_to_shift.c} | 124 ++++------ src/gr_ore_poly/test/t-euler_to_ddx.c | 76 +++--- src/gr_ore_poly/test/t-shift_convert.c | 67 +++--- 17 files changed, 555 insertions(+), 678 deletions(-) rename src/gr_ore_poly/test/{t-differential_shift.c => t-differential_to_shift.c} (61%) diff --git a/doc/source/gr_ore_poly.rst b/doc/source/gr_ore_poly.rst index a15a0a59ba..a422449f62 100644 --- a/doc/source/gr_ore_poly.rst +++ b/doc/source/gr_ore_poly.rst @@ -3,10 +3,6 @@ **gr_ore_poly.h** -- dense univariate Ore polynomials over generic rings =============================================================================== -.. note:: - - This module is under construction. - A :type:`gr_ore_poly_t` represents a univariate Ore polynomial `L \in R[D]` implemented as a dense array of coefficients in a generic ring *R*. The choice of Ore algebra structure (e.g. with `D` being the standard @@ -27,15 +23,15 @@ Ore algebra types Represents one of the following supported Ore algebra types: - .. macro:: ORE_ALGEBRA_CUSTOM + .. enumerator:: ORE_ALGEBRA_CUSTOM Custom Ore polynomials. - .. macro:: ORE_ALGEBRA_COMMUTATIVE + .. enumerator:: ORE_ALGEBRA_COMMUTATIVE Standard polynomials. - .. macro:: ORE_ALGEBRA_DERIVATIVE + .. enumerator:: ORE_ALGEBRA_DERIVATIVE Linear differential operators in the standard derivative. @@ -43,7 +39,7 @@ Ore algebra types `\delta` is the derivative `\frac{d}{dx}` with respect to a generator `x` of the base ring. - .. macro:: ORE_ALGEBRA_EULER_DERIVATIVE + .. enumerator:: ORE_ALGEBRA_EULER_DERIVATIVE Linear differential operators in the Euler derivative. @@ -51,7 +47,7 @@ Ore algebra types `\delta` is the Euler derivative `x\cdot\frac{d}{dx}` with respect to a generator `x` of the base ring. - .. macro:: ORE_ALGEBRA_FORWARD_SHIFT + .. enumerator:: ORE_ALGEBRA_FORWARD_SHIFT Linear difference operators in the standard forward shift. @@ -59,7 +55,7 @@ Ore algebra types to a generator `x` of the base ring, and the `\sigma`-derivation `\delta` is the zero map. - .. macro:: ORE_ALGEBRA_FORWARD_DIFFERENCE + .. enumerator:: ORE_ALGEBRA_FORWARD_DIFFERENCE Linear difference operator in the forward finite difference operator. @@ -67,23 +63,23 @@ Ore algebra types to a generator `x` of the base ring, and the `\sigma`-derivation `\delta` maps `x \mapsto 1`. - .. macro:: ORE_ALGEBRA_BACKWARD_SHIFT + .. enumerator:: ORE_ALGEBRA_BACKWARD_SHIFT Linear difference operators in the standard backward shift. - .. macro:: ORE_ALGEBRA_BACKWARD_DIFFERENCE + .. enumerator:: ORE_ALGEBRA_BACKWARD_DIFFERENCE Linear difference operator in the backward finite difference operator. - .. macro:: ORE_ALGEBRA_Q_SHIFT + .. enumerator:: ORE_ALGEBRA_Q_SHIFT Linear q-difference operators. - .. macro:: ORE_ALGEBRA_MAHLER + .. enumerator:: ORE_ALGEBRA_MAHLER Linear Mahler operators. - .. macro:: ORE_ALGEBRA_FROBENIUS + .. enumerator:: ORE_ALGEBRA_FROBENIUS Ore polynomials over a field twisted by the Frobenius endomorphism. @@ -266,109 +262,112 @@ Action A pointer to a function with the same specification as :func:`gr_ore_poly_sigma_delta`. +.. function:: int gr_ore_poly_apply(gr_ptr res, const gr_ore_poly_t P, gr_srcptr f, gr_ore_poly_ctx_t ctx) + + Sets *res* to the result of applying *P* to the base ring element *f* under + the standard interpretation of *P* as an operator acting on the base ring + (derivative operators differentiate, shift operators shift, etc.). + .. function:: int gr_ore_poly_apply_custom(gr_ptr res, const gr_ore_poly_t P, gr_srcptr f, gr_srcptr d1, gr_ore_poly_ctx_t ctx) Sets *res* to the result of applying *P* to the base ring element *f*, where the generator `D` acts by `g \mapsto \sigma(g) \cdot d1 + \delta(g)` for the - given value ``d1`` of `D(1)`. Any ``d1`` defines a valid action. Aliasing of - *res* with *f* or *d1* is allowed. - -.. function:: int gr_ore_poly_apply(gr_ptr res, const gr_ore_poly_t P, gr_srcptr f, gr_ore_poly_ctx_t ctx) - - As :func:`gr_ore_poly_apply_custom`, but with the standard `D(1)` for the - algebra type (so derivative operators differentiate, shift operators shift, - etc.). Returns ``GR_UNABLE`` when no standard action exists - (`ORE_ALGEBRA_COMMUTATIVE`, `ORE_ALGEBRA_CUSTOM`) or it cannot be computed. + given value *d1* of `D(1)`. Any *d1* defines a valid action. Conversions ------------------------------------------------------------------------------- -These functions rewrite an operator from one Ore algebra into another. The -underscore versions act on raw coefficient arrays over the *coefficient* ring -*ctx* (the base ring of the Ore ring, not the Ore ring itself); *x* denotes its -generator of index *var*, that is, entry *var* of :func:`gr_gens`. Unless stated -otherwise, *res* has the same length *len* as *op* and must not alias it. +The following functions convert between expressions of a linear differential +or difference operator in different bases, represented as Ore polynomials in +different Ore polynomial rings over the same base ring. -A conversion is an exact rewriting up to a unit factor. Since *x* and the -forward shift `S` need not act invertibly in the chosen representation, the -result may be the input premultiplied **on the left** by a power of *x* or `S`, -returned in an output parameter; left multiplication does not change the -solution space of the operator. Conversions that substitute into or -Taylor-shift the variable require a :macro:`GR_CTX_GR_POLY` base ring and return -``GR_UNABLE`` otherwise. +.. function:: int _gr_ore_poly_euler_to_ddx(gr_ptr res, gr_srcptr op, slong len, slong var, gr_ctx_t ctx) -The first pair rewrites between the derivative `D = d/dx` and the Euler -derivative `\theta = x D`. + Rewrites an Ore polynomial *op* of type :enumerator:`ORE_ALGEBRA_EULER_DERIVATIVE` + as an Ore polynomial of type :enumerator:`ORE_ALGEBRA_DERIVATIVE`. + The context *ctx* is the common base ring and *var* is the index of the + generator of *ctx* on which the derivations act. + The output vector *res* has the same length *len* as *op* and must not + alias it. .. function:: int _gr_ore_poly_ddx_to_euler(gr_ptr res, gr_srcptr op, slong len, slong var, gr_ctx_t ctx) - int _gr_ore_poly_euler_to_ddx(gr_ptr res, gr_srcptr op, slong len, slong var, gr_ctx_t ctx) - - Rewrite *op* between `D` and `\theta`. As `x` may not be invertible, - :func:`_gr_ore_poly_ddx_to_euler` returns `x^{len-1} \cdot \text{op}`; the - reverse rewriting needs no factor of `x`. - -The next pair converts among the four shift/difference algebras: the forward and -backward shifts `S`, `S^{-1}` and the forward and backward differences `S-1`, -`1-S^{-1}`. - -.. function:: int _gr_ore_poly_shift_convert(gr_ptr res, slong * epow, gr_srcptr op, slong len, ore_algebra_t src_alg, ore_algebra_t dst_alg, slong var, gr_ctx_t ctx) - - Convert *op* from *src_alg* to *dst_alg*. Each generator is a degree-one - Laurent polynomial in `S`, so on output ``*epow`` `= p` and *res* satisfy - `\text{res}_{dst} = S^{p} \cdot \text{op}_{src}`. A nonzero `p` (hence a - polynomial base ring) is needed only when crossing between the forward side - `S`, `S-1` and the backward side `S^{-1}`, `1-S^{-1}`, where `p = \mp(len-1)`. - Converting between the two difference types uses a single binomial transform. - The variable is the base ring generator of index *var*; only ``var == 0`` - over a univariate polynomial ring is implemented, and other indices return - ``GR_UNABLE``. Returns ``GR_DOMAIN`` if either algebra is not a - shift/difference type. - -The remaining functions bridge differential and shift operators through the -generating-series isomorphism `\theta \mapsto n`, `x \mapsto S^{-1}`: for -`f = \sum_n a_n x^n`, the Euler derivative `\theta = x d/dx` acts on `(a_n)` as -multiplication by `n` and `x` as the backward shift `S^{-1}`. All require a -`GR_CTX_GR_POLY` base ring and exchange the operator order with the coefficient -degree. + + Rewrites an Ore polynomial *op* of type :enumerator:`ORE_ALGEBRA_DERIVATIVE` + as an Ore polynomial *res* of type + :enumerator:`ORE_ALGEBRA_EULER_DERIVATIVE` such that + `\mathit{res} = x^{len-1} \cdot \mathit{op}`, + where `x` is the generator of index *var* of the base ring *ctx*. + The output vector *res* has the same length *len* as *op* and must not + alias it. + +.. function:: int _gr_ore_poly_shift_convert(gr_ptr res, slong * p, gr_srcptr op, slong len, ore_algebra_t src_alg, ore_algebra_t dst_alg, slong var, gr_ctx_t ctx) + + Rewrites an operator *op* from *src_alg* to *dst_alg* where *src_alg* and + *dst_alg* are among the builtin shift and difference algebras (corresponding + to operators written in terms of the forward and backward shifts `S`, + `S^{-1}` and the forward and backward differences `S-1`, `1-S^{-1}`). + The context *ctx* is the common base ring and *var* is the index of the + generator of *ctx* on which `S` acts. Conversions that crosses between the + forward side `S`, `S-1` and the backward side `S^{-1}`, `1-S^{-1}` + additionally require a generic univariate polynomial base ring and otherwise + return ``GR_UNABLE``. + The result satisfies + `S^{\textit{p}} \cdot \textit{res} = \textit{op}`. + The output vector *res* has the same length *len* as *op* and must not + alias it. + This function returns an error status when the source or destination algebra + is not of the required type. + +.. function:: int _gr_ore_poly_shift_convert_difference(gr_ptr res, slong * p, gr_srcptr op, slong len, int to_backward, slong var, gr_ctx_t ctx) + + Specialized version of :func:`_gr_ore_poly_shift_convert` for converting + between :enumerator:`ORE_ALGEBRA_FORWARD_DIFFERENCE` and + :enumerator:`ORE_ALGEBRA_BACKWARD_DIFFERENCE` or back. The *to_backward* + flag indicates the direction of the conversion. + +.. function:: int gr_ore_poly_convert(gr_ore_poly_t res, slong * p, const gr_ore_poly_t op, gr_ore_poly_ctx_t res_ctx, gr_ore_poly_ctx_t op_ctx) + + Convert *op* from *op_ctx* to *res_ctx*. + The meaning of the output parameter *p* is algebra-dependent. + For a conversion within the differential family one has + `x^{p} \cdot \textit{res} = \textit{op}` + where `x` is the generator of the base ring specified in the source context. + For a conversion within the shift/difference family, one has + `S^{p} \cdot \textit{res} = \textit{op}` (a power of the forward shift `S`). + No attempt is currently made to minimize *p* or its absolute value. + +For `f = \sum_n a_n x^n`, the Euler derivative `\theta = x d/dx` acts on +`(a_n)` as multiplication by `n` and `x` acts as the backward shift `S^{-1}`. +The following functions convert between differential and difference operators in +a way compatible with this action, mapping `x d/dx \mapsto n`, +`x \mapsto S^{-1}` and inversely. .. function:: int _gr_ore_poly_euler_to_backshift_univar(gr_ptr res, slong reslen, gr_srcptr op, slong len, gr_ctx_t ctx) int _gr_ore_poly_backshift_to_euler_univar(gr_ptr res, slong reslen, gr_srcptr op, slong len, gr_ctx_t ctx) The two inverse rewritings of the isomorphism between the Euler operator and - the backward shift `S^{-1}`. The caller allocates *res* to *reslen*, one more - than the largest coefficient degree of *op*. - -.. function:: int gr_ore_poly_differential_to_shift(gr_ore_poly_t res, slong * power, const gr_ore_poly_t op, gr_ore_poly_ctx_t res_ctx, gr_ore_poly_ctx_t op_ctx) - int gr_ore_poly_shift_to_differential(gr_ore_poly_t res, slong * power, const gr_ore_poly_t op, gr_ore_poly_ctx_t res_ctx, gr_ore_poly_ctx_t op_ctx) - - Convert a whole operator between a differential algebra - (`ORE_ALGEBRA_DERIVATIVE` or `ORE_ALGEBRA_EULER_DERIVATIVE`) and a - shift/difference algebra, routing through the Euler and backward-shift pivots - above. On output ``*power`` `= p` is the left power absorbed (`S^p` for - :func:`gr_ore_poly_differential_to_shift`, `x^p` for - :func:`gr_ore_poly_shift_to_differential`), so `\text{op}(f) = 0` and - `\text{res} = 0` have the same solutions up to the index shift by `p`. Both - contexts must share a `GR_CTX_GR_POLY` base ring (else ``GR_UNABLE``); returns - ``GR_DOMAIN`` if a context algebra is outside its family. - -.. function:: int gr_ore_poly_convert(gr_ore_poly_t res, slong * power, const gr_ore_poly_t op, gr_ore_poly_ctx_t res_ctx, gr_ore_poly_ctx_t op_ctx) - - Convert *op* from the algebra of *op_ctx* to that of *res_ctx*, provided both - belong to the same family: the differential family - (`ORE_ALGEBRA_DERIVATIVE`, `ORE_ALGEBRA_EULER_DERIVATIVE`) or the - shift/difference family (the four shift and difference types). This - dispatches to the appropriate conversion above and handles memory - management. On output ``*power`` `= p` is the left power absorbed, **but its - meaning depends on the pair of algebras**: for a conversion within the - differential family `\text{res} = x^{p} \cdot \text{op}` (a power of the base - generator `x`), whereas for one within the shift/difference family - `\text{res} = S^{p} \cdot \text{op}` (a power of the forward shift `S`). A - conversion to the same algebra is the identity, with `p = 0`. Returns - ``GR_UNABLE`` for an unsupported pair: one that crosses the - differential/shift boundary (use :func:`gr_ore_poly_differential_to_shift` or - :func:`gr_ore_poly_shift_to_differential` for those) or that involves an - algebra outside the two families, such as `ORE_ALGEBRA_Q_SHIFT`, - `ORE_ALGEBRA_MAHLER`, `ORE_ALGEBRA_FROBENIUS` or `ORE_ALGEBRA_COMMUTATIVE`. + the backward shift `S^{-1}`. + The common base ring *ctx* must be a univariate polynomial ring. + The caller allocates *res* to *reslen*, one more than the largest + coefficient degree of *op*. + +.. function:: int gr_ore_poly_differential_to_shift(gr_ore_poly_t res, slong * p, const gr_ore_poly_t op, gr_ore_poly_ctx_t res_ctx, gr_ore_poly_ctx_t op_ctx) + + Given a differential operator *op* represented as an element of *op_ctx*, + computes a shift/difference operator *res* in *res_ctx* and an integer *p* + such that the above correspondence maps *op* to `S^p \cdot \textit{res}`. + The generators of the base rings specified in the source and + destination contexts play the role of `x` and `n`. + No attempt is currently made to minimize *p* or its absolute value. + +.. function:: int gr_ore_poly_shift_to_differential(gr_ore_poly_t res, slong * p, const gr_ore_poly_t op, gr_ore_poly_ctx_t res_ctx, gr_ore_poly_ctx_t op_ctx) + + Given a shift/difference operator *op* represented as an element of *op_ctx*, + computes a differential operator *res* in *res_ctx* and an integer *p* + such that the above correspondence maps *op* to `x^p \cdot \textit{res}`. + The generators of the base rings specified in the source and + destination contexts play the role of `n` and `x`. + No attempt is currently made to minimize *p* or its absolute value. Arithmetic ------------------------------------------------------------------------------- diff --git a/src/gr_ore_poly.h b/src/gr_ore_poly.h index efeecc84ab..3266fca112 100644 --- a/src/gr_ore_poly.h +++ b/src/gr_ore_poly.h @@ -250,50 +250,14 @@ extern const gr_ore_poly_sigma_delta_t _gr_ore_poly_default_sigma_delta[]; WARN_UNUSED_RESULT int gr_ore_poly_apply_custom(gr_ptr res, const gr_ore_poly_t P, gr_srcptr f, gr_srcptr d1, gr_ore_poly_ctx_t ctx); WARN_UNUSED_RESULT int gr_ore_poly_apply(gr_ptr res, const gr_ore_poly_t P, gr_srcptr f, gr_ore_poly_ctx_t ctx); -/* Conversions between operator types. - - These act on raw coefficient arrays and take the coefficient ring context - (not the Ore polynomial ring context). The variable x is the generator of the - coefficient ring of index var, i.e. entry var of gr_gens(ctx). Since x need - not be invertible, ddx_to_euler multiplies the operator by x^(len-1); - euler_to_ddx leaves it unchanged. The result res has the same length len as - op. */ WARN_UNUSED_RESULT int _gr_ore_poly_ddx_to_euler(gr_ptr res, gr_srcptr op, slong len, slong var, gr_ctx_t ctx); WARN_UNUSED_RESULT int _gr_ore_poly_euler_to_ddx(gr_ptr res, gr_srcptr op, slong len, slong var, gr_ctx_t ctx); - -/* Convert a length-len operator between the shift/difference algebra types - src_alg and dst_alg. Each generator is a degree-one Laurent polynomial in the - forward shift S, so on output *p and res satisfy - S^p * res_in_dst = op_in_src (a power of S appears only when crossing between - the forward side S, S-1 and the backward side S^-1, 1-S^-1; it then needs a - GR_CTX_GR_POLY base ring). Returns GR_DOMAIN if either algebra is not a - shift/difference type. The variable is the base ring generator of index var; - only var == 0 over a univariate polynomial ring is implemented (else - GR_UNABLE). */ WARN_UNUSED_RESULT int _gr_ore_poly_shift_convert(gr_ptr res, slong * p, gr_srcptr op, slong len, ore_algebra_t src_alg, ore_algebra_t dst_alg, slong var, gr_ctx_t ctx); - -/* Exact bridge of the isomorphism theta <-> n, x <-> backward shift S^{-1}, on - raw coefficient vectors over a GR_CTX_GR_POLY base ring (else GR_UNABLE). The - caller allocates res to reslen = 1 + the maximum coefficient degree of op. */ +WARN_UNUSED_RESULT int _gr_ore_poly_shift_convert_difference(gr_ptr res, slong * p, gr_srcptr op, slong len, int to_backward, slong var, gr_ctx_t ctx); WARN_UNUSED_RESULT int _gr_ore_poly_euler_to_backshift_univar(gr_ptr res, slong reslen, gr_srcptr op, slong len, gr_ctx_t ctx); WARN_UNUSED_RESULT int _gr_ore_poly_backshift_to_euler_univar(gr_ptr res, slong reslen, gr_srcptr op, slong len, gr_ctx_t ctx); - -/* Convert between differential operators (op_ctx/res_ctx in {DERIVATIVE, - EULER_DERIVATIVE}) and shift/difference operators (in the four shift types), - via the generating-series isomorphism theta <-> n, x <-> backward shift. - Both contexts must share a GR_CTX_GR_POLY base ring. On output *p is such that - C^p * res corresponds to op under the isomorphism, where C is the companion - generator on the res side (S for *_to_shift, x for *_to_differential). - Returns GR_DOMAIN if a context algebra is outside its family. */ WARN_UNUSED_RESULT int gr_ore_poly_differential_to_shift(gr_ore_poly_t res, slong * p, const gr_ore_poly_t op, gr_ore_poly_ctx_t res_ctx, gr_ore_poly_ctx_t op_ctx); WARN_UNUSED_RESULT int gr_ore_poly_shift_to_differential(gr_ore_poly_t res, slong * p, const gr_ore_poly_t op, gr_ore_poly_ctx_t res_ctx, gr_ore_poly_ctx_t op_ctx); - -/* Convert an operator between two algebras of the same family (both differential - or both shift/difference), dispatching to the conversions above. On output *p - satisfies C^p * res = op, where the companion generator C depends on the pair - (x for differential conversions, the forward shift S for shift ones). - Returns GR_UNABLE for an unsupported pair (crossing the differential/shift - boundary, or involving an algebra outside those two families). */ WARN_UNUSED_RESULT int gr_ore_poly_convert(gr_ore_poly_t res, slong * p, const gr_ore_poly_t op, gr_ore_poly_ctx_t res_ctx, gr_ore_poly_ctx_t op_ctx); /* Arithmetic */ diff --git a/src/gr_ore_poly/apply.c b/src/gr_ore_poly/apply.c index 1a98d1ac9e..d377b492d9 100644 --- a/src/gr_ore_poly/apply.c +++ b/src/gr_ore_poly/apply.c @@ -32,7 +32,7 @@ gr_ore_poly_apply_custom(gr_ptr res, const gr_ore_poly_t P, gr_srcptr f, gr_srcp GR_TMP_INIT5(g, acc, sig, del, term, base); status |= gr_set(g, f, base); - status |= gr_zero(acc, base); + status |= gr_zero(acc, base); /* todo: add underscore version? */ for (slong i = 0; i < len; i++) { @@ -44,10 +44,7 @@ gr_ore_poly_apply_custom(gr_ptr res, const gr_ore_poly_t P, gr_srcptr f, gr_srcp if (i + 1 < len) { if (d1_is_zero == T_TRUE) - { - status |= gr_ore_poly_delta(del, g, ctx); - status |= gr_set(g, del, base); - } + status |= gr_ore_poly_delta(g, g, ctx); else { status |= gr_ore_poly_sigma_delta(sig, del, g, ctx); @@ -74,7 +71,6 @@ gr_ore_poly_apply(gr_ptr res, const gr_ore_poly_t P, gr_srcptr f, gr_ore_poly_ct GR_TMP_INIT(d1, base); - /* Standard value of D(1) for the algebra type. */ switch (GR_ORE_POLY_CTX(ctx)->which_algebra) { case ORE_ALGEBRA_DERIVATIVE: @@ -90,7 +86,7 @@ gr_ore_poly_apply(gr_ptr res, const gr_ore_poly_t P, gr_srcptr f, gr_ore_poly_ct case ORE_ALGEBRA_FROBENIUS: status |= gr_one(d1, base); break; - default: /* ORE_ALGEBRA_COMMUTATIVE, ORE_ALGEBRA_CUSTOM: no standard action */ + default: status = GR_UNABLE; } diff --git a/src/gr_ore_poly/backshift_to_euler_univar.c b/src/gr_ore_poly/backshift_to_euler_univar.c index 8f6ff89cc9..18105d0c7b 100644 --- a/src/gr_ore_poly/backshift_to_euler_univar.c +++ b/src/gr_ore_poly/backshift_to_euler_univar.c @@ -13,11 +13,6 @@ #include "gr_poly.h" #include "gr_ore_poly.h" -/* Bridge from the backward shift operator (K[n]) to the Euler operator - (K[x]), the inverse of _gr_ore_poly_euler_to_backshift: Taylor-shift - column k by +k, then transpose the (S^{-1}-degree, n-degree) coefficient - matrix. The caller must allocate res to reslen, the number of - theta-coefficients (= 1 + max n-degree of op). */ int _gr_ore_poly_backshift_to_euler_univar(gr_ptr res, slong reslen, gr_srcptr op, slong len, gr_ctx_t ctx) { @@ -25,7 +20,7 @@ _gr_ore_poly_backshift_to_euler_univar(gr_ptr res, slong reslen, gr_srcptr op, s gr_ctx_struct * sctx; slong bsz = ctx->sizeof_elem, ssz; slong i, k; - gr_ptr c, cols; + gr_ptr _k, rcoeffs; if (ctx->which_ring != GR_CTX_GR_POLY) return GR_UNABLE; @@ -33,15 +28,15 @@ _gr_ore_poly_backshift_to_euler_univar(gr_ptr res, slong reslen, gr_srcptr op, s sctx = POLYNOMIAL_ELEM_CTX(ctx); ssz = sctx->sizeof_elem; - GR_TMP_INIT(c, sctx); - GR_TMP_INIT_VEC(cols, len, ctx); + GR_TMP_INIT(_k, sctx); + GR_TMP_INIT_VEC(rcoeffs, len, ctx); + for (k = 0; k < len; k++) { - const gr_poly_struct * ok = (const gr_poly_struct *) GR_ENTRY(op, k, bsz); - gr_poly_struct * colk = (gr_poly_struct *) GR_ENTRY(cols, k, bsz); - status |= gr_poly_set(colk, ok, sctx); - status |= gr_set_si(c, k, sctx); - status |= gr_poly_taylor_shift(colk, colk, c, sctx); + gr_poly_struct * rck = (gr_poly_struct *) GR_ENTRY(rcoeffs, k, bsz); + status |= gr_poly_set(rck, GR_ENTRY(op, k, bsz), sctx); + status |= gr_set_si(_k, k, sctx); + status |= gr_poly_taylor_shift(rck, rck, _k, sctx); } for (i = 0; i < reslen; i++) @@ -51,10 +46,10 @@ _gr_ore_poly_backshift_to_euler_univar(gr_ptr res, slong reslen, gr_srcptr op, s gr_poly_fit_length(ri, len, sctx); for (k = 0; k < len; k++) { - gr_poly_struct * colk = (gr_poly_struct *) GR_ENTRY(cols, k, bsz); + gr_poly_struct * rck = (gr_poly_struct *) GR_ENTRY(rcoeffs, k, bsz); gr_ptr dst = GR_ENTRY(ri->coeffs, k, ssz); - if (i < colk->length) - status |= gr_set(dst, GR_ENTRY(colk->coeffs, i, ssz), sctx); + if (i < rck->length) + status |= gr_set(dst, GR_ENTRY(rck->coeffs, i, ssz), sctx); else status |= gr_zero(dst, sctx); } @@ -62,8 +57,8 @@ _gr_ore_poly_backshift_to_euler_univar(gr_ptr res, slong reslen, gr_srcptr op, s _gr_poly_normalise(ri, sctx); } - GR_TMP_CLEAR_VEC(cols, len, ctx); - GR_TMP_CLEAR(c, sctx); + GR_TMP_CLEAR_VEC(rcoeffs, len, ctx); + GR_TMP_CLEAR(_k, sctx); return status; } diff --git a/src/gr_ore_poly/convert.c b/src/gr_ore_poly/convert.c index 03d96a1fcd..df7a1e17fa 100644 --- a/src/gr_ore_poly/convert.c +++ b/src/gr_ore_poly/convert.c @@ -9,10 +9,16 @@ /* generated using Claude Opus 4.8 */ -#include "gr.h" #include "gr_vec.h" #include "gr_ore_poly.h" +static int +is_shift_case(ore_algebra_t a) +{ + return a == ORE_ALGEBRA_FORWARD_SHIFT || a == ORE_ALGEBRA_BACKWARD_SHIFT + || a == ORE_ALGEBRA_FORWARD_DIFFERENCE || a == ORE_ALGEBRA_BACKWARD_DIFFERENCE; +} + int gr_ore_poly_convert(gr_ore_poly_t res, slong * p, const gr_ore_poly_t op, gr_ore_poly_ctx_t res_ctx, gr_ore_poly_ctx_t op_ctx) @@ -21,7 +27,7 @@ gr_ore_poly_convert(gr_ore_poly_t res, slong * p, const gr_ore_poly_t op, ore_algebra_t sa = GR_ORE_POLY_CTX(op_ctx)->which_algebra; ore_algebra_t da = GR_ORE_POLY_CTX(res_ctx)->which_algebra; int status = GR_SUCCESS; - slong len, sp = 0; + slong len; *p = 0; @@ -32,45 +38,35 @@ gr_ore_poly_convert(gr_ore_poly_t res, slong * p, const gr_ore_poly_t op, if (len == 0) return gr_ore_poly_zero(res, res_ctx); + gr_ore_poly_ore_data_t * od = GR_ORE_POLY_ORE_DATA(op_ctx); + slong var = (od != NULL) ? od->base_var : -1; + gr_ore_poly_fit_length(res, len, res_ctx); - if ((sa == ORE_ALGEBRA_DERIVATIVE || sa == ORE_ALGEBRA_EULER_DERIVATIVE) - && (da == ORE_ALGEBRA_DERIVATIVE || da == ORE_ALGEBRA_EULER_DERIVATIVE)) + if (sa == da) + status |= _gr_vec_set(res->coeffs, op->coeffs, len, base); + else if (sa == ORE_ALGEBRA_DERIVATIVE && da == ORE_ALGEBRA_EULER_DERIVATIVE) { - /* differential family {d/dx, theta} */ - slong var = GR_ORE_POLY_ORE_DATA(op_ctx)->base_var; - - if (sa == da) - status |= _gr_vec_set(res->coeffs, op->coeffs, len, base); - else if (sa == ORE_ALGEBRA_DERIVATIVE) /* d/dx -> theta */ - { - status |= _gr_ore_poly_ddx_to_euler(res->coeffs, op->coeffs, len, var, base); - /* res = x^(len-1) * op, i.e. x^p * res = op with p = -(len-1) */ - *p = -(len - 1); - } - else /* theta -> d/dx */ - status |= _gr_ore_poly_euler_to_ddx(res->coeffs, op->coeffs, len, var, base); + status |= _gr_ore_poly_ddx_to_euler(res->coeffs, op->coeffs, len, var, base); + *p = -(len - 1); } - else - { - /* _gr_ore_poly_shift_convert handles the shift/difference family (the - same-algebra case included) and returns GR_DOMAIN for anything that is - not a shift/difference type, which is exactly the remaining unsupported - case: a boundary-crossing pair, or an algebra such as Q_SHIFT, MAHLER - or COMMUTATIVE (whose identity we deliberately do not claim, since its - parameters could differ between the two contexts). */ - /* a custom algebra has no ore_data; shift_convert then returns GR_DOMAIN - regardless of var, so the placeholder 0 is harmless */ - gr_ore_poly_ore_data_t * od = GR_ORE_POLY_ORE_DATA(op_ctx); - slong var = (od != NULL) ? od->base_var : 0; - - status = _gr_ore_poly_shift_convert(res->coeffs, &sp, op->coeffs, len, + else if (sa == ORE_ALGEBRA_EULER_DERIVATIVE && da == ORE_ALGEBRA_DERIVATIVE) + status |= _gr_ore_poly_euler_to_ddx(res->coeffs, op->coeffs, len, var, base); + else if (sa == ORE_ALGEBRA_FORWARD_DIFFERENCE + && da == ORE_ALGEBRA_BACKWARD_DIFFERENCE) + status |= _gr_ore_poly_shift_convert_difference(res->coeffs, p, + op->coeffs, len, 1, var, + base); + else if (sa == ORE_ALGEBRA_BACKWARD_DIFFERENCE + && da == ORE_ALGEBRA_FORWARD_DIFFERENCE) + status |= _gr_ore_poly_shift_convert_difference(res->coeffs, p, + op->coeffs, len, 0, var, + base); + else if (is_shift_case(sa) && is_shift_case(da)) + status = _gr_ore_poly_shift_convert(res->coeffs, p, op->coeffs, len, sa, da, var, base); - if (status == GR_DOMAIN) - status = GR_UNABLE; - else - *p = sp; - } + else + return GR_UNABLE; if (status == GR_SUCCESS) { diff --git a/src/gr_ore_poly/ddx_to_euler.c b/src/gr_ore_poly/ddx_to_euler.c index f6f8cf61b3..bc751f7cb4 100644 --- a/src/gr_ore_poly/ddx_to_euler.c +++ b/src/gr_ore_poly/ddx_to_euler.c @@ -9,38 +9,35 @@ /* generated using Claude Opus 4.8 */ -#include "fmpz.h" #include "gr.h" #include "gr_mat.h" #include "gr_vec.h" #include "gr_ore_poly.h" -/* x^(len-1) * sum_i op[i] D^i, expanded via x^i D^i = sum_k s(i,k) theta^k with - s(i,k) the signed Stirling numbers of the first kind. Here x is the generator - of index var. */ + int _gr_ore_poly_ddx_to_euler(gr_ptr res, gr_srcptr op, slong len, slong var, gr_ctx_t ctx) { int status = GR_SUCCESS; slong sz = ctx->sizeof_elem; slong i, k; - gr_ctx_t zz; - gr_mat_t S; + gr_ctx_t ZZ; + gr_mat_t stirling; gr_vec_t gens; gr_ptr x, xp, c; if (len <= 0) return GR_SUCCESS; - gr_ctx_init_fmpz(zz); - gr_mat_init(S, len, len, zz); - /* s(i,k) = 0 for k > i, so only the lower triangle of S is used below */ - status |= gr_mat_stirling(S, 1, zz); + gr_ctx_init_fmpz(ZZ); + gr_mat_init(stirling, len, len, ZZ); + status |= gr_mat_stirling(stirling, 1, ZZ); GR_TMP_INIT3(x, xp, c, ctx); status |= _gr_vec_zero(res, len, ctx); + /* todo: gr_nth_gen? */ gr_vec_init(gens, 0, ctx); status |= gr_gens(gens, ctx); if (var >= 0 && var < gens->length) @@ -53,17 +50,19 @@ _gr_ore_poly_ddx_to_euler(gr_ptr res, gr_srcptr op, slong len, slong var, gr_ctx for (i = len - 1; i >= 0; i--) { - status |= gr_mul(c, xp, GR_ENTRY(op, i, sz), ctx); for (k = 0; k <= i; k++) - status |= gr_addmul_fmpz(GR_ENTRY(res, k, sz), c, - gr_mat_entry_srcptr(S, i, k, zz), ctx); + { + status |= gr_mul_fmpz(c, GR_ENTRY(op, i, sz), + gr_mat_entry_srcptr(stirling, i, k, ZZ), ctx); + status |= gr_addmul(GR_ENTRY(res, k, sz), c, xp, ctx); + } if (i > 0) status |= gr_mul(xp, xp, x, ctx); } GR_TMP_CLEAR3(x, xp, c, ctx); - gr_mat_clear(S, zz); - gr_ctx_clear(zz); + gr_mat_clear(stirling, ZZ); + gr_ctx_clear(ZZ); return status; } diff --git a/src/gr_ore_poly/differential_to_shift.c b/src/gr_ore_poly/differential_to_shift.c index ff736a776e..604f9b1b80 100644 --- a/src/gr_ore_poly/differential_to_shift.c +++ b/src/gr_ore_poly/differential_to_shift.c @@ -10,8 +10,6 @@ /* generated using Claude Opus 4.8 */ #include "gr.h" -#include "gr_poly.h" -#include "gr_vec.h" #include "gr_ore_poly.h" int @@ -23,27 +21,23 @@ gr_ore_poly_differential_to_shift(gr_ore_poly_t res, slong * p, const gr_ore_pol ore_algebra_t da = GR_ORE_POLY_CTX(res_ctx)->which_algebra; slong var = GR_ORE_POLY_ORE_DATA(op_ctx)->base_var; slong bsz = base->sizeof_elem; - int status = GR_SUCCESS; - slong len, rlen = 0, a = 0, sp = 0, i; + slong a, sp = 0; gr_srcptr eul; - gr_ptr cur = NULL, rec = NULL; + gr_ptr tmp0 = NULL, tmp1 = NULL; + + int status = GR_SUCCESS; *p = 0; if (GR_ORE_POLY_ELEM_CTX(res_ctx) != base || base->which_ring != GR_CTX_GR_POLY) return GR_UNABLE; - len = op->length; - if (len == 0) - return gr_ore_poly_zero(res, res_ctx); - - /* differential -> Euler (the Euler operator is op itself when sa is already - the Euler derivative) */ + slong len = op->length; if (sa == ORE_ALGEBRA_DERIVATIVE) { - GR_TMP_INIT_VEC(cur, len, base); - status |= _gr_ore_poly_ddx_to_euler(cur, op->coeffs, len, var, base); - eul = cur; + GR_TMP_INIT_VEC(tmp0, len, base); + status |= _gr_ore_poly_ddx_to_euler(tmp0, op->coeffs, len, var, base); + eul = tmp0; a = len - 1; } else if (sa == ORE_ALGEBRA_EULER_DERIVATIVE) @@ -52,30 +46,28 @@ gr_ore_poly_differential_to_shift(gr_ore_poly_t res, slong * p, const gr_ore_pol a = 0; } else - { - status = GR_DOMAIN; - goto cleanup; - } + return GR_UNABLE; - for (i = 0; i < len; i++) + slong rlen = 0; + for (slong i = 0; i < len; i++) { const gr_poly_struct * q = (const gr_poly_struct *) GR_ENTRY(eul, i, bsz); if (q->length > rlen) rlen = q->length; } - if (rlen == 0) - { - status |= gr_ore_poly_zero(res, res_ctx); - goto cleanup; - } - - /* Euler -> backward shift, then backward shift -> destination shift type */ - GR_TMP_INIT_VEC(rec, rlen, base); - status |= _gr_ore_poly_euler_to_backshift_univar(rec, rlen, eul, len, base); gr_ore_poly_fit_length(res, rlen, res_ctx); - status |= _gr_ore_poly_shift_convert(res->coeffs, &sp, rec, rlen, - ORE_ALGEBRA_BACKWARD_SHIFT, da, var, base); + if (da == ORE_ALGEBRA_BACKWARD_SHIFT) + status |= _gr_ore_poly_euler_to_backshift_univar(res->coeffs, rlen, eul, len, base); + else + { + GR_TMP_INIT_VEC(tmp1, rlen, base); + status |= _gr_ore_poly_euler_to_backshift_univar(tmp1, rlen, eul, len, base); + /* checks that the output algebra is of an appropriate type */ + status |= _gr_ore_poly_shift_convert(res->coeffs, &sp, tmp1, rlen, + ORE_ALGEBRA_BACKWARD_SHIFT, da, var, base); + GR_TMP_CLEAR_VEC(tmp1, rlen, base); + } if (status == GR_SUCCESS) { _gr_ore_poly_set_length(res, rlen, res_ctx); @@ -83,11 +75,8 @@ gr_ore_poly_differential_to_shift(gr_ore_poly_t res, slong * p, const gr_ore_pol *p = sp + a; } -cleanup: - if (cur != NULL) - GR_TMP_CLEAR_VEC(cur, len, base); - if (rec != NULL) - GR_TMP_CLEAR_VEC(rec, rlen, base); + if (tmp0 != NULL) + GR_TMP_CLEAR_VEC(tmp0, len, base); return status; } diff --git a/src/gr_ore_poly/euler_to_backshift_univar.c b/src/gr_ore_poly/euler_to_backshift_univar.c index 555b8c04aa..5a15379187 100644 --- a/src/gr_ore_poly/euler_to_backshift_univar.c +++ b/src/gr_ore_poly/euler_to_backshift_univar.c @@ -13,13 +13,6 @@ #include "gr_poly.h" #include "gr_ore_poly.h" -/* Bridge from the Euler operator (K[x], theta = x d/dx) to the backward - shift operator (K[n], S^{-1} a(n) = a(n-1)) under the exact - isomorphism theta <-> n, x <-> S^{-1}. Writing op = sum_k x^k Q_k(theta), the - image is sum_k Q_k(n-k) S^{-k}: transpose the (theta-degree, x-degree) - coefficient matrix and Taylor-shift column k by -k. The caller must allocate - res to reslen, the number of S^{-1}-coefficients (= 1 + max x-degree of - op). */ int _gr_ore_poly_euler_to_backshift_univar(gr_ptr res, slong reslen, gr_srcptr op, slong len, gr_ctx_t ctx) { @@ -27,7 +20,7 @@ _gr_ore_poly_euler_to_backshift_univar(gr_ptr res, slong reslen, gr_srcptr op, s gr_ctx_struct * sctx; slong bsz = ctx->sizeof_elem, ssz; slong i, k; - gr_ptr c; + gr_ptr negk; if (ctx->which_ring != GR_CTX_GR_POLY) return GR_UNABLE; @@ -35,7 +28,7 @@ _gr_ore_poly_euler_to_backshift_univar(gr_ptr res, slong reslen, gr_srcptr op, s sctx = POLYNOMIAL_ELEM_CTX(ctx); ssz = sctx->sizeof_elem; - GR_TMP_INIT(c, sctx); + GR_TMP_INIT(negk, sctx); for (k = 0; k < reslen; k++) { gr_poly_struct * rk = (gr_poly_struct *) GR_ENTRY(res, k, bsz); @@ -53,10 +46,10 @@ _gr_ore_poly_euler_to_backshift_univar(gr_ptr res, slong reslen, gr_srcptr op, s _gr_poly_set_length(rk, len, sctx); _gr_poly_normalise(rk, sctx); - status |= gr_set_si(c, -k, sctx); - status |= gr_poly_taylor_shift(rk, rk, c, sctx); + status |= gr_set_si(negk, -k, sctx); + status |= gr_poly_taylor_shift(rk, rk, negk, sctx); } - GR_TMP_CLEAR(c, sctx); + GR_TMP_CLEAR(negk, sctx); return status; } diff --git a/src/gr_ore_poly/euler_to_ddx.c b/src/gr_ore_poly/euler_to_ddx.c index f7e21581a5..2237a7ed6d 100644 --- a/src/gr_ore_poly/euler_to_ddx.c +++ b/src/gr_ore_poly/euler_to_ddx.c @@ -9,36 +9,32 @@ /* generated using Claude Opus 4.8 */ -#include "fmpz.h" #include "gr.h" #include "gr_mat.h" #include "gr_vec.h" #include "gr_ore_poly.h" -/* sum_k op[k] theta^k expanded via theta^k = sum_i S(k,i) x^i D^i with S(k,i) - the Stirling numbers of the second kind. Here x is the generator of index - var. */ int _gr_ore_poly_euler_to_ddx(gr_ptr res, gr_srcptr op, slong len, slong var, gr_ctx_t ctx) { int status = GR_SUCCESS; slong sz = ctx->sizeof_elem; slong i, k; - gr_ctx_t zz; - gr_mat_t S; + gr_ctx_t ZZ; + gr_mat_t stirling; gr_vec_t gens; gr_ptr x, xp, t; - if (len <= 0) + if (len == 0) return GR_SUCCESS; - gr_ctx_init_fmpz(zz); - gr_mat_init(S, len, len, zz); - /* S(k,i) = 0 for i > k, so only the lower triangle of S is used below */ - status |= gr_mat_stirling(S, 2, zz); + gr_ctx_init_fmpz(ZZ); + gr_mat_init(stirling, len, len, ZZ); + status |= gr_mat_stirling(stirling, 2, ZZ); GR_TMP_INIT3(x, xp, t, ctx); + /* todo: gr_nth_gen? */ gr_vec_init(gens, 0, ctx); status |= gr_gens(gens, ctx); if (var >= 0 && var < gens->length) @@ -49,19 +45,21 @@ _gr_ore_poly_euler_to_ddx(gr_ptr res, gr_srcptr op, slong len, slong var, gr_ctx status |= gr_one(xp, ctx); + /* theta^k = sum_i S(k,i) x^i D^i, so res[i] = (sum_k S(k, i) op[k]) x^i */ for (i = 0; i < len; i++) { status |= gr_zero(t, ctx); for (k = i; k < len; k++) status |= gr_addmul_fmpz(t, GR_ENTRY(op, k, sz), - gr_mat_entry_srcptr(S, k, i, zz), ctx); + gr_mat_entry_srcptr(stirling, k, i, ZZ), + ctx); status |= gr_mul(GR_ENTRY(res, i, sz), xp, t, ctx); status |= gr_mul(xp, xp, x, ctx); } GR_TMP_CLEAR3(x, xp, t, ctx); - gr_mat_clear(S, zz); - gr_ctx_clear(zz); + gr_mat_clear(stirling, ZZ); + gr_ctx_clear(ZZ); return status; } diff --git a/src/gr_ore_poly/shift_convert.c b/src/gr_ore_poly/shift_convert.c index ed2f3032c3..33c2ab9237 100644 --- a/src/gr_ore_poly/shift_convert.c +++ b/src/gr_ore_poly/shift_convert.c @@ -16,61 +16,48 @@ #include "gr_ore_poly.h" static int -shift_alg_is_shift(ore_algebra_t a) +is_shift_case(ore_algebra_t a) { return a == ORE_ALGEBRA_FORWARD_SHIFT || a == ORE_ALGEBRA_BACKWARD_SHIFT || a == ORE_ALGEBRA_FORWARD_DIFFERENCE || a == ORE_ALGEBRA_BACKWARD_DIFFERENCE; } static int -shift_alg_is_backward(ore_algebra_t a) +is_backward(ore_algebra_t a) { return a == ORE_ALGEBRA_BACKWARD_SHIFT || a == ORE_ALGEBRA_BACKWARD_DIFFERENCE; } static int -shift_alg_is_difference(ore_algebra_t a) +is_difference(ore_algebra_t a) { return a == ORE_ALGEBRA_FORWARD_DIFFERENCE || a == ORE_ALGEBRA_BACKWARD_DIFFERENCE; } -/* In-place change of basis between a shift generator and the associated - finite-difference generator, computing - vec[m] <- sum_{k >= m} w(k, m) * binomial(k, m) * vec[k], - where the sign w is 1 for sign == 0, (-1)^(k-m) for sign == +1, and (-1)^m - for sign == -1. The binomial coefficients along each row are built up - incrementally, as in _arb_poly_binomial_transform_basecase. The update reads - only entries of index >= m, so writing vec[m] in place is safe. */ +/* vec[m] <- sum_{k >= m} w(k, m) * binomial(k, m) * vec[k], where the sign w is + 1 for sign == 0, (-1)^(k-m) for sign == +1, and (-1)^m for sign == -1 */ static int binomial_transform(gr_ptr vec, slong len, int sign, gr_ctx_t ctx) { int status = GR_SUCCESS; slong sz = ctx->sizeof_elem, m, k; - gr_ptr acc; fmpz_t t; - GR_TMP_INIT(acc, ctx); fmpz_init(t); for (m = 0; m < len; m++) { - status |= gr_zero(acc, ctx); - for (k = m; k < len; k++) + fmpz_set_si(t, (sign == -1 && (m & 1)) ? -1 : 1); + status |= gr_mul_fmpz(GR_ENTRY(vec, m, sz), GR_ENTRY(vec, m, sz), t, ctx); + for (k = m + 1; k < len; k++) { - if (k == m) - fmpz_set_si(t, (sign == -1 && (m & 1)) ? -1 : 1); - else - { - fmpz_mul_si(t, t, (sign == 1) ? -k : k); - fmpz_divexact_ui(t, t, (ulong) (k - m)); - } - status |= gr_addmul_fmpz(acc, GR_ENTRY(vec, k, sz), t, ctx); + fmpz_mul_si(t, t, (sign == 1) ? -k : k); + fmpz_divexact_ui(t, t, (ulong) (k - m)); + status |= gr_addmul_fmpz(GR_ENTRY(vec, m, sz), GR_ENTRY(vec, k, sz), t, ctx); } - status |= gr_set(GR_ENTRY(vec, m, sz), acc, ctx); } fmpz_clear(t); - GR_TMP_CLEAR(acc, ctx); return status; } @@ -84,7 +71,6 @@ reverse(gr_ptr vec, slong len, gr_ctx_t ctx) gr_swap(GR_ENTRY(vec, k, sz), GR_ENTRY(vec, len - 1 - k, sz), ctx); } -/* Taylor-shift every coefficient by p; requires a polynomial base ring. */ static int taylor_shift_all(gr_ptr vec, slong len, slong p, gr_ctx_t ctx) { @@ -115,66 +101,78 @@ _gr_ore_poly_shift_convert(gr_ptr res, slong * p, gr_srcptr op, slong len, *p = 0; - if (!shift_alg_is_shift(src_alg) || !shift_alg_is_shift(dst_alg)) - return GR_DOMAIN; + if (!is_shift_case(src_alg) || !is_shift_case(dst_alg)) + return GR_UNABLE; if (len <= 0) return GR_SUCCESS; - /* The implementation is limited to a univariate polynomial base ring, whose - only generator has index 0. */ - if (var != 0) - return GR_UNABLE; + status |= _gr_vec_set(res, op, len, ctx); - if (src_alg == dst_alg) - return _gr_vec_set(res, op, len, ctx); + if (len == 1 || src_alg == dst_alg) + return status; - src_back = shift_alg_is_backward(src_alg); - dst_back = shift_alg_is_backward(dst_alg); + src_back = is_backward(src_alg); + dst_back = is_backward(dst_alg); crossing = (src_back != dst_back); - /* res is computed as S^s * op, so the reported p (satisfying S^p * res = op) - is -s. */ s = dst_back ? -(len - 1) : (len - 1); - /* Crossing between the forward side S, S - 1 and the backward side S^{-1}, - 1 - S^{-1} conjugates by a power of S, realized as a Taylor shift of the - coefficients; this needs a polynomial base ring once len >= 2. */ - if (crossing && len >= 2 && ctx->which_ring != GR_CTX_GR_POLY) - return GR_UNABLE; - - status |= _gr_vec_set(res, op, len, ctx); - - /* The two difference generators sit on opposite sides, so converting between - them always crosses. Because S - 1 = S (1 - S^{-1}), a single binomial - transform suffices: reversing the coefficients turns the combination into - the standard one handled by binomial_transform. */ - if (shift_alg_is_difference(src_alg) && shift_alg_is_difference(dst_alg)) + if (crossing && len >= 2) /* cases requiring a Taylor shift */ { - reverse(res, len, ctx); - status |= binomial_transform(res, len, dst_back ? 1 : 0, ctx); - reverse(res, len, ctx); - if (s != 0) - status |= taylor_shift_all(res, len, s, ctx); - *p = -s; - return status; + if (ctx->which_ring != GR_CTX_GR_POLY) + return GR_UNABLE; + else + FLINT_ASSERT(var == 0); } - /* General pipeline: source difference -> its shift, cross sides if needed, - then shift -> destination difference. At most one binomial transform runs - (the difference <-> difference case above is the only one needing two). */ - if (shift_alg_is_difference(src_alg)) + if (is_difference(src_alg)) status |= binomial_transform(res, len, src_back ? -1 : 1, ctx); if (crossing) { - if (s != 0) - status |= taylor_shift_all(res, len, s, ctx); + status |= taylor_shift_all(res, len, s, ctx); reverse(res, len, ctx); *p = -s; } - if (shift_alg_is_difference(dst_alg)) + if (is_difference(dst_alg)) status |= binomial_transform(res, len, dst_back ? -1 : 0, ctx); return status; } + + +int +_gr_ore_poly_shift_convert_difference(gr_ptr res, slong * p, gr_srcptr op, + slong len, int to_backward, slong var, + gr_ctx_t ctx) +{ + int status = GR_SUCCESS; + + *p = 0; + + if (len <= 0) + return GR_SUCCESS; + + status |= _gr_vec_set(res, op, len, ctx); + + if (len == 1) + return status; + + if (len >= 2) + { + if (ctx->which_ring != GR_CTX_GR_POLY) + return GR_UNABLE; + else + FLINT_ASSERT(var == 0); + } + + *p = to_backward ? (len - 1) : -(len - 1); + + reverse(res, len, ctx); + status |= binomial_transform(res, len, to_backward ? 1 : 0, ctx); + reverse(res, len, ctx); + status |= taylor_shift_all(res, len, -*p, ctx); + + return status; +} diff --git a/src/gr_ore_poly/shift_to_differential.c b/src/gr_ore_poly/shift_to_differential.c index 35413ff6cc..3dd3ea12c5 100644 --- a/src/gr_ore_poly/shift_to_differential.c +++ b/src/gr_ore_poly/shift_to_differential.c @@ -10,45 +10,50 @@ /* generated using Claude Opus 4.8 */ #include "gr.h" -#include "gr_poly.h" -#include "gr_vec.h" #include "gr_ore_poly.h" int -gr_ore_poly_shift_to_differential(gr_ore_poly_t res, slong * p, const gr_ore_poly_t op, - gr_ore_poly_ctx_t res_ctx, gr_ore_poly_ctx_t op_ctx) +gr_ore_poly_shift_to_differential(gr_ore_poly_t res, slong * p, + const gr_ore_poly_t op, + gr_ore_poly_ctx_t res_ctx, + gr_ore_poly_ctx_t op_ctx) { - gr_ctx_struct * base = GR_ORE_POLY_ELEM_CTX(op_ctx); - ore_algebra_t sa = GR_ORE_POLY_CTX(op_ctx)->which_algebra; - ore_algebra_t da = GR_ORE_POLY_CTX(res_ctx)->which_algebra; - slong var = GR_ORE_POLY_ORE_DATA(res_ctx)->base_var; - slong bsz = base->sizeof_elem; - int status = GR_SUCCESS; - slong len, elen = 0, sp = 0, k; - gr_ptr rec = NULL, eul = NULL; + gr_ore_poly_ctx_t rec_ctx, eul_ctx; + gr_ore_poly_t rec, eul; - *p = 0; + int status = GR_SUCCESS; + gr_ctx_struct * base = GR_ORE_POLY_ELEM_CTX(op_ctx); if (GR_ORE_POLY_ELEM_CTX(res_ctx) != base || base->which_ring != GR_CTX_GR_POLY) return GR_UNABLE; + ore_algebra_t da = GR_ORE_POLY_CTX(res_ctx)->which_algebra; if (da != ORE_ALGEBRA_DERIVATIVE && da != ORE_ALGEBRA_EULER_DERIVATIVE) - return GR_DOMAIN; + return GR_UNABLE; - len = op->length; - if (len == 0) - return gr_ore_poly_zero(res, res_ctx); + slong len = op->length; + slong var = GR_ORE_POLY_ORE_DATA(res_ctx)->base_var; - /* source shift type -> backward shift (the algebra is validated here) */ - GR_TMP_INIT_VEC(rec, len, base); - status |= _gr_ore_poly_shift_convert(rec, &sp, op->coeffs, len, - sa, ORE_ALGEBRA_BACKWARD_SHIFT, var, base); - if (status != GR_SUCCESS) + gr_ore_poly_ctx_init(rec_ctx, base, var, ORE_ALGEBRA_BACKWARD_SHIFT); + gr_ore_poly_ctx_init(eul_ctx, base, var, ORE_ALGEBRA_EULER_DERIVATIVE); + + gr_ore_poly_init(rec, rec_ctx); + gr_ore_poly_init(eul, eul_ctx); + + /* this is where the input algebra type is validated */ + status |= gr_ore_poly_convert(rec, p, op, rec_ctx, op_ctx); + *p = -*p; + + if (len == 0) + { + status |= gr_ore_poly_zero(res, res_ctx); goto cleanup; + } - for (k = 0; k < len; k++) + slong elen = 0; + for (slong k = 0; k < len; k++) { - const gr_poly_struct * q = (const gr_poly_struct *) GR_ENTRY(rec, k, bsz); + const gr_poly_struct * q = GR_ENTRY(rec->coeffs, k, base->sizeof_elem); if (q->length > elen) elen = q->length; } @@ -58,28 +63,24 @@ gr_ore_poly_shift_to_differential(gr_ore_poly_t res, slong * p, const gr_ore_pol goto cleanup; } - /* backward shift -> Euler -> destination differential type. For the Euler - destination, backshift_to_euler_univar writes the result in place. */ gr_ore_poly_fit_length(res, elen, res_ctx); if (da == ORE_ALGEBRA_EULER_DERIVATIVE) - { - status |= _gr_ore_poly_backshift_to_euler_univar(res->coeffs, elen, rec, len, base); - } + status |= _gr_ore_poly_backshift_to_euler_univar(res->coeffs, elen, rec->coeffs, len, base); else { - GR_TMP_INIT_VEC(eul, elen, base); - status |= _gr_ore_poly_backshift_to_euler_univar(eul, elen, rec, len, base); - status |= _gr_ore_poly_euler_to_ddx(res->coeffs, eul, elen, var, base); + gr_ore_poly_fit_length(eul, elen, res_ctx); + status |= _gr_ore_poly_backshift_to_euler_univar(eul->coeffs, elen, rec->coeffs, len, base); + status |= _gr_ore_poly_euler_to_ddx(res->coeffs, eul->coeffs, elen, var, base); } _gr_ore_poly_set_length(res, elen, res_ctx); _gr_ore_poly_normalise(res, res_ctx); - *p = -sp; cleanup: - if (rec != NULL) - GR_TMP_CLEAR_VEC(rec, len, base); - if (eul != NULL) - GR_TMP_CLEAR_VEC(eul, elen, base); + + gr_ore_poly_clear(rec, rec_ctx); + gr_ore_poly_clear(eul, eul_ctx); + gr_ore_poly_ctx_clear(rec_ctx); + gr_ore_poly_ctx_clear(eul_ctx); return status; } diff --git a/src/gr_ore_poly/test/main.c b/src/gr_ore_poly/test/main.c index edbf9648a3..b13ee8cc9b 100644 --- a/src/gr_ore_poly/test/main.c +++ b/src/gr_ore_poly/test/main.c @@ -20,7 +20,7 @@ #include "t-ddx_to_euler.c" #include "t-euler_to_ddx.c" #include "t-shift_convert.c" -#include "t-differential_shift.c" +#include "t-differential_to_shift.c" #include "t-convert.c" /* Array of test functions ***************************************************/ diff --git a/src/gr_ore_poly/test/t-apply.c b/src/gr_ore_poly/test/t-apply.c index f498e55964..af34738dd5 100644 --- a/src/gr_ore_poly/test/t-apply.c +++ b/src/gr_ore_poly/test/t-apply.c @@ -26,8 +26,7 @@ check_gen_action(gr_srcptr x, gr_srcptr expected, gr_ore_poly_ctx_t octx, gr_ctx if (status != GR_SUCCESS || gr_equal(got, expected, cctx) != T_TRUE) { - flint_printf("FAIL: got %{gr}, expected %{gr}, in %{gr_ctx}\n", - got, cctx, expected, cctx, cctx); + flint_printf("FAIL: D(gen)\n\noctx = %{gr_ctx}\n", octx); flint_abort(); } @@ -35,106 +34,102 @@ check_gen_action(gr_srcptr x, gr_srcptr expected, gr_ore_poly_ctx_t octx, gr_ctx gr_ore_poly_clear(D, octx); } -TEST_GR_FUNCTION_START(gr_ore_poly_apply, state, count_success, count_domain, count_unable) +static void +check_gen_actions(flint_rand_t state) { - { - int status = GR_SUCCESS; - gr_ctx_t zctx, cctx; - gr_ore_poly_ctx_t octx; - gr_ptr x, q, expected; - - gr_ctx_init_fmpz(zctx); - gr_ctx_init_gr_poly(cctx, zctx); - x = gr_heap_init(cctx); - q = gr_heap_init(cctx); - expected = gr_heap_init(cctx); - status |= gr_gen(x, cctx); - - gr_ore_poly_ctx_init(octx, cctx, 0, ORE_ALGEBRA_DERIVATIVE); - status |= gr_one(expected, cctx); - check_gen_action(x, expected, octx, cctx); - gr_ore_poly_ctx_clear(octx); - - gr_ore_poly_ctx_init(octx, cctx, 0, ORE_ALGEBRA_EULER_DERIVATIVE); - status |= gr_set(expected, x, cctx); - check_gen_action(x, expected, octx, cctx); - gr_ore_poly_ctx_clear(octx); - - gr_ore_poly_ctx_init(octx, cctx, 0, ORE_ALGEBRA_FORWARD_SHIFT); - status |= gr_add_si(expected, x, 1, cctx); - check_gen_action(x, expected, octx, cctx); - gr_ore_poly_ctx_clear(octx); + int status = GR_SUCCESS; + gr_ctx_t zctx, cctx; + gr_ore_poly_ctx_t octx; + gr_ptr x, q, expected; + + gr_ctx_init_fmpz(zctx); + gr_ctx_init_gr_poly(cctx, zctx); + x = gr_heap_init(cctx); + q = gr_heap_init(cctx); + expected = gr_heap_init(cctx); + status |= gr_gen(x, cctx); + + gr_ore_poly_ctx_init(octx, cctx, 0, ORE_ALGEBRA_DERIVATIVE); + status |= gr_one(expected, cctx); + check_gen_action(x, expected, octx, cctx); + gr_ore_poly_ctx_clear(octx); + + gr_ore_poly_ctx_init(octx, cctx, 0, ORE_ALGEBRA_EULER_DERIVATIVE); + status |= gr_set(expected, x, cctx); + check_gen_action(x, expected, octx, cctx); + gr_ore_poly_ctx_clear(octx); + + gr_ore_poly_ctx_init(octx, cctx, 0, ORE_ALGEBRA_FORWARD_SHIFT); + status |= gr_add_si(expected, x, 1, cctx); + check_gen_action(x, expected, octx, cctx); + gr_ore_poly_ctx_clear(octx); + + gr_ore_poly_ctx_init(octx, cctx, 0, ORE_ALGEBRA_FORWARD_DIFFERENCE); + status |= gr_one(expected, cctx); + check_gen_action(x, expected, octx, cctx); + gr_ore_poly_ctx_clear(octx); + + gr_ore_poly_ctx_init(octx, cctx, 0, ORE_ALGEBRA_BACKWARD_SHIFT); + status |= gr_add_si(expected, x, -1, cctx); + check_gen_action(x, expected, octx, cctx); + gr_ore_poly_ctx_clear(octx); + + gr_ore_poly_ctx_init(octx, cctx, 0, ORE_ALGEBRA_BACKWARD_DIFFERENCE); + status |= gr_one(expected, cctx); + check_gen_action(x, expected, octx, cctx); + gr_ore_poly_ctx_clear(octx); + + status |= gr_set_si(q, 3, cctx); + status |= gr_ore_poly_ctx_init_q_shift(octx, cctx, 0, q); + status |= gr_mul(expected, q, x, cctx); + check_gen_action(x, expected, octx, cctx); + gr_ore_poly_ctx_clear(octx); + + status |= gr_ore_poly_ctx_init_mahler(octx, cctx, 0, 3); + status |= gr_pow_ui(expected, x, 3, cctx); + check_gen_action(x, expected, octx, cctx); + gr_ore_poly_ctx_clear(octx); + + gr_heap_clear(x, cctx); + gr_heap_clear(q, cctx); + gr_heap_clear(expected, cctx); + gr_ctx_clear(cctx); + gr_ctx_clear(zctx); - gr_ore_poly_ctx_init(octx, cctx, 0, ORE_ALGEBRA_FORWARD_DIFFERENCE); - status |= gr_one(expected, cctx); - check_gen_action(x, expected, octx, cctx); - gr_ore_poly_ctx_clear(octx); + { + fmpz_t p; + gr_ctx_t fctx; + gr_ptr fx, fexpected; - gr_ore_poly_ctx_init(octx, cctx, 0, ORE_ALGEBRA_BACKWARD_SHIFT); - status |= gr_add_si(expected, x, -1, cctx); - check_gen_action(x, expected, octx, cctx); - gr_ore_poly_ctx_clear(octx); + fmpz_init_set_ui(p, 17); + gr_ctx_init_fq(fctx, p, 3, NULL); + fmpz_clear(p); - gr_ore_poly_ctx_init(octx, cctx, 0, ORE_ALGEBRA_BACKWARD_DIFFERENCE); - status |= gr_one(expected, cctx); - check_gen_action(x, expected, octx, cctx); - gr_ore_poly_ctx_clear(octx); + fx = gr_heap_init(fctx); + fexpected = gr_heap_init(fctx); - status |= gr_set_si(q, 3, cctx); - status |= gr_ore_poly_ctx_init_q_shift(octx, cctx, 0, q); - status |= gr_mul(expected, q, x, cctx); - check_gen_action(x, expected, octx, cctx); - gr_ore_poly_ctx_clear(octx); + status |= gr_gen(fx, fctx); + status |= gr_pow_ui(fexpected, fx, 17, fctx); - status |= gr_ore_poly_ctx_init_mahler(octx, cctx, 0, 3); - status |= gr_pow_ui(expected, x, 3, cctx); - check_gen_action(x, expected, octx, cctx); + gr_ore_poly_ctx_init(octx, fctx, 0, ORE_ALGEBRA_FROBENIUS); + check_gen_action(fx, fexpected, octx, fctx); gr_ore_poly_ctx_clear(octx); - if (status != GR_SUCCESS) - { - flint_printf("FAIL: unexpected failure in %{gr_ctx}\n", cctx, cctx); - flint_abort(); - } - - gr_heap_clear(x, cctx); - gr_heap_clear(q, cctx); - gr_heap_clear(expected, cctx); - gr_ctx_clear(cctx); - gr_ctx_clear(zctx); + gr_heap_clear(fx, fctx); + gr_heap_clear(fexpected, fctx); + gr_ctx_clear(fctx); + } - { - fmpz_t p; - gr_ctx_t fctx; - gr_ptr fx, fexpected; - - fmpz_init_set_ui(p, 17); - gr_ctx_init_fq(fctx, p, 3, NULL); - fmpz_clear(p); - - fx = gr_heap_init(fctx); - fexpected = gr_heap_init(fctx); - - status = GR_SUCCESS; - status |= gr_gen(fx, fctx); - status |= gr_pow_ui(fexpected, fx, 17, fctx); - - gr_ore_poly_ctx_init(octx, fctx, 0, ORE_ALGEBRA_FROBENIUS); - check_gen_action(fx, fexpected, octx, fctx); - gr_ore_poly_ctx_clear(octx); - - if (status != GR_SUCCESS) - { - flint_printf("FAIL: unexpected failure in %{gr_ctx}\n", fctx, fctx); - flint_abort(); - } - - gr_heap_clear(fx, fctx); - gr_heap_clear(fexpected, fctx); - gr_ctx_clear(fctx); - } + if (status != GR_SUCCESS) + { + flint_printf("FAIL: unexpected failure\n"); + flint_abort(); } +} +TEST_GR_FUNCTION_START(gr_ore_poly_apply, state, count_success, count_domain, count_unable) +{ + check_gen_actions(state); for (slong iter = 0; iter < 1000 * flint_test_multiplier(); iter++) { @@ -144,8 +139,6 @@ TEST_GR_FUNCTION_START(gr_ore_poly_apply, state, count_success, count_domain, co gr_ore_poly_ctx_init_randtest2(cctx, ctx, state); - /* Repeated application of D blows up the degree for these, so keep the - operators short. */ if (GR_ORE_POLY_CTX(ctx)->which_algebra == ORE_ALGEBRA_MAHLER || GR_ORE_POLY_CTX(ctx)->which_algebra == ORE_ALGEBRA_Q_SHIFT) maxlen = 2; @@ -175,11 +168,27 @@ TEST_GR_FUNCTION_START(gr_ore_poly_apply, state, count_success, count_domain, co status |= gr_ore_poly_randtest(P, state, 1 + n_randint(state, maxlen), ctx); status |= gr_ore_poly_randtest(Q, state, 1 + n_randint(state, maxlen), ctx); status |= gr_ore_poly_gen(D, ctx); - status |= gr_randtest(f, state, cctx); + if (n_randint(state, 8)) + status |= gr_randtest_not_zero(f, state, cctx); + else + status |= gr_randtest(f, state, cctx); status |= gr_randtest(g, state, cctx); - status |= gr_randtest(d1, state, cctx); - /* D(f) = sigma(f)*d1 + delta(f) */ + status |= gr_one(one, cctx); + status |= gr_ore_poly_apply(c, D, one, ctx); + status |= gr_ore_poly_apply(lhs, P, f, ctx); + status |= gr_ore_poly_apply_custom(rhs, P, f, c, ctx); + if (status == GR_SUCCESS && gr_equal(lhs, rhs, cctx) == T_FALSE) + { + flint_printf("FAIL: apply(P, f) = apply_custom(P, f, D(1))\n"); + flint_abort(); + } + + if (status == GR_SUCCESS && n_randint(state, 2)) + status |= gr_set(d1, c, cctx); + else + status |= gr_randtest(d1, state, cctx); + status |= gr_ore_poly_apply_custom(u, D, f, d1, ctx); status |= gr_ore_poly_sigma_delta(sf, df, f, ctx); status |= gr_mul(sf, sf, d1, cctx); @@ -190,7 +199,6 @@ TEST_GR_FUNCTION_START(gr_ore_poly_apply, state, count_success, count_domain, co flint_abort(); } - /* (P*Q)(f) = P(Q(f)) */ status |= gr_ore_poly_mul(PQ, P, Q, ctx); status |= gr_ore_poly_apply_custom(lhs, PQ, f, d1, ctx); status |= gr_ore_poly_apply_custom(u, Q, f, d1, ctx); @@ -201,7 +209,6 @@ TEST_GR_FUNCTION_START(gr_ore_poly_apply, state, count_success, count_domain, co flint_abort(); } - /* (P+Q)(f) = P(f) + Q(f) */ status |= gr_ore_poly_add(PpQ, P, Q, ctx); status |= gr_ore_poly_apply_custom(lhs, PpQ, f, d1, ctx); status |= gr_ore_poly_apply_custom(u, P, f, d1, ctx); @@ -213,7 +220,6 @@ TEST_GR_FUNCTION_START(gr_ore_poly_apply, state, count_success, count_domain, co flint_abort(); } - /* P(f+g) = P(f) + P(g) */ status |= gr_add(w, f, g, cctx); status |= gr_ore_poly_apply_custom(lhs, P, w, d1, ctx); status |= gr_ore_poly_apply_custom(u, P, f, d1, ctx); @@ -225,7 +231,6 @@ TEST_GR_FUNCTION_START(gr_ore_poly_apply, state, count_success, count_domain, co flint_abort(); } - /* aliasing res == f */ status |= gr_ore_poly_apply_custom(u, P, f, d1, ctx); status |= gr_set(v, f, cctx); status |= gr_ore_poly_apply_custom(v, P, v, d1, ctx); @@ -235,7 +240,6 @@ TEST_GR_FUNCTION_START(gr_ore_poly_apply, state, count_success, count_domain, co flint_abort(); } - /* aliasing res == d1 */ status |= gr_ore_poly_apply_custom(u, P, f, d1, ctx); status |= gr_set(v, d1, cctx); status |= gr_ore_poly_apply_custom(v, P, f, v, ctx); @@ -245,17 +249,6 @@ TEST_GR_FUNCTION_START(gr_ore_poly_apply, state, count_success, count_domain, co flint_abort(); } - /* standard apply matches custom with c = D(1) = apply(D, 1) */ - status |= gr_one(one, cctx); - status |= gr_ore_poly_apply(c, D, one, ctx); - status |= gr_ore_poly_apply(lhs, P, f, ctx); - status |= gr_ore_poly_apply_custom(rhs, P, f, c, ctx); - if (status == GR_SUCCESS && gr_equal(lhs, rhs, cctx) == T_FALSE) - { - flint_printf("FAIL: apply(P, f) = apply_custom(P, f, D(1))\n"); - flint_abort(); - } - count_success += (status == GR_SUCCESS); count_domain += ((status & GR_DOMAIN) != 0); count_unable += ((status & GR_UNABLE) != 0); diff --git a/src/gr_ore_poly/test/t-ddx_to_euler.c b/src/gr_ore_poly/test/t-ddx_to_euler.c index 04b581f377..b1e6fec31e 100644 --- a/src/gr_ore_poly/test/t-ddx_to_euler.c +++ b/src/gr_ore_poly/test/t-ddx_to_euler.c @@ -24,15 +24,12 @@ TEST_GR_FUNCTION_START(gr_ore_poly_ddx_to_euler, state, count_success, count_dom slong len, i, j, sz, ngens, var; int status = GR_SUCCESS; - /* occasionally use a multivariate base ring so that the conversions are - exercised with a generator other than the default one */ switch (n_randint(state, 8)) { case 0: gr_ctx_init_random(cctx, state); break; case 1: - case 2: gr_ctx_init_random_mpoly(cctx, state); break; default: @@ -55,6 +52,14 @@ TEST_GR_FUNCTION_START(gr_ore_poly_ddx_to_euler, state, count_success, count_dom gr_ore_poly_init(P_d2, ctx_d); gr_ore_poly_init(scaled_P, ctx_d); + gr_ptr xpow = gr_heap_init(cctx); + gr_ptr f = gr_heap_init(cctx); + gr_ptr g_d = gr_heap_init(cctx); + gr_ptr g_e = gr_heap_init(cctx); + gr_ptr scaled = gr_heap_init(cctx); + + /* main run */ + status |= gr_ore_poly_randtest(P_d, state, 1 + n_randint(state, 5), ctx_d); len = P_d->length; @@ -63,18 +68,22 @@ TEST_GR_FUNCTION_START(gr_ore_poly_ddx_to_euler, state, count_success, count_dom _gr_ore_poly_set_length(P_e, len, ctx_e); _gr_ore_poly_normalise(P_e, ctx_e); - gr_ore_poly_fit_length(P_d2, len, ctx_d); - status |= _gr_ore_poly_euler_to_ddx(P_d2->coeffs, P_e->coeffs, len, var, cctx); - _gr_ore_poly_set_length(P_d2, len, ctx_d); - _gr_ore_poly_normalise(P_d2, ctx_d); + int expect_success = (var == 0 && cctx->which_ring == GR_CTX_GR_POLY + && (POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_FMPZ + || POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_CC_ACB + || POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_NMOD)); + if (expect_success && status != GR_SUCCESS) + { + flint_printf("FAIL: unexpected failure\n"); + flint_abort(); + } - gr_ptr xpow = gr_heap_init(cctx); - gr_ptr f = gr_heap_init(cctx); - gr_ptr g_d = gr_heap_init(cctx); - gr_ptr g_e = gr_heap_init(cctx); - gr_ptr scaled = gr_heap_init(cctx); + count_success += (status == GR_SUCCESS); + count_domain += ((status & GR_DOMAIN) != 0); + count_unable += ((status & GR_UNABLE) != 0); + + /* correcting factor */ - /* x is the generator of index var */ gr_vec_init(gens, 0, cctx); status |= gr_gens(gens, cctx); if (var < gens->length) @@ -84,7 +93,13 @@ TEST_GR_FUNCTION_START(gr_ore_poly_ddx_to_euler, state, count_success, count_dom gr_vec_clear(gens, cctx); status |= gr_pow_ui(xpow, xpow, (len >= 1) ? (ulong) (len - 1) : 0, cctx); - /* round trip euler_to_ddx(ddx_to_euler(P_d)) = x^(len-1) P_d */ + /* test round trip */ + + gr_ore_poly_fit_length(P_d2, len, ctx_d); + status |= _gr_ore_poly_euler_to_ddx(P_d2->coeffs, P_e->coeffs, len, var, cctx); + _gr_ore_poly_set_length(P_d2, len, ctx_d); + _gr_ore_poly_normalise(P_d2, ctx_d); + gr_ore_poly_fit_length(scaled_P, len, ctx_d); for (i = 0; i < len; i++) status |= gr_mul(GR_ENTRY(scaled_P->coeffs, i, sz), xpow, GR_ENTRY(P_d->coeffs, i, sz), cctx); @@ -93,40 +108,26 @@ TEST_GR_FUNCTION_START(gr_ore_poly_ddx_to_euler, state, count_success, count_dom if (status == GR_SUCCESS && gr_ore_poly_equal(P_d2, scaled_P, ctx_d) == T_FALSE) { - flint_printf("FAIL: ddx_to_euler round trip\n"); + flint_printf("FAIL: round trip\n"); flint_abort(); } + /* test action */ + for (j = 0; j < 4; j++) { - status |= gr_randtest(f, state, cctx); + status |= gr_randtest_not_zero(f, state, cctx); status |= gr_ore_poly_apply(g_d, P_d, f, ctx_d); status |= gr_ore_poly_apply(g_e, P_e, f, ctx_e); - /* P_e(f) = x^(len-1) P_d(f) */ status |= gr_mul(scaled, xpow, g_d, cctx); if (status == GR_SUCCESS && gr_equal(g_e, scaled, cctx) == T_FALSE) { - flint_printf("FAIL: ddx_to_euler application identity\n"); + flint_printf("FAIL: application identity\n"); flint_abort(); } } - /* these representative base rings never fail for the implemented var 0 */ - int expect_success = (var == 0 && cctx->which_ring == GR_CTX_GR_POLY - && (POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_FMPZ - || POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_CC_ACB - || POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_NMOD)); - if (expect_success && status != GR_SUCCESS) - { - flint_printf("FAIL: ddx_to_euler unexpected failure\n"); - flint_abort(); - } - - count_success += (status == GR_SUCCESS); - count_domain += ((status & GR_DOMAIN) != 0); - count_unable += ((status & GR_UNABLE) != 0); - gr_heap_clear(xpow, cctx); gr_heap_clear(f, cctx); gr_heap_clear(g_d, cctx); diff --git a/src/gr_ore_poly/test/t-differential_shift.c b/src/gr_ore_poly/test/t-differential_to_shift.c similarity index 61% rename from src/gr_ore_poly/test/t-differential_shift.c rename to src/gr_ore_poly/test/t-differential_to_shift.c index 68ac9411c1..75171ab22d 100644 --- a/src/gr_ore_poly/test/t-differential_shift.c +++ b/src/gr_ore_poly/test/t-differential_to_shift.c @@ -54,29 +54,29 @@ check_scaled(const gr_ore_poly_t a, const gr_ore_poly_t b, slong net, return ok; } -/* Apply a forward-shift recurrence R = sum_k q_k(nu) S^k to the coefficients - (a_m = coeff_m(f)) of a polynomial: s_m = sum_k q_k(m) a_{m+k}. */ +/* Applies a forward-shift recurrence R = sum_k q_k(nu) S^k to the coefficients + of a polynomial. */ static int apply_forward_recurrence(gr_poly_t s, const gr_ore_poly_t R, const gr_poly_t f, gr_ctx_t cctx) { gr_ctx_struct * sctx = POLYNOMIAL_ELEM_CTX(cctx); slong bsz = cctx->sizeof_elem, ssz = sctx->sizeof_elem; - slong df = f->length, m, k; + slong df = f->length; int status = GR_SUCCESS; gr_ptr mval, qval, term; GR_TMP_INIT3(mval, qval, term, sctx); - gr_poly_fit_length(s, FLINT_MAX(df, 1), sctx); - for (m = 0; m < df; m++) + gr_poly_fit_length(s, df, sctx); + for (slong n = 0; n < df; n++) { - gr_ptr sm = GR_ENTRY(s->coeffs, m, ssz); + gr_ptr sm = GR_ENTRY(s->coeffs, n, ssz); status |= gr_zero(sm, sctx); - status |= gr_set_si(mval, m, sctx); - for (k = 0; k < R->length && m + k < df; k++) + status |= gr_set_si(mval, n, sctx); + for (slong k = 0; k < R->length && n + k < df; k++) { const gr_poly_struct * qk = (const gr_poly_struct *) GR_ENTRY(R->coeffs, k, bsz); status |= gr_poly_evaluate(qval, qk, mval, sctx); - status |= gr_mul(term, qval, GR_ENTRY(f->coeffs, m + k, ssz), sctx); + status |= gr_mul(term, qval, GR_ENTRY(f->coeffs, n + k, ssz), sctx); status |= gr_add(sm, sm, term, sctx); } } @@ -91,21 +91,17 @@ TEST_GR_FUNCTION_START(gr_ore_poly_differential_to_shift, state, count_success, for (slong iter = 0; iter < 1000 * flint_test_multiplier(); iter++) { gr_ctx_t cctx; - gr_ore_poly_ctx_t ctx_d, ctx_s, ctx_fwd, ctx_eul, ctx_back; + gr_ore_poly_ctx_t ctx_d, ctx_s; ore_algebra_t diff_alg, shift_alg; - slong p1, p2, pf, pL, pM, pLM, ngens, var; + slong p1, p2, ngens, var; int status = GR_SUCCESS; - /* occasionally use a multivariate base ring with a non-default - generator; the bridges are only implemented over a univariate - polynomial base and return GR_UNABLE otherwise */ switch (n_randint(state, 8)) { case 0: gr_ctx_init_random(cctx, state); break; case 1: - case 2: gr_ctx_init_random_mpoly(cctx, state); break; default: @@ -115,16 +111,13 @@ TEST_GR_FUNCTION_START(gr_ore_poly_differential_to_shift, state, count_success, ngens = 0; status |= gr_ctx_ngens(&ngens, cctx); - var = (ngens >= 1) ? (slong) n_randint(state, ngens) : 0; + var = n_randint(state, ngens); diff_alg = ds_diff_algs[n_randint(state, 2)]; shift_alg = ds_shift_algs[n_randint(state, 4)]; gr_ore_poly_ctx_init(ctx_d, cctx, var, diff_alg); gr_ore_poly_ctx_init(ctx_s, cctx, var, shift_alg); - gr_ore_poly_ctx_init(ctx_fwd, cctx, var, ORE_ALGEBRA_FORWARD_SHIFT); - gr_ore_poly_ctx_init(ctx_eul, cctx, var, ORE_ALGEBRA_EULER_DERIVATIVE); - gr_ore_poly_ctx_init(ctx_back, cctx, var, ORE_ALGEBRA_BACKWARD_SHIFT); gr_ore_poly_t op, res_s, op2; gr_ore_poly_init(op, ctx_d); @@ -133,108 +126,69 @@ TEST_GR_FUNCTION_START(gr_ore_poly_differential_to_shift, state, count_success, status |= gr_ore_poly_randtest(op, state, 1 + n_randint(state, 5), ctx_d); - /* round trip recovers x^(p1-p2) * op (the corrections are left powers) */ status |= gr_ore_poly_differential_to_shift(res_s, &p1, op, ctx_s, ctx_d); status |= gr_ore_poly_shift_to_differential(op2, &p2, res_s, ctx_d, ctx_s); + + int expect_success = (var == 0 && cctx->which_ring == GR_CTX_GR_POLY + && (POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_FMPZ + || POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_CC_ACB + || POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_NMOD)); + if (expect_success && status != GR_SUCCESS) + { + flint_printf("FAIL: differential<->shift unexpected failure\n"); + flint_abort(); + } + + count_success += (status == GR_SUCCESS); + count_domain += ((status & GR_DOMAIN) != 0); + count_unable += ((status & GR_UNABLE) != 0); + + /* round trip recovers x^(p1-p2) * op */ + if (status == GR_SUCCESS && !check_scaled(op2, op, p1 - p2, ctx_d, cctx)) { - flint_printf("FAIL: differential<->shift round trip\n"); + flint_printf("FAIL: round trip\n"); flint_abort(); } - /* series action: a forward-shift recurrence R from op satisfies + /* a forward-shift recurrence R from op satisfies (R a)_m = coeff_{m-pf}(op . f) for the coefficient sequence of f */ - if (cctx->which_ring == GR_CTX_GR_POLY) + + if (cctx->which_ring == GR_CTX_GR_POLY && shift_alg == ORE_ALGEBRA_FORWARD_SHIFT) { gr_ctx_struct * sctx = POLYNOMIAL_ELEM_CTX(cctx); - gr_ore_poly_t Rf; gr_ptr f = gr_heap_init(cctx); gr_ptr g = gr_heap_init(cctx); gr_poly_t s, eg; - gr_ore_poly_init(Rf, ctx_fwd); gr_poly_init(s, sctx); gr_poly_init(eg, sctx); - status |= gr_ore_poly_differential_to_shift(Rf, &pf, op, ctx_fwd, ctx_d); status |= gr_randtest(f, state, cctx); status |= gr_ore_poly_apply(g, op, f, ctx_d); - status |= apply_forward_recurrence(s, Rf, (gr_poly_struct *) f, cctx); - if (pf >= 0) - status |= gr_poly_shift_left(eg, (gr_poly_struct *) g, pf, sctx); + status |= apply_forward_recurrence(s, res_s, (gr_poly_struct *) f, cctx); + if (p1 >= 0) + status |= gr_poly_shift_left(eg, (gr_poly_struct *) g, p1, sctx); else - status |= gr_poly_shift_right(eg, (gr_poly_struct *) g, -pf, sctx); + status |= gr_poly_shift_right(eg, (gr_poly_struct *) g, -p1, sctx); + if (status == GR_SUCCESS && gr_poly_equal(s, eg, sctx) == T_FALSE) { - flint_printf("FAIL: differential->shift series action\n"); + flint_printf("FAIL: action\n"); flint_abort(); } - gr_ore_poly_clear(Rf, ctx_fwd); gr_poly_clear(s, sctx); gr_poly_clear(eg, sctx); gr_heap_clear(f, cctx); gr_heap_clear(g, cctx); } - /* the Euler -> backward shift bridge is an exact algebra homomorphism */ - gr_ore_poly_t L, M, LM, Lb, Mb, LMb, LbMb; - gr_ore_poly_init(L, ctx_eul); - gr_ore_poly_init(M, ctx_eul); - gr_ore_poly_init(LM, ctx_eul); - gr_ore_poly_init(Lb, ctx_back); - gr_ore_poly_init(Mb, ctx_back); - gr_ore_poly_init(LMb, ctx_back); - gr_ore_poly_init(LbMb, ctx_back); - - status |= gr_ore_poly_randtest(L, state, 1 + n_randint(state, 3), ctx_eul); - status |= gr_ore_poly_randtest(M, state, 1 + n_randint(state, 3), ctx_eul); - status |= gr_ore_poly_mul(LM, L, M, ctx_eul); - status |= gr_ore_poly_differential_to_shift(Lb, &pL, L, ctx_back, ctx_eul); - status |= gr_ore_poly_differential_to_shift(Mb, &pM, M, ctx_back, ctx_eul); - status |= gr_ore_poly_differential_to_shift(LMb, &pLM, LM, ctx_back, ctx_eul); - status |= gr_ore_poly_mul(LbMb, Lb, Mb, ctx_back); - if (status == GR_SUCCESS && (pL != 0 || pM != 0 || pLM != 0)) - { - flint_printf("FAIL: Euler -> backward shift not exact\n"); - flint_abort(); - } - if (status == GR_SUCCESS && gr_ore_poly_equal(LMb, LbMb, ctx_back) == T_FALSE) - { - flint_printf("FAIL: bridge not multiplicative\n"); - flint_abort(); - } - - /* these representative base rings never fail for the implemented var 0 */ - int expect_success = (var == 0 && cctx->which_ring == GR_CTX_GR_POLY - && (POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_FMPZ - || POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_CC_ACB - || POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_NMOD)); - if (expect_success && status != GR_SUCCESS) - { - flint_printf("FAIL: differential<->shift unexpected failure\n"); - flint_abort(); - } - - count_success += (status == GR_SUCCESS); - count_domain += ((status & GR_DOMAIN) != 0); - count_unable += ((status & GR_UNABLE) != 0); - gr_ore_poly_clear(op, ctx_d); gr_ore_poly_clear(res_s, ctx_s); gr_ore_poly_clear(op2, ctx_d); - gr_ore_poly_clear(L, ctx_eul); - gr_ore_poly_clear(M, ctx_eul); - gr_ore_poly_clear(LM, ctx_eul); - gr_ore_poly_clear(Lb, ctx_back); - gr_ore_poly_clear(Mb, ctx_back); - gr_ore_poly_clear(LMb, ctx_back); - gr_ore_poly_clear(LbMb, ctx_back); gr_ore_poly_ctx_clear(ctx_d); gr_ore_poly_ctx_clear(ctx_s); - gr_ore_poly_ctx_clear(ctx_fwd); - gr_ore_poly_ctx_clear(ctx_eul); - gr_ore_poly_ctx_clear(ctx_back); gr_ctx_clear(cctx); } diff --git a/src/gr_ore_poly/test/t-euler_to_ddx.c b/src/gr_ore_poly/test/t-euler_to_ddx.c index c5693537fe..61317be147 100644 --- a/src/gr_ore_poly/test/t-euler_to_ddx.c +++ b/src/gr_ore_poly/test/t-euler_to_ddx.c @@ -24,15 +24,12 @@ TEST_GR_FUNCTION_START(gr_ore_poly_euler_to_ddx, state, count_success, count_dom slong len, i, j, sz, ngens, var; int status = GR_SUCCESS; - /* occasionally use a multivariate base ring so that the conversions are - exercised with a generator other than the default one */ switch (n_randint(state, 8)) { case 0: gr_ctx_init_random(cctx, state); break; case 1: - case 2: gr_ctx_init_random_mpoly(cctx, state); break; default: @@ -44,7 +41,12 @@ TEST_GR_FUNCTION_START(gr_ore_poly_euler_to_ddx, state, count_success, count_dom ngens = 0; status |= gr_ctx_ngens(&ngens, cctx); - var = (ngens >= 1) ? (slong) n_randint(state, ngens) : 0; + if (ngens < 1) + { + gr_ctx_clear(cctx); + continue; + } + var = n_randint(state, ngens); gr_ore_poly_ctx_init(ctx_d, cctx, var, ORE_ALGEBRA_DERIVATIVE); gr_ore_poly_ctx_init(ctx_e, cctx, var, ORE_ALGEBRA_EULER_DERIVATIVE); @@ -55,6 +57,13 @@ TEST_GR_FUNCTION_START(gr_ore_poly_euler_to_ddx, state, count_success, count_dom gr_ore_poly_init(P_e2, ctx_e); gr_ore_poly_init(scaled_P, ctx_e); + gr_ptr xpow = gr_heap_init(cctx); + gr_ptr f = gr_heap_init(cctx); + gr_ptr g_e = gr_heap_init(cctx); + gr_ptr g_d = gr_heap_init(cctx); + + gr_vec_init(gens, 0, cctx); + status |= gr_ore_poly_randtest(P_e, state, 1 + n_randint(state, 5), ctx_e); len = P_e->length; @@ -63,27 +72,34 @@ TEST_GR_FUNCTION_START(gr_ore_poly_euler_to_ddx, state, count_success, count_dom _gr_ore_poly_set_length(P_d, len, ctx_d); _gr_ore_poly_normalise(P_d, ctx_d); + int expect_success = (cctx->which_ring == GR_CTX_GR_POLY + && (POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_FMPZ + || POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_CC_ACB + || POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_NMOD)); + if (expect_success && status != GR_SUCCESS) + { + flint_printf("FAIL: unexpected failure\n"); + flint_abort(); + } + + count_success += (status == GR_SUCCESS); + count_domain += ((status & GR_DOMAIN) != 0); + count_unable += ((status & GR_UNABLE) != 0); + + if (status != GR_SUCCESS) + goto cleanup; + + /* test round trip */ + gr_ore_poly_fit_length(P_e2, len, ctx_e); status |= _gr_ore_poly_ddx_to_euler(P_e2->coeffs, P_d->coeffs, len, var, cctx); _gr_ore_poly_set_length(P_e2, len, ctx_e); _gr_ore_poly_normalise(P_e2, ctx_e); - gr_ptr xpow = gr_heap_init(cctx); - gr_ptr f = gr_heap_init(cctx); - gr_ptr g_e = gr_heap_init(cctx); - gr_ptr g_d = gr_heap_init(cctx); - - /* x is the generator of index var */ - gr_vec_init(gens, 0, cctx); status |= gr_gens(gens, cctx); - if (var < gens->length) - status |= gr_set(xpow, gr_vec_entry_srcptr(gens, var, cctx), cctx); - else - status |= GR_UNABLE; - gr_vec_clear(gens, cctx); + status |= gr_set(xpow, gr_vec_entry_srcptr(gens, var, cctx), cctx); status |= gr_pow_ui(xpow, xpow, (len >= 1) ? (ulong) (len - 1) : 0, cctx); - /* round trip ddx_to_euler(euler_to_ddx(P_e)) = x^(len-1) P_e */ gr_ore_poly_fit_length(scaled_P, len, ctx_e); for (i = 0; i < len; i++) status |= gr_mul(GR_ENTRY(scaled_P->coeffs, i, sz), xpow, GR_ENTRY(P_e->coeffs, i, sz), cctx); @@ -92,39 +108,26 @@ TEST_GR_FUNCTION_START(gr_ore_poly_euler_to_ddx, state, count_success, count_dom if (status == GR_SUCCESS && gr_ore_poly_equal(P_e2, scaled_P, ctx_e) == T_FALSE) { - flint_printf("FAIL: euler_to_ddx round trip\n"); + flint_printf("FAIL: round trip\n"); flint_abort(); } + /* test action */ + for (j = 0; j < 4; j++) { - status |= gr_randtest(f, state, cctx); + status |= gr_randtest_not_zero(f, state, cctx); status |= gr_ore_poly_apply(g_e, P_e, f, ctx_e); status |= gr_ore_poly_apply(g_d, P_d, f, ctx_d); - /* euler_to_ddx leaves the operator unchanged: P_d(f) = P_e(f) */ if (status == GR_SUCCESS && gr_equal(g_d, g_e, cctx) == T_FALSE) { - flint_printf("FAIL: euler_to_ddx application identity\n"); + flint_printf("FAIL: application identity\n"); flint_abort(); } } - /* these representative base rings never fail for the implemented var 0 */ - int expect_success = (var == 0 && cctx->which_ring == GR_CTX_GR_POLY - && (POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_FMPZ - || POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_CC_ACB - || POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_NMOD)); - if (expect_success && status != GR_SUCCESS) - { - flint_printf("FAIL: euler_to_ddx unexpected failure\n"); - flint_abort(); - } - - count_success += (status == GR_SUCCESS); - count_domain += ((status & GR_DOMAIN) != 0); - count_unable += ((status & GR_UNABLE) != 0); - +cleanup: gr_heap_clear(xpow, cctx); gr_heap_clear(f, cctx); gr_heap_clear(g_e, cctx); @@ -135,6 +138,7 @@ TEST_GR_FUNCTION_START(gr_ore_poly_euler_to_ddx, state, count_success, count_dom gr_ore_poly_clear(scaled_P, ctx_e); gr_ore_poly_ctx_clear(ctx_d); gr_ore_poly_ctx_clear(ctx_e); + gr_vec_clear(gens, cctx); gr_ctx_clear(cctx); } diff --git a/src/gr_ore_poly/test/t-shift_convert.c b/src/gr_ore_poly/test/t-shift_convert.c index b30b7979dc..a43b049963 100644 --- a/src/gr_ore_poly/test/t-shift_convert.c +++ b/src/gr_ore_poly/test/t-shift_convert.c @@ -19,6 +19,12 @@ static const ore_algebra_t shift_algs[4] = { ORE_ALGEBRA_FORWARD_DIFFERENCE, ORE_ALGEBRA_BACKWARD_DIFFERENCE }; +static int +is_difference(ore_algebra_t a) +{ + return a == ORE_ALGEBRA_FORWARD_DIFFERENCE || a == ORE_ALGEBRA_BACKWARD_DIFFERENCE; +} + TEST_GR_FUNCTION_START(gr_ore_poly_shift_convert, state, count_success, count_domain, count_unable) { for (slong iter = 0; iter < 1000 * flint_test_multiplier(); iter++) @@ -26,19 +32,15 @@ TEST_GR_FUNCTION_START(gr_ore_poly_shift_convert, state, count_success, count_do gr_ctx_t cctx; gr_ore_poly_ctx_t ctx_src, ctx_dst; ore_algebra_t src_alg, dst_alg; - slong len, j, epow, epow2, ngens, var; + slong len, j, sh, sh2, ngens, var; int status = GR_SUCCESS; - /* occasionally use a multivariate base ring so that the conversion is - exercised with a generator other than the default one (which is not - yet implemented and returns GR_UNABLE) */ switch (n_randint(state, 8)) { case 0: gr_ctx_init_random(cctx, state); break; case 1: - case 2: gr_ctx_init_random_mpoly(cctx, state); break; default: @@ -48,7 +50,7 @@ TEST_GR_FUNCTION_START(gr_ore_poly_shift_convert, state, count_success, count_do ngens = 0; status |= gr_ctx_ngens(&ngens, cctx); - var = (ngens >= 1) ? (slong) n_randint(state, ngens) : 0; + var = n_randint(state, ngens); src_alg = shift_algs[n_randint(state, 4)]; dst_alg = shift_algs[n_randint(state, 4)]; @@ -65,83 +67,78 @@ TEST_GR_FUNCTION_START(gr_ore_poly_shift_convert, state, count_success, count_do len = P->length; gr_ore_poly_fit_length(Pd, len, ctx_dst); - status |= _gr_ore_poly_shift_convert(Pd->coeffs, &epow, P->coeffs, len, src_alg, dst_alg, var, cctx); + if (is_difference(src_alg) && is_difference(dst_alg) + && src_alg != dst_alg && n_randint(state, 2)) + status |= _gr_ore_poly_shift_convert_difference(Pd->coeffs, &sh, + P->coeffs, len, dst_alg == ORE_ALGEBRA_BACKWARD_DIFFERENCE, + var, cctx); + else + status |= _gr_ore_poly_shift_convert(Pd->coeffs, &sh, P->coeffs, + len, src_alg, dst_alg, var, + cctx); _gr_ore_poly_set_length(Pd, len, ctx_dst); _gr_ore_poly_normalise(Pd, ctx_dst); gr_ore_poly_fit_length(P2, len, ctx_src); - status |= _gr_ore_poly_shift_convert(P2->coeffs, &epow2, Pd->coeffs, len, dst_alg, src_alg, var, cctx); + status |= _gr_ore_poly_shift_convert(P2->coeffs, &sh2, Pd->coeffs, len, + dst_alg, src_alg, var, cctx); _gr_ore_poly_set_length(P2, len, ctx_src); _gr_ore_poly_normalise(P2, ctx_src); - /* these representative base rings never fail for the implemented var 0 */ + count_success += (status == GR_SUCCESS); + count_domain += ((status & GR_DOMAIN) != 0); + count_unable += ((status & GR_UNABLE) != 0); + int expect_success = (var == 0 && cctx->which_ring == GR_CTX_GR_POLY && (POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_FMPZ || POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_CC_ACB || POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_NMOD)); if (expect_success && status != GR_SUCCESS) { - flint_printf("FAIL: shift_convert unexpected failure\n"); + flint_printf("FAIL: unexpected failure\n"); flint_abort(); } - /* the round trip recovers the operator: the powers of S cancel */ - if (status == GR_SUCCESS && epow + epow2 != 0) + if (status == GR_SUCCESS && sh + sh2 != 0) { - flint_printf("FAIL: shift_convert round-trip power\n"); + flint_printf("FAIL: round-trip power\n"); flint_abort(); } if (status == GR_SUCCESS && gr_ore_poly_equal(P2, P, ctx_src) == T_FALSE) { - flint_printf("FAIL: shift_convert round trip\n"); + flint_printf("FAIL: round trip\n"); flint_abort(); } - /* application identity dst(f) = S^{-epow}(src(f)), over polynomial base - rings (S^epow * dst = src as operators, so the correction is S^{-epow} - applied on the left of the value) */ if (cctx->which_ring == GR_CTX_GR_POLY) { gr_ctx_struct * sctx = POLYNOMIAL_ELEM_CTX(cctx); gr_ptr f = gr_heap_init(cctx); - gr_ptr shifted = gr_heap_init(cctx); gr_ptr g1 = gr_heap_init(cctx); gr_ptr g2 = gr_heap_init(cctx); gr_ptr c = gr_heap_init(sctx); - status |= gr_set_si(c, -epow, sctx); + status |= gr_set_si(c, -sh, sctx); for (j = 0; j < 4; j++) { - status |= gr_randtest(f, state, cctx); + status |= gr_randtest_not_zero(f, state, cctx); status |= gr_ore_poly_apply(g1, Pd, f, ctx_dst); status |= gr_ore_poly_apply(g2, P, f, ctx_src); - status |= gr_poly_taylor_shift((gr_poly_struct *) shifted, (gr_poly_struct *) g2, c, sctx); - if (status == GR_SUCCESS && gr_equal(g1, shifted, cctx) == T_FALSE) + status |= gr_poly_taylor_shift((gr_poly_struct *) g2, (gr_poly_struct *) g2, c, sctx); + if (status == GR_SUCCESS && gr_equal(g1, g2, cctx) == T_FALSE) { - flint_printf("FAIL: shift_convert application identity\n"); + flint_printf("FAIL: application identity\n"); flint_abort(); } } gr_heap_clear(f, cctx); - gr_heap_clear(shifted, cctx); gr_heap_clear(g1, cctx); gr_heap_clear(g2, cctx); gr_heap_clear(c, sctx); } - /* a non-shift algebra type is rejected cleanly */ - if (_gr_ore_poly_shift_convert(P2->coeffs, &epow2, P->coeffs, len, ORE_ALGEBRA_DERIVATIVE, dst_alg, var, cctx) != GR_DOMAIN) - { - flint_printf("FAIL: shift_convert expected GR_DOMAIN\n"); - flint_abort(); - } - - count_success += (status == GR_SUCCESS); - count_domain += ((status & GR_DOMAIN) != 0); - count_unable += ((status & GR_UNABLE) != 0); - gr_ore_poly_clear(P, ctx_src); gr_ore_poly_clear(Pd, ctx_dst); gr_ore_poly_clear(P2, ctx_src); From 87152192503f671fd6c268bdadca2e485d14c978 Mon Sep 17 00:00:00 2001 From: Marc Mezzarobba Date: Mon, 6 Jul 2026 10:10:02 +0200 Subject: [PATCH 4/7] gr_ore_poly: test shift_to_differential This commit contains both AI-generated code and human changes. Co-Authored-By: Claude Opus 4.8 --- src/gr_ore_poly/test/main.c | 3 +- src/gr_ore_poly/test/t-differential_shift.c | 382 ++++++++++++++++++ .../test/t-differential_to_shift.c | 196 --------- 3 files changed, 384 insertions(+), 197 deletions(-) create mode 100644 src/gr_ore_poly/test/t-differential_shift.c delete mode 100644 src/gr_ore_poly/test/t-differential_to_shift.c diff --git a/src/gr_ore_poly/test/main.c b/src/gr_ore_poly/test/main.c index b13ee8cc9b..f1d2f39563 100644 --- a/src/gr_ore_poly/test/main.c +++ b/src/gr_ore_poly/test/main.c @@ -20,7 +20,7 @@ #include "t-ddx_to_euler.c" #include "t-euler_to_ddx.c" #include "t-shift_convert.c" -#include "t-differential_to_shift.c" +#include "t-differential_shift.c" #include "t-convert.c" /* Array of test functions ***************************************************/ @@ -37,6 +37,7 @@ test_struct tests[] = TEST_FUNCTION(gr_ore_poly_euler_to_ddx), TEST_FUNCTION(gr_ore_poly_shift_convert), TEST_FUNCTION(gr_ore_poly_differential_to_shift), + TEST_FUNCTION(gr_ore_poly_shift_to_differential), TEST_FUNCTION(gr_ore_poly_convert) }; diff --git a/src/gr_ore_poly/test/t-differential_shift.c b/src/gr_ore_poly/test/t-differential_shift.c new file mode 100644 index 0000000000..2d4c398b0e --- /dev/null +++ b/src/gr_ore_poly/test/t-differential_shift.c @@ -0,0 +1,382 @@ +/* + 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 . +*/ + +/* generated using Claude Opus 4.8 */ + +#include "test_helpers.h" +#include "gr.h" +#include "gr_poly.h" +#include "gr_ore_poly.h" + +static const ore_algebra_t ds_diff_algs[2] = { + ORE_ALGEBRA_DERIVATIVE, ORE_ALGEBRA_EULER_DERIVATIVE +}; + +static const ore_algebra_t ds_shift_algs[4] = { + ORE_ALGEBRA_FORWARD_SHIFT, ORE_ALGEBRA_BACKWARD_SHIFT, + ORE_ALGEBRA_FORWARD_DIFFERENCE, ORE_ALGEBRA_BACKWARD_DIFFERENCE +}; + +/* Returns 0 only if a == x^net * b is provably false */ +static int +check_scaled(const gr_ore_poly_t a, const gr_ore_poly_t b, slong net, + gr_ore_poly_ctx_t octx) +{ + int ok = 1, status = GR_SUCCESS; + gr_ctx_struct * cctx = GR_ORE_POLY_ELEM_CTX(octx); + slong sz = cctx->sizeof_elem, i; + const gr_ore_poly_struct * lo = (net >= 0) ? b : a; + const gr_ore_poly_struct * hi = (net >= 0) ? a : b; + gr_ptr xp; + gr_ore_poly_t scaled; + + GR_TMP_INIT(xp, cctx); + gr_ore_poly_init(scaled, octx); + + status |= gr_gen(xp, cctx); + status |= gr_pow_ui(xp, xp, (ulong) FLINT_ABS(net), cctx); + gr_ore_poly_fit_length(scaled, lo->length, octx); + for (i = 0; i < lo->length; i++) + status |= gr_mul(GR_ENTRY(scaled->coeffs, i, sz), xp, GR_ENTRY(lo->coeffs, i, sz), cctx); + _gr_ore_poly_set_length(scaled, lo->length, octx); + _gr_ore_poly_normalise(scaled, octx); + + if (status == GR_SUCCESS && gr_ore_poly_equal(scaled, hi, octx) == T_FALSE) + ok = 0; + + gr_ore_poly_clear(scaled, octx); + GR_TMP_CLEAR(xp, cctx); + return ok; +} + +/* Returns 0 only if a == S^net * b is provably false, where S is the forward + shift operator. Since S is not an element of octx when it is a backward + algebra, we compare in an auxiliary algebra in which S is the generator: we + convert both operands to it and clear the S^p offsets returned by convert. */ +static int +check_shifted(const gr_ore_poly_t a, const gr_ore_poly_t b, slong net, + gr_ore_poly_ctx_t octx) +{ + int ok = 1, status = GR_SUCCESS; + slong var = GR_ORE_POLY_ORE_DATA(octx)->base_var; + gr_ctx_struct * base = GR_ORE_POLY_ELEM_CTX(octx); + gr_ore_poly_ctx_t fs_ctx; + gr_ore_poly_t af, bf, s, lhs, rhs; + slong ea, eb, la, lb, c; + + gr_ore_poly_ctx_init(fs_ctx, base, var, ORE_ALGEBRA_FORWARD_SHIFT); + gr_ore_poly_init(af, fs_ctx); + gr_ore_poly_init(bf, fs_ctx); + gr_ore_poly_init(s, fs_ctx); + gr_ore_poly_init(lhs, fs_ctx); + gr_ore_poly_init(rhs, fs_ctx); + + /* S^ea * af = a and S^eb * bf = b */ + status |= gr_ore_poly_convert(af, &ea, a, fs_ctx, octx); + status |= gr_ore_poly_convert(bf, &eb, b, fs_ctx, octx); + + /* a == S^net * b becomes S^(ea+c) af == S^(net+eb+c) bf, with c shifting + both exponents to be nonnegative so the powers of S are representable */ + la = ea; + lb = net + eb; + c = FLINT_MAX(0, -FLINT_MIN(la, lb)); + + status |= gr_ore_poly_gen(s, fs_ctx); + status |= gr_pow_ui(lhs, s, (ulong) (la + c), fs_ctx); + status |= gr_ore_poly_mul(lhs, lhs, af, fs_ctx); + status |= gr_pow_ui(rhs, s, (ulong) (lb + c), fs_ctx); + status |= gr_ore_poly_mul(rhs, rhs, bf, fs_ctx); + + if (status == GR_SUCCESS && gr_ore_poly_equal(lhs, rhs, fs_ctx) == T_FALSE) + ok = 0; + + gr_ore_poly_clear(af, fs_ctx); + gr_ore_poly_clear(bf, fs_ctx); + gr_ore_poly_clear(s, fs_ctx); + gr_ore_poly_clear(lhs, fs_ctx); + gr_ore_poly_clear(rhs, fs_ctx); + gr_ore_poly_ctx_clear(fs_ctx); + return ok; +} + +/* Applies a forward-shift recurrence R = sum_k q_k(nu) S^k to the coefficients + of a polynomial. */ +static int +apply_forward_recurrence(gr_poly_t s, const gr_ore_poly_t R, const gr_poly_t f, gr_ctx_t cctx) +{ + gr_ctx_struct * sctx = POLYNOMIAL_ELEM_CTX(cctx); + slong bsz = cctx->sizeof_elem, ssz = sctx->sizeof_elem; + slong df = f->length; + int status = GR_SUCCESS; + gr_ptr mval, qval, term; + + GR_TMP_INIT3(mval, qval, term, sctx); + gr_poly_fit_length(s, df, sctx); + for (slong n = 0; n < df; n++) + { + gr_ptr sm = GR_ENTRY(s->coeffs, n, ssz); + status |= gr_zero(sm, sctx); + status |= gr_set_si(mval, n, sctx); + for (slong k = 0; k < R->length && n + k < df; k++) + { + const gr_poly_struct * qk = (const gr_poly_struct *) GR_ENTRY(R->coeffs, k, bsz); + status |= gr_poly_evaluate(qval, qk, mval, sctx); + status |= gr_mul(term, qval, GR_ENTRY(f->coeffs, n + k, ssz), sctx); + status |= gr_add(sm, sm, term, sctx); + } + } + _gr_poly_set_length(s, df, sctx); + _gr_poly_normalise(s, sctx); + GR_TMP_CLEAR3(mval, qval, term, sctx); + return status; +} + +TEST_GR_FUNCTION_START(gr_ore_poly_differential_to_shift, state, count_success, count_domain, count_unable) +{ + for (slong iter = 0; iter < 1000 * flint_test_multiplier(); iter++) + { + gr_ctx_t cctx; + gr_ore_poly_ctx_t ctx_d, ctx_s; + ore_algebra_t diff_alg, shift_alg; + slong p1, p2, ngens, var; + int status = GR_SUCCESS; + int expect_success = 1; + + switch (n_randint(state, 8)) + { + case 0: + gr_ctx_init_random(cctx, state); + break; + case 1: + gr_ctx_init_random_mpoly(cctx, state); + break; + default: + gr_ctx_init_random_poly(cctx, state); + break; + } + + ngens = 0; + status |= gr_ctx_ngens(&ngens, cctx); + var = n_randint(state, ngens); + + diff_alg = ds_diff_algs[n_randint(state, 2)]; + shift_alg = ds_shift_algs[n_randint(state, 4)]; + + if (n_randint(state, 16)) + { + diff_alg = ORE_ALGEBRA_COMMUTATIVE; + expect_success = 0; + } + if (n_randint(state, 16)) + { + shift_alg = ORE_ALGEBRA_COMMUTATIVE; + expect_success = 0; + } + + gr_ore_poly_ctx_init(ctx_d, cctx, var, diff_alg); + gr_ore_poly_ctx_init(ctx_s, cctx, var, shift_alg); + + gr_ore_poly_t op, res_s, op2; + gr_ore_poly_init(op, ctx_d); + gr_ore_poly_init(res_s, ctx_s); + gr_ore_poly_init(op2, ctx_d); + + status |= gr_ore_poly_randtest(op, state, 1 + n_randint(state, 5), ctx_d); + + status |= gr_ore_poly_differential_to_shift(res_s, &p1, op, ctx_s, ctx_d); + status |= gr_ore_poly_shift_to_differential(op2, &p2, res_s, ctx_d, ctx_s); + + expect_success = (expect_success && var == 0 + && cctx->which_ring == GR_CTX_GR_POLY + && (POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_FMPZ + || POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_CC_ACB + || POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_NMOD)); + if (expect_success && status != GR_SUCCESS) + { + flint_printf("FAIL: unexpected failure\n"); + flint_abort(); + } + + count_success += (status == GR_SUCCESS); + count_domain += ((status & GR_DOMAIN) != 0); + count_unable += ((status & GR_UNABLE) != 0); + + /* round trip recovers x^(p1-p2) * op */ + + if (status == GR_SUCCESS && !check_scaled(op2, op, p1 - p2, ctx_d)) + { + flint_printf("FAIL: round trip\n"); + flint_abort(); + } + + /* a forward-shift recurrence R from op satisfies + (R a)_m = coeff_{m-pf}(op . f) for the coefficient sequence of f */ + + if (cctx->which_ring == GR_CTX_GR_POLY && shift_alg == ORE_ALGEBRA_FORWARD_SHIFT) + { + gr_ctx_struct * sctx = POLYNOMIAL_ELEM_CTX(cctx); + gr_ptr f = gr_heap_init(cctx); + gr_ptr g = gr_heap_init(cctx); + gr_poly_t s, eg; + + gr_poly_init(s, sctx); + gr_poly_init(eg, sctx); + + status |= gr_randtest(f, state, cctx); + status |= gr_ore_poly_apply(g, op, f, ctx_d); + status |= apply_forward_recurrence(s, res_s, (gr_poly_struct *) f, cctx); + if (p1 >= 0) + status |= gr_poly_shift_left(eg, (gr_poly_struct *) g, p1, sctx); + else + status |= gr_poly_shift_right(eg, (gr_poly_struct *) g, -p1, sctx); + + if (status == GR_SUCCESS && gr_poly_equal(s, eg, sctx) == T_FALSE) + { + flint_printf("FAIL: action\n"); + flint_abort(); + } + + gr_poly_clear(s, sctx); + gr_poly_clear(eg, sctx); + gr_heap_clear(f, cctx); + gr_heap_clear(g, cctx); + } + + gr_ore_poly_clear(op, ctx_d); + gr_ore_poly_clear(res_s, ctx_s); + gr_ore_poly_clear(op2, ctx_d); + gr_ore_poly_ctx_clear(ctx_d); + gr_ore_poly_ctx_clear(ctx_s); + gr_ctx_clear(cctx); + } + + TEST_GR_FUNCTION_END(state, count_success, count_domain, count_unable); +} + +TEST_GR_FUNCTION_START(gr_ore_poly_shift_to_differential, state, count_success, count_domain, count_unable) +{ + for (slong iter = 0; iter < 1000 * flint_test_multiplier(); iter++) + { + gr_ctx_t cctx; + gr_ore_poly_ctx_t ctx_d, ctx_s; + ore_algebra_t diff_alg, shift_alg; + slong p1, p2, ngens, var; + int status = GR_SUCCESS; + int expect_success = 1; + + switch (n_randint(state, 8)) + { + case 0: + gr_ctx_init_random(cctx, state); + break; + case 1: + gr_ctx_init_random_mpoly(cctx, state); + break; + default: + gr_ctx_init_random_poly(cctx, state); + break; + } + + ngens = 0; + status |= gr_ctx_ngens(&ngens, cctx); + var = n_randint(state, ngens); + + diff_alg = ds_diff_algs[n_randint(state, 2)]; + shift_alg = ds_shift_algs[n_randint(state, 4)]; + + + if (n_randint(state, 16)) + { + diff_alg = ORE_ALGEBRA_COMMUTATIVE; + expect_success = 0; + } + if (n_randint(state, 16)) + { + shift_alg = ORE_ALGEBRA_COMMUTATIVE; + expect_success = 0; + } + + gr_ore_poly_ctx_init(ctx_d, cctx, var, diff_alg); + gr_ore_poly_ctx_init(ctx_s, cctx, var, shift_alg); + + gr_ore_poly_t op, res_d, op2; + gr_ore_poly_init(op, ctx_s); + gr_ore_poly_init(res_d, ctx_d); + gr_ore_poly_init(op2, ctx_s); + + status |= gr_ore_poly_randtest(op, state, 1 + n_randint(state, 5), ctx_s); + + status |= gr_ore_poly_shift_to_differential(res_d, &p1, op, ctx_d, ctx_s); + status |= gr_ore_poly_differential_to_shift(op2, &p2, res_d, ctx_s, ctx_d); + + expect_success = (expect_success && var == 0 + && cctx->which_ring == GR_CTX_GR_POLY + && (POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_FMPZ + || POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_CC_ACB + || POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_NMOD)); + if (expect_success && status != GR_SUCCESS) + { + flint_printf("FAIL: unexpected failure\n"); + flint_abort(); + } + + count_success += (status == GR_SUCCESS); + count_domain += ((status & GR_DOMAIN) != 0); + count_unable += ((status & GR_UNABLE) != 0); + + /* round trip recovers S^(p1-p2) * op */ + + if (status == GR_SUCCESS && !check_shifted(op2, op, p1 - p2, ctx_s)) + { + flint_printf("FAIL: round trip\n"); + flint_abort(); + } + + if (cctx->which_ring == GR_CTX_GR_POLY && shift_alg == ORE_ALGEBRA_FORWARD_SHIFT) + { + gr_ctx_struct * sctx = POLYNOMIAL_ELEM_CTX(cctx); + gr_ptr g = gr_heap_init(cctx); + gr_poly_t s, eg; + + gr_poly_init(s, sctx); + gr_poly_init(eg, sctx); + + slong flen = op->length + FLINT_ABS(p1) + 2 + n_randint(state, 5); + gr_ptr f = gr_heap_init(cctx); + status |= gr_poly_randtest((gr_poly_struct *) f, state, flen, sctx); + + status |= gr_ore_poly_apply(g, res_d, f, ctx_d); + status |= apply_forward_recurrence(s, op, (gr_poly_struct *) f, cctx); + if (p1 >= 0) + status |= gr_poly_shift_left(eg, (gr_poly_struct *) g, p1, sctx); + else + status |= gr_poly_shift_right(eg, (gr_poly_struct *) g, -p1, sctx); + + if (status == GR_SUCCESS && gr_poly_equal(s, eg, sctx) == T_FALSE) + { + flint_printf("FAIL: action\n"); + flint_abort(); + } + + gr_poly_clear(s, sctx); + gr_poly_clear(eg, sctx); + gr_heap_clear(f, cctx); + gr_heap_clear(g, cctx); + } + + gr_ore_poly_clear(op, ctx_s); + gr_ore_poly_clear(res_d, ctx_d); + gr_ore_poly_clear(op2, ctx_s); + gr_ore_poly_ctx_clear(ctx_d); + gr_ore_poly_ctx_clear(ctx_s); + gr_ctx_clear(cctx); + } + + TEST_GR_FUNCTION_END(state, count_success, count_domain, count_unable); +} diff --git a/src/gr_ore_poly/test/t-differential_to_shift.c b/src/gr_ore_poly/test/t-differential_to_shift.c deleted file mode 100644 index 75171ab22d..0000000000 --- a/src/gr_ore_poly/test/t-differential_to_shift.c +++ /dev/null @@ -1,196 +0,0 @@ -/* - 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 . -*/ - -/* generated using Claude Opus 4.8 */ - -#include "test_helpers.h" -#include "gr.h" -#include "gr_poly.h" -#include "gr_ore_poly.h" - -static const ore_algebra_t ds_diff_algs[2] = { - ORE_ALGEBRA_DERIVATIVE, ORE_ALGEBRA_EULER_DERIVATIVE -}; - -static const ore_algebra_t ds_shift_algs[4] = { - ORE_ALGEBRA_FORWARD_SHIFT, ORE_ALGEBRA_BACKWARD_SHIFT, - ORE_ALGEBRA_FORWARD_DIFFERENCE, ORE_ALGEBRA_BACKWARD_DIFFERENCE -}; - -/* Returns 0 only if a == x^net * b is provably false (x = gen of cctx). */ -static int -check_scaled(const gr_ore_poly_t a, const gr_ore_poly_t b, slong net, - gr_ore_poly_ctx_t octx, gr_ctx_t cctx) -{ - int ok = 1, status = GR_SUCCESS; - slong sz = cctx->sizeof_elem, i; - const gr_ore_poly_struct * lo = (net >= 0) ? b : a; - const gr_ore_poly_struct * hi = (net >= 0) ? a : b; - gr_ptr xp; - gr_ore_poly_t scaled; - - GR_TMP_INIT(xp, cctx); - gr_ore_poly_init(scaled, octx); - - status |= gr_gen(xp, cctx); - status |= gr_pow_ui(xp, xp, (ulong) FLINT_ABS(net), cctx); - gr_ore_poly_fit_length(scaled, lo->length, octx); - for (i = 0; i < lo->length; i++) - status |= gr_mul(GR_ENTRY(scaled->coeffs, i, sz), xp, GR_ENTRY(lo->coeffs, i, sz), cctx); - _gr_ore_poly_set_length(scaled, lo->length, octx); - _gr_ore_poly_normalise(scaled, octx); - - if (status == GR_SUCCESS && gr_ore_poly_equal(scaled, hi, octx) == T_FALSE) - ok = 0; - - gr_ore_poly_clear(scaled, octx); - GR_TMP_CLEAR(xp, cctx); - return ok; -} - -/* Applies a forward-shift recurrence R = sum_k q_k(nu) S^k to the coefficients - of a polynomial. */ -static int -apply_forward_recurrence(gr_poly_t s, const gr_ore_poly_t R, const gr_poly_t f, gr_ctx_t cctx) -{ - gr_ctx_struct * sctx = POLYNOMIAL_ELEM_CTX(cctx); - slong bsz = cctx->sizeof_elem, ssz = sctx->sizeof_elem; - slong df = f->length; - int status = GR_SUCCESS; - gr_ptr mval, qval, term; - - GR_TMP_INIT3(mval, qval, term, sctx); - gr_poly_fit_length(s, df, sctx); - for (slong n = 0; n < df; n++) - { - gr_ptr sm = GR_ENTRY(s->coeffs, n, ssz); - status |= gr_zero(sm, sctx); - status |= gr_set_si(mval, n, sctx); - for (slong k = 0; k < R->length && n + k < df; k++) - { - const gr_poly_struct * qk = (const gr_poly_struct *) GR_ENTRY(R->coeffs, k, bsz); - status |= gr_poly_evaluate(qval, qk, mval, sctx); - status |= gr_mul(term, qval, GR_ENTRY(f->coeffs, n + k, ssz), sctx); - status |= gr_add(sm, sm, term, sctx); - } - } - _gr_poly_set_length(s, df, sctx); - _gr_poly_normalise(s, sctx); - GR_TMP_CLEAR3(mval, qval, term, sctx); - return status; -} - -TEST_GR_FUNCTION_START(gr_ore_poly_differential_to_shift, state, count_success, count_domain, count_unable) -{ - for (slong iter = 0; iter < 1000 * flint_test_multiplier(); iter++) - { - gr_ctx_t cctx; - gr_ore_poly_ctx_t ctx_d, ctx_s; - ore_algebra_t diff_alg, shift_alg; - slong p1, p2, ngens, var; - int status = GR_SUCCESS; - - switch (n_randint(state, 8)) - { - case 0: - gr_ctx_init_random(cctx, state); - break; - case 1: - gr_ctx_init_random_mpoly(cctx, state); - break; - default: - gr_ctx_init_random_poly(cctx, state); - break; - } - - ngens = 0; - status |= gr_ctx_ngens(&ngens, cctx); - var = n_randint(state, ngens); - - diff_alg = ds_diff_algs[n_randint(state, 2)]; - shift_alg = ds_shift_algs[n_randint(state, 4)]; - - gr_ore_poly_ctx_init(ctx_d, cctx, var, diff_alg); - gr_ore_poly_ctx_init(ctx_s, cctx, var, shift_alg); - - gr_ore_poly_t op, res_s, op2; - gr_ore_poly_init(op, ctx_d); - gr_ore_poly_init(res_s, ctx_s); - gr_ore_poly_init(op2, ctx_d); - - status |= gr_ore_poly_randtest(op, state, 1 + n_randint(state, 5), ctx_d); - - status |= gr_ore_poly_differential_to_shift(res_s, &p1, op, ctx_s, ctx_d); - status |= gr_ore_poly_shift_to_differential(op2, &p2, res_s, ctx_d, ctx_s); - - int expect_success = (var == 0 && cctx->which_ring == GR_CTX_GR_POLY - && (POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_FMPZ - || POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_CC_ACB - || POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_NMOD)); - if (expect_success && status != GR_SUCCESS) - { - flint_printf("FAIL: differential<->shift unexpected failure\n"); - flint_abort(); - } - - count_success += (status == GR_SUCCESS); - count_domain += ((status & GR_DOMAIN) != 0); - count_unable += ((status & GR_UNABLE) != 0); - - /* round trip recovers x^(p1-p2) * op */ - - if (status == GR_SUCCESS && !check_scaled(op2, op, p1 - p2, ctx_d, cctx)) - { - flint_printf("FAIL: round trip\n"); - flint_abort(); - } - - /* a forward-shift recurrence R from op satisfies - (R a)_m = coeff_{m-pf}(op . f) for the coefficient sequence of f */ - - if (cctx->which_ring == GR_CTX_GR_POLY && shift_alg == ORE_ALGEBRA_FORWARD_SHIFT) - { - gr_ctx_struct * sctx = POLYNOMIAL_ELEM_CTX(cctx); - gr_ptr f = gr_heap_init(cctx); - gr_ptr g = gr_heap_init(cctx); - gr_poly_t s, eg; - - gr_poly_init(s, sctx); - gr_poly_init(eg, sctx); - - status |= gr_randtest(f, state, cctx); - status |= gr_ore_poly_apply(g, op, f, ctx_d); - status |= apply_forward_recurrence(s, res_s, (gr_poly_struct *) f, cctx); - if (p1 >= 0) - status |= gr_poly_shift_left(eg, (gr_poly_struct *) g, p1, sctx); - else - status |= gr_poly_shift_right(eg, (gr_poly_struct *) g, -p1, sctx); - - if (status == GR_SUCCESS && gr_poly_equal(s, eg, sctx) == T_FALSE) - { - flint_printf("FAIL: action\n"); - flint_abort(); - } - - gr_poly_clear(s, sctx); - gr_poly_clear(eg, sctx); - gr_heap_clear(f, cctx); - gr_heap_clear(g, cctx); - } - - gr_ore_poly_clear(op, ctx_d); - gr_ore_poly_clear(res_s, ctx_s); - gr_ore_poly_clear(op2, ctx_d); - gr_ore_poly_ctx_clear(ctx_d); - gr_ore_poly_ctx_clear(ctx_s); - gr_ctx_clear(cctx); - } - - TEST_GR_FUNCTION_END(state, count_success, count_domain, count_unable); -} From da8b43303f521b927dc28095231cea5bba9872e3 Mon Sep 17 00:00:00 2001 From: "Marc Mezzarobba (AI)" Date: Sun, 5 Jul 2026 18:57:30 +0200 Subject: [PATCH 5/7] gr_ore_poly: fix convert The just-added implementation was incorrect for Ore polynomial rings of the same algebra type but with unequal parameters. Co-Authored-By: Claude Fable 5 --- src/gr_ore_poly/convert.c | 20 +++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/src/gr_ore_poly/convert.c b/src/gr_ore_poly/convert.c index df7a1e17fa..8c45f98220 100644 --- a/src/gr_ore_poly/convert.c +++ b/src/gr_ore_poly/convert.c @@ -12,6 +12,12 @@ #include "gr_vec.h" #include "gr_ore_poly.h" +static int +is_diff_case(ore_algebra_t a) +{ + return a == ORE_ALGEBRA_DERIVATIVE || a == ORE_ALGEBRA_EULER_DERIVATIVE; +} + static int is_shift_case(ore_algebra_t a) { @@ -39,11 +45,23 @@ gr_ore_poly_convert(gr_ore_poly_t res, slong * p, const gr_ore_poly_t op, return gr_ore_poly_zero(res, res_ctx); gr_ore_poly_ore_data_t * od = GR_ORE_POLY_ORE_DATA(op_ctx); + gr_ore_poly_ore_data_t * rd = GR_ORE_POLY_ORE_DATA(res_ctx); slong var = (od != NULL) ? od->base_var : -1; gr_ore_poly_fit_length(res, len, res_ctx); - if (sa == da) + if (res_ctx == op_ctx) + status |= _gr_vec_set(res->coeffs, op->coeffs, len, base); + /* Distinct contexts with equal which_algebra need not define the same + algebra (e.g., two q-shift algebras with different q, or operators + acting on different variables of the base ring). Only the differential + and shift/difference algebras are fully determined by which_algebra and + base_var. */ + else if (!(is_diff_case(sa) || is_shift_case(sa)) + || !(is_diff_case(da) || is_shift_case(da)) + || od->base_var != rd->base_var) + return GR_UNABLE; + else if (sa == da) status |= _gr_vec_set(res->coeffs, op->coeffs, len, base); else if (sa == ORE_ALGEBRA_DERIVATIVE && da == ORE_ALGEBRA_EULER_DERIVATIVE) { From eeda48a56a049c7b295600e60e56ff8d32dd731c Mon Sep 17 00:00:00 2001 From: Marc Mezzarobba Date: Sun, 5 Jul 2026 18:10:05 +0200 Subject: [PATCH 6/7] gr_ore_poly: revise AI-generated code --- doc/source/gr_ore_poly.rst | 4 +- src/gr_ore_poly.h | 2 + src/gr_ore_poly/convert.c | 20 ++++--- src/gr_ore_poly/ctx.c | 3 - src/gr_ore_poly/test/t-convert.c | 63 +++++++-------------- src/gr_ore_poly/test/t-differential_shift.c | 21 +++---- src/gr_ore_poly/test/t-shift_convert.c | 1 + 7 files changed, 44 insertions(+), 70 deletions(-) diff --git a/doc/source/gr_ore_poly.rst b/doc/source/gr_ore_poly.rst index a422449f62..5bf5380b0a 100644 --- a/doc/source/gr_ore_poly.rst +++ b/doc/source/gr_ore_poly.rst @@ -307,9 +307,9 @@ different Ore polynomial rings over the same base ring. to operators written in terms of the forward and backward shifts `S`, `S^{-1}` and the forward and backward differences `S-1`, `1-S^{-1}`). The context *ctx* is the common base ring and *var* is the index of the - generator of *ctx* on which `S` acts. Conversions that crosses between the + generator of *ctx* on which `S` acts. Conversions that cross between the forward side `S`, `S-1` and the backward side `S^{-1}`, `1-S^{-1}` - additionally require a generic univariate polynomial base ring and otherwise + currently require a generic univariate polynomial base ring and otherwise return ``GR_UNABLE``. The result satisfies `S^{\textit{p}} \cdot \textit{res} = \textit{op}`. diff --git a/src/gr_ore_poly.h b/src/gr_ore_poly.h index 3266fca112..3a3799c431 100644 --- a/src/gr_ore_poly.h +++ b/src/gr_ore_poly.h @@ -250,6 +250,8 @@ extern const gr_ore_poly_sigma_delta_t _gr_ore_poly_default_sigma_delta[]; WARN_UNUSED_RESULT int gr_ore_poly_apply_custom(gr_ptr res, const gr_ore_poly_t P, gr_srcptr f, gr_srcptr d1, gr_ore_poly_ctx_t ctx); WARN_UNUSED_RESULT int gr_ore_poly_apply(gr_ptr res, const gr_ore_poly_t P, gr_srcptr f, gr_ore_poly_ctx_t ctx); +/* Conversions */ + WARN_UNUSED_RESULT int _gr_ore_poly_ddx_to_euler(gr_ptr res, gr_srcptr op, slong len, slong var, gr_ctx_t ctx); WARN_UNUSED_RESULT int _gr_ore_poly_euler_to_ddx(gr_ptr res, gr_srcptr op, slong len, slong var, gr_ctx_t ctx); WARN_UNUSED_RESULT int _gr_ore_poly_shift_convert(gr_ptr res, slong * p, gr_srcptr op, slong len, ore_algebra_t src_alg, ore_algebra_t dst_alg, slong var, gr_ctx_t ctx); diff --git a/src/gr_ore_poly/convert.c b/src/gr_ore_poly/convert.c index 8c45f98220..f66d438a42 100644 --- a/src/gr_ore_poly/convert.c +++ b/src/gr_ore_poly/convert.c @@ -51,17 +51,19 @@ gr_ore_poly_convert(gr_ore_poly_t res, slong * p, const gr_ore_poly_t op, gr_ore_poly_fit_length(res, len, res_ctx); if (res_ctx == op_ctx) + { status |= _gr_vec_set(res->coeffs, op->coeffs, len, base); - /* Distinct contexts with equal which_algebra need not define the same - algebra (e.g., two q-shift algebras with different q, or operators - acting on different variables of the base ring). Only the differential - and shift/difference algebras are fully determined by which_algebra and - base_var. */ - else if (!(is_diff_case(sa) || is_shift_case(sa)) - || !(is_diff_case(da) || is_shift_case(da)) - || od->base_var != rd->base_var) + return GR_SUCCESS; + } + + if (!(is_diff_case(sa) || is_shift_case(sa)) + || !(is_diff_case(da) || is_shift_case(da))) return GR_UNABLE; - else if (sa == da) + + if (od->base_var != rd->base_var) + return GR_UNABLE; + + if (sa == da) status |= _gr_vec_set(res->coeffs, op->coeffs, len, base); else if (sa == ORE_ALGEBRA_DERIVATIVE && da == ORE_ALGEBRA_EULER_DERIVATIVE) { diff --git a/src/gr_ore_poly/ctx.c b/src/gr_ore_poly/ctx.c index 3dc997aa9c..42bdb7493d 100644 --- a/src/gr_ore_poly/ctx.c +++ b/src/gr_ore_poly/ctx.c @@ -416,9 +416,6 @@ gr_ore_poly_ctx_init_randtest(gr_ore_poly_ctx_t ctx, flint_rand_t state, if (gr_ctx_ngens(&ngens, base_ring) != GR_SUCCESS) ngens = 0; - /* an invalid index when the base ring has no generator, so that a spurious - reliance on a default generator surfaces as GR_UNABLE rather than silently - using generator 0 */ base_var = (ngens > 0) ? (slong) n_randint(state, ngens) : -1; if (alg == ORE_ALGEBRA_CUSTOM) diff --git a/src/gr_ore_poly/test/t-convert.c b/src/gr_ore_poly/test/t-convert.c index 282fe765bb..13ac1ddf28 100644 --- a/src/gr_ore_poly/test/t-convert.c +++ b/src/gr_ore_poly/test/t-convert.c @@ -23,11 +23,8 @@ TEST_GR_FUNCTION_START(gr_ore_poly_convert, state, count_success, count_domain, gr_ore_poly_ctx_t ctx_src, ctx_dst; ore_algebra_t sa, da; slong power = 0, j; - int independent, cstatus, status = GR_SUCCESS; + int independent, status = GR_SUCCESS; - /* Most of the time convert between two algebras over a shared base ring; - occasionally use unrelated base rings to check that the function copes - sensibly (it returns GR_UNABLE) rather than misbehaving. */ independent = (n_randint(state, 4) == 0); gr_ore_poly_ctx_init_randtest2(cctx, ctx_src, state); @@ -45,12 +42,8 @@ TEST_GR_FUNCTION_START(gr_ore_poly_convert, state, count_success, count_domain, status |= gr_ore_poly_randtest(op, state, 1 + n_randint(state, 5), ctx_src); - cstatus = gr_ore_poly_convert(res, &power, op, ctx_dst, ctx_src); - status |= cstatus; + status = gr_ore_poly_convert(res, &power, op, ctx_dst, ctx_src); - /* Over these representative base rings the conversion never fails, as - long as both algebras stay within one family (differential or - shift/difference) and the two contexts share the base ring. */ int sa_diff = (sa == ORE_ALGEBRA_DERIVATIVE || sa == ORE_ALGEBRA_EULER_DERIVATIVE); int da_diff = (da == ORE_ALGEBRA_DERIVATIVE || da == ORE_ALGEBRA_EULER_DERIVATIVE); int sa_shift = (sa == ORE_ALGEBRA_FORWARD_SHIFT || sa == ORE_ALGEBRA_BACKWARD_SHIFT @@ -62,72 +55,56 @@ TEST_GR_FUNCTION_START(gr_ore_poly_convert, state, count_success, count_domain, && (POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_FMPZ || POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_CC_ACB || POLYNOMIAL_ELEM_CTX(cctx)->which_ring == GR_CTX_NMOD)); - if (expect_success && cstatus != GR_SUCCESS) + if (expect_success && status != GR_SUCCESS) { - flint_printf("FAIL: convert unexpected failure\n"); + flint_printf("FAIL: unexpected failure\n"); flint_printf("sa = %d, da = %d\n", sa, da); flint_abort(); } - /* When the conversion succeeds, the converted operator must act on the - base ring consistently with the original, up to the reported left - power: x^power for the differential family, the forward shift S^power - (a Taylor shift) for the shift/difference family. */ - if (cstatus == GR_SUCCESS) - { - int differential = (sa == ORE_ALGEBRA_DERIVATIVE || sa == ORE_ALGEBRA_EULER_DERIVATIVE) - && (da == ORE_ALGEBRA_DERIVATIVE || da == ORE_ALGEBRA_EULER_DERIVATIVE); + /* the converted operator must act on the base ring consistently with + the original */ + if (status == GR_SUCCESS) + { gr_ptr f = gr_heap_init(cctx); gr_ptr g_src = gr_heap_init(cctx); gr_ptr g_dst = gr_heap_init(cctx); gr_ptr corrected = gr_heap_init(cctx); - for (j = 0; j < 2; j++) + for (j = 0; j < 4; j++) { - int s1, s2; - status |= gr_randtest(f, state, cctx); - s1 = gr_ore_poly_apply(g_src, op, f, ctx_src); - s2 = gr_ore_poly_apply(g_dst, res, f, ctx_dst); - status |= s1 | s2; + status |= gr_ore_poly_apply(g_src, op, f, ctx_src); + status |= gr_ore_poly_apply(g_dst, res, f, ctx_dst); - if (s1 != GR_SUCCESS || s2 != GR_SUCCESS) + if (status != GR_SUCCESS) continue; - /* power == 0 is the identity correction and is the only case - reachable over a non-polynomial base ring (apply succeeds there - only for a constant operator, whose conversion has power 0). A - nonzero power implies a GR_CTX_GR_POLY base ring. */ if (power == 0) { status |= gr_set(corrected, g_src, cctx); } - else if (differential) + else if (sa_diff && da_diff) { - /* x^power * res = op, so corrected = x^{-power} * g_src, x the - generator of index var */ + /* for currently implemented conversions */ + FLINT_ASSERT(power <= 0); + slong var = GR_ORE_POLY_ORE_DATA(ctx_src)->base_var; gr_vec_t gens; - gr_ptr xp = gr_heap_init(cctx); gr_vec_init(gens, 0, cctx); status |= gr_gens(gens, cctx); if (var < gens->length) - status |= gr_set(xp, gr_vec_entry_srcptr(gens, var, cctx), cctx); + status |= gr_pow_si(corrected, gr_vec_entry_srcptr(gens, var, cctx), -power, cctx); else status |= GR_UNABLE; gr_vec_clear(gens, cctx); - status |= gr_pow_ui(xp, xp, (ulong) (-power), cctx); - status |= gr_mul(corrected, xp, g_src, cctx); - - gr_heap_clear(xp, cctx); + status |= gr_mul(corrected, corrected, g_src, cctx); } - else + else if (sa_shift && da_shift) { - /* S^power * res = op, so corrected = S^{-power}(g_src), a - Taylor shift of the base ring element (a polynomial) */ gr_ctx_struct * sctx = POLYNOMIAL_ELEM_CTX(cctx); gr_ptr c = gr_heap_init(sctx); @@ -139,7 +116,7 @@ TEST_GR_FUNCTION_START(gr_ore_poly_convert, state, count_success, count_domain, if (status == GR_SUCCESS && gr_equal(g_dst, corrected, cctx) == T_FALSE) { - flint_printf("FAIL: convert application inconsistent\n"); + flint_printf("FAIL: application\n"); flint_printf("sa = %d, da = %d, power = %wd\n", sa, da, power); flint_abort(); } diff --git a/src/gr_ore_poly/test/t-differential_shift.c b/src/gr_ore_poly/test/t-differential_shift.c index 2d4c398b0e..f1b6cdf496 100644 --- a/src/gr_ore_poly/test/t-differential_shift.c +++ b/src/gr_ore_poly/test/t-differential_shift.c @@ -56,9 +56,7 @@ check_scaled(const gr_ore_poly_t a, const gr_ore_poly_t b, slong net, } /* Returns 0 only if a == S^net * b is provably false, where S is the forward - shift operator. Since S is not an element of octx when it is a backward - algebra, we compare in an auxiliary algebra in which S is the generator: we - convert both operands to it and clear the S^p offsets returned by convert. */ + shift operator. */ static int check_shifted(const gr_ore_poly_t a, const gr_ore_poly_t b, slong net, gr_ore_poly_ctx_t octx) @@ -77,12 +75,9 @@ check_shifted(const gr_ore_poly_t a, const gr_ore_poly_t b, slong net, gr_ore_poly_init(lhs, fs_ctx); gr_ore_poly_init(rhs, fs_ctx); - /* S^ea * af = a and S^eb * bf = b */ status |= gr_ore_poly_convert(af, &ea, a, fs_ctx, octx); status |= gr_ore_poly_convert(bf, &eb, b, fs_ctx, octx); - /* a == S^net * b becomes S^(ea+c) af == S^(net+eb+c) bf, with c shifting - both exponents to be nonnegative so the powers of S are representable */ la = ea; lb = net + eb; c = FLINT_MAX(0, -FLINT_MIN(la, lb)); @@ -105,8 +100,8 @@ check_shifted(const gr_ore_poly_t a, const gr_ore_poly_t b, slong net, return ok; } -/* Applies a forward-shift recurrence R = sum_k q_k(nu) S^k to the coefficients - of a polynomial. */ +/* Applies a recurrence R = sum_k q_k(nu) S^k to the coefficients of a + polynomial. */ static int apply_forward_recurrence(gr_poly_t s, const gr_ore_poly_t R, const gr_poly_t f, gr_ctx_t cctx) { @@ -168,12 +163,12 @@ TEST_GR_FUNCTION_START(gr_ore_poly_differential_to_shift, state, count_success, diff_alg = ds_diff_algs[n_randint(state, 2)]; shift_alg = ds_shift_algs[n_randint(state, 4)]; - if (n_randint(state, 16)) + if (!n_randint(state, 16)) { diff_alg = ORE_ALGEBRA_COMMUTATIVE; expect_success = 0; } - if (n_randint(state, 16)) + if (!n_randint(state, 16)) { shift_alg = ORE_ALGEBRA_COMMUTATIVE; expect_success = 0; @@ -207,7 +202,7 @@ TEST_GR_FUNCTION_START(gr_ore_poly_differential_to_shift, state, count_success, count_domain += ((status & GR_DOMAIN) != 0); count_unable += ((status & GR_UNABLE) != 0); - /* round trip recovers x^(p1-p2) * op */ + /* round trip */ if (status == GR_SUCCESS && !check_scaled(op2, op, p1 - p2, ctx_d)) { @@ -291,12 +286,12 @@ TEST_GR_FUNCTION_START(gr_ore_poly_shift_to_differential, state, count_success, shift_alg = ds_shift_algs[n_randint(state, 4)]; - if (n_randint(state, 16)) + if (!n_randint(state, 16)) { diff_alg = ORE_ALGEBRA_COMMUTATIVE; expect_success = 0; } - if (n_randint(state, 16)) + if (!n_randint(state, 16)) { shift_alg = ORE_ALGEBRA_COMMUTATIVE; expect_success = 0; diff --git a/src/gr_ore_poly/test/t-shift_convert.c b/src/gr_ore_poly/test/t-shift_convert.c index a43b049963..7f83a7fe11 100644 --- a/src/gr_ore_poly/test/t-shift_convert.c +++ b/src/gr_ore_poly/test/t-shift_convert.c @@ -99,6 +99,7 @@ TEST_GR_FUNCTION_START(gr_ore_poly_shift_convert, state, count_success, count_do flint_abort(); } + /* not guaranteed by the documentation */ if (status == GR_SUCCESS && sh + sh2 != 0) { flint_printf("FAIL: round-trip power\n"); From 5cce146d2ba8746521b34c6b751f24d4e3282201 Mon Sep 17 00:00:00 2001 From: "Marc Mezzarobba (AI)" Date: Wed, 1 Jul 2026 13:50:09 +0200 Subject: [PATCH 7/7] gr_ore_poly: minor preexisting documentation issues Co-Authored-By: Claude Opus 4.8 --- doc/source/gr_ore_poly.rst | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/doc/source/gr_ore_poly.rst b/doc/source/gr_ore_poly.rst index 5bf5380b0a..5933a3da7c 100644 --- a/doc/source/gr_ore_poly.rst +++ b/doc/source/gr_ore_poly.rst @@ -57,7 +57,7 @@ Ore algebra types .. enumerator:: ORE_ALGEBRA_FORWARD_DIFFERENCE - Linear difference operator in the forward finite difference operator. + Linear difference operators in the forward finite difference. The endomorphism `\sigma` is the shift `x \mapsto x + 1` with respect to a generator `x` of the base ring, and the `\sigma`-derivation @@ -69,7 +69,7 @@ Ore algebra types .. enumerator:: ORE_ALGEBRA_BACKWARD_DIFFERENCE - Linear difference operator in the backward finite difference operator. + Linear difference operators in the backward finite difference. .. enumerator:: ORE_ALGEBRA_Q_SHIFT @@ -136,8 +136,7 @@ Context object methods specific initialization function is listed below. .. function:: int gr_ore_poly_ctx_init_q_shift(gr_ore_poly_ctx_t ctx, gr_ctx_t base_ring, slong base_var, gr_srcptr q) - int gr_ore_poly_ctx_init_q_difference(gr_ore_poly_ctx_t ctx, gr_ctx_t base_ring, slong base_var, gr_srcptr q) - int gr_ore_poly_ctx_init_mahler(gr_ore_poly_ctx_t ctx, gr_ctx_t base_ring, slong base_var, long mahler_base) + int gr_ore_poly_ctx_init_mahler(gr_ore_poly_ctx_t ctx, gr_ctx_t base_ring, slong base_var, slong mahler_base) Like :func:`gr_ore_poly_ctx_init` for predefined Ore polynomial types where `\sigma` and `\delta` depend on parameters.