diff --git a/doc/source/gr_ore_poly.rst b/doc/source/gr_ore_poly.rst index 47c4c2fb09..5933a3da7c 100644 --- a/doc/source/gr_ore_poly.rst +++ b/doc/source/gr_ore_poly.rst @@ -3,11 +3,6 @@ **gr_ore_poly.h** -- dense univariate Ore polynomials over generic rings =============================================================================== -.. note:: - - This module is under construction. Functionality is currently limited to - memory management, additive arithmetic, and multiplication. - 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 @@ -28,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. @@ -44,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. @@ -52,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. @@ -60,31 +55,31 @@ 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. + 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 `\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. + Linear difference operators in the backward finite difference. - .. 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. @@ -141,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. @@ -267,6 +261,113 @@ 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. + +Conversions +------------------------------------------------------------------------------- + +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. + +.. function:: int _gr_ore_poly_euler_to_ddx(gr_ptr res, gr_srcptr op, slong len, slong var, gr_ctx_t ctx) + + 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) + + 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 cross between the + forward side `S`, `S-1` and the backward side `S^{-1}`, `1-S^{-1}` + currently 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 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 cf189b74df..3a3799c431 100644 --- a/src/gr_ore_poly.h +++ b/src/gr_ore_poly.h @@ -247,6 +247,21 @@ 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); + +/* 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); +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); +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); +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/apply.c b/src/gr_ore_poly/apply.c new file mode 100644 index 0000000000..d377b492d9 --- /dev/null +++ b/src/gr_ore_poly/apply.c @@ -0,0 +1,99 @@ +/* + 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); /* todo: add underscore version? */ + + 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(g, g, ctx); + 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); + + 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: + 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/backshift_to_euler_univar.c b/src/gr_ore_poly/backshift_to_euler_univar.c new file mode 100644 index 0000000000..18105d0c7b --- /dev/null +++ b/src/gr_ore_poly/backshift_to_euler_univar.c @@ -0,0 +1,64 @@ +/* + 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" + +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 _k, rcoeffs; + + if (ctx->which_ring != GR_CTX_GR_POLY) + return GR_UNABLE; + + sctx = POLYNOMIAL_ELEM_CTX(ctx); + ssz = sctx->sizeof_elem; + + GR_TMP_INIT(_k, sctx); + GR_TMP_INIT_VEC(rcoeffs, len, ctx); + + for (k = 0; k < len; k++) + { + 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++) + { + 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 * rck = (gr_poly_struct *) GR_ENTRY(rcoeffs, k, bsz); + gr_ptr dst = GR_ENTRY(ri->coeffs, k, ssz); + if (i < rck->length) + status |= gr_set(dst, GR_ENTRY(rck->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(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 new file mode 100644 index 0000000000..f66d438a42 --- /dev/null +++ b/src/gr_ore_poly/convert.c @@ -0,0 +1,98 @@ +/* + 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_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) +{ + 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) +{ + 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; + + *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_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 (res_ctx == op_ctx) + { + status |= _gr_vec_set(res->coeffs, op->coeffs, len, base); + return GR_SUCCESS; + } + + if (!(is_diff_case(sa) || is_shift_case(sa)) + || !(is_diff_case(da) || is_shift_case(da))) + return GR_UNABLE; + + 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) + { + status |= _gr_ore_poly_ddx_to_euler(res->coeffs, op->coeffs, len, var, base); + *p = -(len - 1); + } + 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); + else + return GR_UNABLE; + + 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..42bdb7493d 100644 --- a/src/gr_ore_poly/ctx.c +++ b/src/gr_ore_poly/ctx.c @@ -409,11 +409,15 @@ 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; + 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..bc751f7cb4 --- /dev/null +++ b/src/gr_ore_poly/ddx_to_euler.c @@ -0,0 +1,68 @@ +/* + 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_mat.h" +#include "gr_vec.h" +#include "gr_ore_poly.h" + + +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 stirling; + gr_vec_t gens; + gr_ptr x, xp, c; + + if (len <= 0) + return GR_SUCCESS; + + 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) + 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--) + { + for (k = 0; k <= i; k++) + { + 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(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 new file mode 100644 index 0000000000..604f9b1b80 --- /dev/null +++ b/src/gr_ore_poly/differential_to_shift.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_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; + slong a, sp = 0; + gr_srcptr eul; + 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; + + slong len = op->length; + if (sa == ORE_ALGEBRA_DERIVATIVE) + { + 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) + { + eul = op->coeffs; + a = 0; + } + else + return GR_UNABLE; + + 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; + } + + gr_ore_poly_fit_length(res, rlen, res_ctx); + 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); + _gr_ore_poly_normalise(res, res_ctx); + *p = sp + a; + } + + 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 new file mode 100644 index 0000000000..5a15379187 --- /dev/null +++ b/src/gr_ore_poly/euler_to_backshift_univar.c @@ -0,0 +1,55 @@ +/* + 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" + +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 negk; + + if (ctx->which_ring != GR_CTX_GR_POLY) + return GR_UNABLE; + + sctx = POLYNOMIAL_ELEM_CTX(ctx); + ssz = sctx->sizeof_elem; + + GR_TMP_INIT(negk, 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(negk, -k, sctx); + status |= gr_poly_taylor_shift(rk, rk, negk, 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 new file mode 100644 index 0000000000..2237a7ed6d --- /dev/null +++ b/src/gr_ore_poly/euler_to_ddx.c @@ -0,0 +1,65 @@ +/* + 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_mat.h" +#include "gr_vec.h" +#include "gr_ore_poly.h" + +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 stirling; + gr_vec_t gens; + gr_ptr x, xp, t; + + if (len == 0) + return GR_SUCCESS; + + 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) + 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); + + /* 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(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(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 new file mode 100644 index 0000000000..33c2ab9237 --- /dev/null +++ b/src/gr_ore_poly/shift_convert.c @@ -0,0 +1,178 @@ +/* + 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 +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 +is_backward(ore_algebra_t a) +{ + return a == ORE_ALGEBRA_BACKWARD_SHIFT || a == ORE_ALGEBRA_BACKWARD_DIFFERENCE; +} + +static int +is_difference(ore_algebra_t a) +{ + return a == ORE_ALGEBRA_FORWARD_DIFFERENCE || a == ORE_ALGEBRA_BACKWARD_DIFFERENCE; +} + +/* 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; + fmpz_t t; + + fmpz_init(t); + + for (m = 0; m < len; m++) + { + 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++) + { + 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); + } + } + + fmpz_clear(t); + + 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); +} + +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 (!is_shift_case(src_alg) || !is_shift_case(dst_alg)) + return GR_UNABLE; + + if (len <= 0) + return GR_SUCCESS; + + status |= _gr_vec_set(res, op, len, ctx); + + if (len == 1 || src_alg == dst_alg) + return status; + + src_back = is_backward(src_alg); + dst_back = is_backward(dst_alg); + crossing = (src_back != dst_back); + s = dst_back ? -(len - 1) : (len - 1); + + if (crossing && len >= 2) /* cases requiring a Taylor shift */ + { + if (ctx->which_ring != GR_CTX_GR_POLY) + return GR_UNABLE; + else + FLINT_ASSERT(var == 0); + } + + if (is_difference(src_alg)) + status |= binomial_transform(res, len, src_back ? -1 : 1, ctx); + + if (crossing) + { + status |= taylor_shift_all(res, len, s, ctx); + reverse(res, len, ctx); + *p = -s; + } + + 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 new file mode 100644 index 0000000000..3dd3ea12c5 --- /dev/null +++ b/src/gr_ore_poly/shift_to_differential.c @@ -0,0 +1,86 @@ +/* + 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_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_ctx_t rec_ctx, eul_ctx; + gr_ore_poly_t rec, eul; + + 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_UNABLE; + + slong len = op->length; + slong var = GR_ORE_POLY_ORE_DATA(res_ctx)->base_var; + + 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; + } + + slong elen = 0; + for (slong k = 0; k < len; k++) + { + const gr_poly_struct * q = GR_ENTRY(rec->coeffs, k, base->sizeof_elem); + if (q->length > elen) + elen = q->length; + } + if (elen == 0) + { + status |= gr_ore_poly_zero(res, res_ctx); + goto cleanup; + } + + 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->coeffs, len, base); + else + { + 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); + +cleanup: + + 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 79e983f851..f1d2f39563 100644 --- a/src/gr_ore_poly/test/main.c +++ b/src/gr_ore_poly/test/main.c @@ -16,6 +16,12 @@ #include "t-sigma_delta.c" #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 ***************************************************/ @@ -25,7 +31,14 @@ 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), + 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_shift_to_differential), + TEST_FUNCTION(gr_ore_poly_convert) }; /* 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..af34738dd5 --- /dev/null +++ b/src/gr_ore_poly/test/t-apply.c @@ -0,0 +1,278 @@ +/* + 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: D(gen)\n\noctx = %{gr_ctx}\n", octx); + flint_abort(); + } + + gr_heap_clear(got, cctx); + gr_ore_poly_clear(D, octx); +} + +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); + + 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); + + { + 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_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); + + 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++) + { + gr_ctx_t cctx, ctx; + slong maxlen; + int status = GR_SUCCESS; + + gr_ore_poly_ctx_init_randtest2(cctx, ctx, state); + + 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); + 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_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); + 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(); + } + + 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(); + } + + 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(); + } + + 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(); + } + + 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(); + } + + 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(); + } + + 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); +} 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..13ac1ddf28 --- /dev/null +++ b/src/gr_ore_poly/test/t-convert.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_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, status = GR_SUCCESS; + + 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); + + status = gr_ore_poly_convert(res, &power, op, ctx_dst, ctx_src); + + 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 && status != GR_SUCCESS) + { + flint_printf("FAIL: unexpected failure\n"); + flint_printf("sa = %d, da = %d\n", sa, da); + flint_abort(); + } + + /* 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 < 4; j++) + { + status |= gr_randtest(f, state, cctx); + status |= gr_ore_poly_apply(g_src, op, f, ctx_src); + status |= gr_ore_poly_apply(g_dst, res, f, ctx_dst); + + if (status != GR_SUCCESS) + continue; + + if (power == 0) + { + status |= gr_set(corrected, g_src, cctx); + } + else if (sa_diff && da_diff) + { + /* for currently implemented conversions */ + FLINT_ASSERT(power <= 0); + + slong var = GR_ORE_POLY_ORE_DATA(ctx_src)->base_var; + gr_vec_t gens; + + gr_vec_init(gens, 0, cctx); + status |= gr_gens(gens, cctx); + if (var < gens->length) + 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_mul(corrected, corrected, g_src, cctx); + } + else if (sa_shift && da_shift) + { + 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: application\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..b1e6fec31e --- /dev/null +++ b/src/gr_ore_poly/test/t-ddx_to_euler.c @@ -0,0 +1,146 @@ +/* + 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; + + 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; + } + + 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); + + 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; + + 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); + + 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(); + } + + count_success += (status == GR_SUCCESS); + count_domain += ((status & GR_DOMAIN) != 0); + count_unable += ((status & GR_UNABLE) != 0); + + /* correcting factor */ + + 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); + + /* 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); + _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: round trip\n"); + flint_abort(); + } + + /* test action */ + + for (j = 0; j < 4; j++) + { + 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); + + status |= gr_mul(scaled, xpow, g_d, cctx); + if (status == GR_SUCCESS && gr_equal(g_e, scaled, cctx) == T_FALSE) + { + flint_printf("FAIL: application identity\n"); + flint_abort(); + } + } + + 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..f1b6cdf496 --- /dev/null +++ b/src/gr_ore_poly/test/t-differential_shift.c @@ -0,0 +1,377 @@ +/* + 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. */ +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); + + status |= gr_ore_poly_convert(af, &ea, a, fs_ctx, octx); + status |= gr_ore_poly_convert(bf, &eb, b, fs_ctx, octx); + + 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 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 */ + + 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-euler_to_ddx.c b/src/gr_ore_poly/test/t-euler_to_ddx.c new file mode 100644 index 0000000000..61317be147 --- /dev/null +++ b/src/gr_ore_poly/test/t-euler_to_ddx.c @@ -0,0 +1,146 @@ +/* + 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; + + 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; + } + + sz = cctx->sizeof_elem; + + ngens = 0; + status |= gr_ctx_ngens(&ngens, cctx); + 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); + + 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); + + 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; + + 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); + + 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); + + status |= gr_gens(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); + + 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: round trip\n"); + flint_abort(); + } + + /* test action */ + + for (j = 0; j < 4; j++) + { + 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); + + if (status == GR_SUCCESS && gr_equal(g_d, g_e, cctx) == T_FALSE) + { + flint_printf("FAIL: application identity\n"); + flint_abort(); + } + } + +cleanup: + 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_vec_clear(gens, cctx); + 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..7f83a7fe11 --- /dev/null +++ b/src/gr_ore_poly/test/t-shift_convert.c @@ -0,0 +1,152 @@ +/* + 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 +}; + +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++) + { + gr_ctx_t cctx; + gr_ore_poly_ctx_t ctx_src, ctx_dst; + ore_algebra_t src_alg, dst_alg; + slong len, j, sh, sh2, 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); + + 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); + 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, &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); + + 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: unexpected failure\n"); + flint_abort(); + } + + /* not guaranteed by the documentation */ + if (status == GR_SUCCESS && sh + sh2 != 0) + { + 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: round trip\n"); + flint_abort(); + } + + 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 g1 = gr_heap_init(cctx); + gr_ptr g2 = gr_heap_init(cctx); + gr_ptr c = gr_heap_init(sctx); + + status |= gr_set_si(c, -sh, sctx); + + for (j = 0; j < 4; j++) + { + 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 *) g2, (gr_poly_struct *) g2, c, sctx); + if (status == GR_SUCCESS && gr_equal(g1, g2, cctx) == T_FALSE) + { + flint_printf("FAIL: application identity\n"); + flint_abort(); + } + } + + gr_heap_clear(f, cctx); + gr_heap_clear(g1, cctx); + gr_heap_clear(g2, cctx); + gr_heap_clear(c, sctx); + } + + 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); +}