From 7997e6b213103f95bca45a1fa361e4e705fd00be Mon Sep 17 00:00:00 2001 From: Ricardo Buring Date: Sat, 12 Apr 2025 12:33:58 +0200 Subject: [PATCH 1/3] Add generic Ore polynomial module This contains only the basic structure so far, such as memory management, additive arithmetic, and multiplication from the left by an element of the base ring. --- CMakeLists.txt | 3 +- Makefile.in | 3 +- doc/source/gr_domains.rst | 11 + doc/source/gr_ore_poly.rst | 233 +++++++++++++++++++ doc/source/index.rst | 1 + doc/source/index_generic.rst | 1 + src/generic_files/io.c | 9 + src/gr.h | 1 + src/gr/ore_poly.c | 18 ++ src/gr_ore_poly.h | 218 +++++++++++++++++ src/gr_ore_poly/ctx.c | 261 +++++++++++++++++++++ src/gr_ore_poly/inlines.c | 3 + src/gr_ore_poly/ring.c | 385 +++++++++++++++++++++++++++++++ src/gr_ore_poly/test/main.c | 27 +++ src/gr_ore_poly/test/t-ring.c | 69 ++++++ src/gr_ore_poly/test/t-set_str.c | 85 +++++++ 16 files changed, 1326 insertions(+), 2 deletions(-) create mode 100644 doc/source/gr_ore_poly.rst create mode 100644 src/gr/ore_poly.c create mode 100644 src/gr_ore_poly.h create mode 100644 src/gr_ore_poly/ctx.c create mode 100644 src/gr_ore_poly/inlines.c create mode 100644 src/gr_ore_poly/ring.c create mode 100644 src/gr_ore_poly/test/main.c create mode 100644 src/gr_ore_poly/test/t-ring.c create mode 100644 src/gr_ore_poly/test/t-set_str.c diff --git a/CMakeLists.txt b/CMakeLists.txt index 916aea9cbd..13bb1dd6ae 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -261,7 +261,8 @@ set(_BUILD_DIRS acb_theta dirichlet bernoulli hypgeom gr gr_generic gr_vec gr_mat - gr_poly gr_mpoly gr_series gr_special + gr_poly gr_mpoly gr_ore_poly gr_series + gr_special calcium fmpz_mpoly_q diff --git a/Makefile.in b/Makefile.in index 5b153a1a16..5c6258a978 100644 --- a/Makefile.in +++ b/Makefile.in @@ -210,7 +210,8 @@ HEADER_DIRS := \ acb_theta dirichlet bernoulli hypgeom \ \ gr gr_generic gr_vec gr_mat \ - gr_poly gr_mpoly gr_series gr_special \ + gr_poly gr_mpoly gr_ore_poly gr_series \ + gr_special \ \ calcium \ fmpz_mpoly_q \ diff --git a/doc/source/gr_domains.rst b/doc/source/gr_domains.rst index 45b371a878..d1f3deda18 100644 --- a/doc/source/gr_domains.rst +++ b/doc/source/gr_domains.rst @@ -318,6 +318,17 @@ Polynomial rings with monomial ordering *ord*. Elements have type :type:`gr_mpoly_struct`. +Ore polynomials +------------------------------------------------------------------------------- + +.. function:: void gr_ctx_init_gr_ore_poly(gr_ctx_t ctx, gr_ctx_t base_ring, slong base_var, const ore_algebra_t which_algebra) + + Initializes *ctx* to a ring of densely represented Ore polynomials over the + given *base_ring*, with the choice of Ore algebra structure given by + *which_algebra*. The Ore algebra structure may refer to a distinguished + generator of *base_ring*; this will be the generator of index *base_var*. + Elements have type :type:`gr_ore_poly_struct`. + Power series ------------------------------------------------------------------------------- diff --git a/doc/source/gr_ore_poly.rst b/doc/source/gr_ore_poly.rst new file mode 100644 index 0000000000..8db8cb6d9e --- /dev/null +++ b/doc/source/gr_ore_poly.rst @@ -0,0 +1,233 @@ +.. _gr-ore-poly: + +**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 from the left + by an element of the base ring. + +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 +derivative or Euler derivative) is stored in the context object +:type:`gr_ore_poly_ctx_t`. + +Most functions are provided in two versions: an underscore method which +operates directly on pre-allocated arrays of coefficients and generally +has some restrictions (often requiring the lengths to be nonzero +and not supporting aliasing of the input and output arrays), +and a non-underscore method which performs automatic memory +management and handles degenerate cases. + +Ore algebra types +-------------------------------------------------------------------------------- +.. type:: ore_algebra_t + + Represents one of the following supported Ore algebra types: + + .. macro:: ORE_ALGEBRA_DERIVATIVE + + The endomorphism `\sigma` is the identity, and the `\sigma`-derivation + `\delta` is the derivative `\frac{d}{dx}` with respect to a generator + `x` of the base ring. + + .. macro:: ORE_ALGEBRA_EULER_DERIVATIVE + + The endomorphism `\sigma` is the identity, and the `\sigma`-derivation + `\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 + + 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` is the zero map. + + .. macro:: ORE_ALGEBRA_FORWARD_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`. + +.. function:: ore_algebra_t ore_algebra_randtest(flint_rand_t state) + + Return a random Ore algebra type. + +Type compatibility +------------------------------------------------------------------------------- + +The ``gr_ore_poly`` type has the same data layout as ``gr_poly``. +Methods of ``gr_poly`` can therefore be used for linear and container +operations on a ``gr_ore_poly``, given that one is careful about providing +the right context object. + +Weak normalization +------------------------------------------------------------------------------- + +A :type:`gr_ore_poly_t` is always normalised by removing leading zeros. +For rings without decidable equality (e.g. rings with inexact +representation), only coefficients that are provably zero will be +removed, and there can thus be spurious leading zeros in the +internal representation. +Methods that depend on knowing the exact degree of a polynomial +will act appropriately, typically by returning ``GR_UNABLE`` +when it is unknown whether the leading stored coefficient is nonzero. + +Types, macros and constants +------------------------------------------------------------------------------- + +.. type:: gr_ore_poly_struct + +.. type:: gr_ore_poly_t + + Contains a pointer to an array of coefficients (``coeffs``), the used + length (``length``), and the allocated size of the array (``alloc``). + + A ``gr_ore_poly_t`` is defined as an array of length one of type + ``gr_ore_poly_struct``, permitting a ``gr_ore_poly_t`` to + be passed by reference. + +Context object methods +------------------------------------------------------------------------------- + +.. function:: void gr_ore_poly_ctx_init(gr_ore_poly_ctx_t ctx, gr_ctx_t base_ring, slong base_var, const ore_algebra_t which_algebra) + + Initializes ``ctx`` to a ring of densely represented Ore polynomials over + the given ``base_ring``, with the choice of Ore algebra structure given by + ``which_algebra``. The Ore algebra structure may refer to a distinguished + generator of ``base_ring``; this will be the generator of index + ``base_var``. + +.. function:: void gr_ore_poly_ctx_clear(gr_ore_poly_ctx_t ctx) + + Clears the context object ``ctx``. + +.. function:: void gr_ore_poly_ctx_init_rand(gr_ore_poly_ctx_t ctx, flint_rand_t state, gr_ctx_t base_ring) + + Initializes ``ctx`` with a random Ore algebra structure. + +The following methods implement parts of the standard interface +for ``gr`` context objects. + +.. function:: int _gr_ore_poly_ctx_set_gen_name(gr_ctx_t ctx, const char * s) + int _gr_ore_poly_ctx_set_gen_names(gr_ctx_t ctx, const char ** s) + + Sets the name of the generator to the string in ``s``, respectively the + first string in ``s``. + +.. function:: int gr_ore_poly_ctx_write(gr_stream_t out, gr_ore_poly_ctx_t ctx) + truth_t gr_ore_poly_ctx_is_ring(gr_ore_poly_ctx_t ctx) + truth_t gr_ore_poly_ctx_is_zero_ring(gr_ore_poly_ctx_t ctx) + truth_t gr_ore_poly_ctx_is_commutative_ring(gr_ore_poly_ctx_t ctx) + truth_t gr_ore_poly_ctx_is_integral_domain(gr_ore_poly_ctx_t ctx) + truth_t gr_ore_poly_ctx_is_threadsafe(gr_ore_poly_ctx_t ctx) + +Memory management +------------------------------------------------------------------------------- + +.. function:: void gr_ore_poly_init(gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx) + +.. function:: void gr_ore_poly_init2(gr_ore_poly_t poly, slong len, gr_ore_poly_ctx_t ctx) + +.. function:: void gr_ore_poly_clear(gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx) + +.. function:: gr_ptr gr_ore_poly_coeff_ptr(gr_ore_poly_t poly, slong i, gr_ore_poly_ctx_t ctx) + gr_srcptr gr_ore_poly_coeff_srcptr(const gr_ore_poly_t poly, slong i, gr_ore_poly_ctx_t ctx) + +.. function:: slong gr_ore_poly_length(const gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx) + +.. function:: void gr_ore_poly_swap(gr_ore_poly_t poly1, gr_ore_poly_t poly2, gr_ore_poly_ctx_t ctx) + +.. function:: void gr_ore_poly_fit_length(gr_ore_poly_t poly, slong len, gr_ore_poly_ctx_t ctx) + +.. function:: void _gr_ore_poly_set_length(gr_ore_poly_t poly, slong len, gr_ore_poly_ctx_t ctx) + +Basic manipulation +------------------------------------------------------------------------------- + +.. function:: void _gr_ore_poly_normalise(gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx) + +.. function:: int gr_ore_poly_set(gr_ore_poly_t res, const gr_ore_poly_t src, gr_ore_poly_ctx_t ctx) + +.. function:: int gr_ore_poly_truncate(gr_ore_poly_t res, const gr_ore_poly_t poly, slong newlen, gr_ore_poly_ctx_t ctx) + +.. function:: int gr_ore_poly_zero(gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx) + int gr_ore_poly_one(gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx) + int gr_ore_poly_neg_one(gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx) + int gr_ore_poly_gen(gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx) + +.. function:: int gr_ore_poly_write(gr_stream_t out, const gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx) + int _gr_ore_poly_write(gr_stream_t out, gr_srcptr poly, slong n, gr_ore_poly_ctx_t ctx) + int _gr_ore_poly_get_str(char ** res, const gr_ore_poly_t f, gr_ore_poly_ctx_t ctx) + int gr_ore_poly_get_str(char ** res, const gr_ore_poly_t f, gr_ore_poly_ctx_t ctx) + int gr_ore_poly_print(const gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx) + +.. function:: int _gr_ore_poly_set_str(gr_ptr res, const char * s, slong len, gr_ore_poly_ctx_t ctx) + int gr_ore_poly_set_str(gr_ore_poly_t res, const char * s, gr_ore_poly_ctx_t ctx) + + Parse Ore polynomial from an expression string, assuming the name of the + generator is the one set in *ctx*. The underscore method zero-pads the + result if the length of the parsed polynomial is shorter than *len*, + and returns ``GR_UNABLE`` if the length of the parsed polynomial exceeds + *len*. Intermediate terms are allowed to be longer than *len*. + +.. function:: int gr_ore_poly_randtest(gr_ore_poly_t poly, flint_rand_t state, slong len, gr_ore_poly_ctx_t ctx) + +.. function:: truth_t _gr_ore_poly_equal(gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, gr_ore_poly_ctx_t ctx) + truth_t gr_ore_poly_equal(const gr_ore_poly_t poly1, const gr_ore_poly_t poly2, gr_ore_poly_ctx_t ctx) + +.. function:: truth_t gr_ore_poly_is_zero(const gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx) + truth_t gr_ore_poly_is_one(const gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx) + truth_t gr_ore_poly_is_gen(const gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx) + +.. function:: int gr_ore_poly_set_si(gr_ore_poly_t poly, slong c, gr_ore_poly_ctx_t ctx) + int gr_ore_poly_set_ui(gr_ore_poly_t poly, ulong c, gr_ore_poly_ctx_t ctx) + int gr_ore_poly_set_fmpz(gr_ore_poly_t poly, const fmpz_t c, gr_ore_poly_ctx_t ctx) + int gr_ore_poly_set_fmpq(gr_ore_poly_t poly, const fmpq_t c, gr_ore_poly_ctx_t ctx) + int gr_ore_poly_set_other(gr_ore_poly_t poly, gr_srcptr x, gr_ctx_t x_ctx, gr_ore_poly_ctx_t ctx) + +Arithmetic +------------------------------------------------------------------------------- + +.. function:: int gr_ore_poly_neg(gr_ore_poly_t res, const gr_ore_poly_t src, gr_ore_poly_ctx_t ctx) + +.. function:: int _gr_ore_poly_add(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, gr_ore_poly_ctx_t ctx) + int gr_ore_poly_add(gr_ore_poly_t res, const gr_ore_poly_t poly1, const gr_ore_poly_t poly2, gr_ore_poly_ctx_t ctx) + +.. function:: int _gr_ore_poly_sub(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, gr_ore_poly_ctx_t ctx) + int gr_ore_poly_sub(gr_ore_poly_t res, const gr_ore_poly_t poly1, const gr_ore_poly_t poly2, gr_ore_poly_ctx_t ctx) + +.. function:: int gr_ore_poly_add_ui(gr_ore_poly_t res, const gr_ore_poly_t poly, ulong c, gr_ore_poly_ctx_t ctx) + int gr_ore_poly_add_si(gr_ore_poly_t res, const gr_ore_poly_t poly, slong c, gr_ore_poly_ctx_t ctx) + int gr_ore_poly_add_fmpz(gr_ore_poly_t res, const gr_ore_poly_t poly, const fmpz c, gr_ore_poly_ctx_t ctx) + int gr_ore_poly_add_fmpq(gr_ore_poly_t res, const gr_ore_poly_t poly, const fmpq c, gr_ore_poly_ctx_t ctx) + int gr_ore_poly_add_other(gr_ore_poly_t res, const gr_ore_poly_t poly, gr_srcptr x, gr_ctx_t x_ctx, gr_ore_poly_ctx_t ctx) + + Sets *res* to *poly* plus the scalar *c* which must be + an element of or coercible to the coefficient ring. + +.. function:: int gr_ore_poly_sub_ui(gr_ore_poly_t res, const gr_ore_poly_t poly, ulong c, gr_ore_poly_ctx_t ctx) + int gr_ore_poly_sub_si(gr_ore_poly_t res, const gr_ore_poly_t poly, slong c, gr_ore_poly_ctx_t ctx) + int gr_ore_poly_sub_fmpz(gr_ore_poly_t res, const gr_ore_poly_t poly, const fmpz c, gr_ore_poly_ctx_t ctx) + int gr_ore_poly_sub_fmpq(gr_ore_poly_t res, const gr_ore_poly_t poly, const fmpq c, gr_ore_poly_ctx_t ctx) + int gr_ore_poly_sub_other(gr_ore_poly_t res, const gr_ore_poly_t poly, gr_srcptr x, gr_ctx_t x_ctx, gr_ore_poly_ctx_t ctx) + + Sets *res* to *poly* minus *c* which must be + an element of or coercible to the coefficient ring. + +.. function:: int gr_ore_poly_mul_ui(gr_ore_poly_t res, const gr_ore_poly_t poly, ulong c, gr_ore_poly_ctx_t ctx) + int gr_ore_poly_mul_si(gr_ore_poly_t res, const gr_ore_poly_t poly, slong c, gr_ore_poly_ctx_t ctx) + int gr_ore_poly_mul_fmpz(gr_ore_poly_t res, const gr_ore_poly_t poly, const fmpz c, gr_ore_poly_ctx_t ctx) + int gr_ore_poly_mul_fmpq(gr_ore_poly_t res, const gr_ore_poly_t poly, const fmpq c, gr_ore_poly_ctx_t ctx) + int gr_ore_poly_other_mul(gr_ore_poly_t res, gr_srcptr x, gr_ctx_t x_ctx, const gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx) + + Sets *res* to *poly* multiplied by *c* (or *x* multiplied by *poly*) + which must be an element of or coercible to the coefficient ring. + +.. raw:: latex + + \newpage + diff --git a/doc/source/index.rst b/doc/source/index.rst index 94f8fa66da..f9cea12d61 100644 --- a/doc/source/index.rst +++ b/doc/source/index.rst @@ -60,6 +60,7 @@ Generic rings gr_mat.rst gr_poly.rst gr_mpoly.rst + gr_ore_poly.rst gr_series.rst .. only:: not latex diff --git a/doc/source/index_generic.rst b/doc/source/index_generic.rst index 8903dc74c6..ce80504476 100644 --- a/doc/source/index_generic.rst +++ b/doc/source/index_generic.rst @@ -16,5 +16,6 @@ gr_mat.rst gr_poly.rst gr_mpoly.rst + gr_ore_poly.rst gr_series.rst diff --git a/src/generic_files/io.c b/src/generic_files/io.c index f680ddc2d6..efb71bd6ba 100644 --- a/src/generic_files/io.c +++ b/src/generic_files/io.c @@ -26,6 +26,7 @@ #include "gr.h" #include "gr_vec.h" #include "gr_poly.h" +#include "gr_ore_poly.h" #include "gr_mat.h" int _gr_mat_write(gr_stream_t out, const gr_mat_t mat, int linebreaks, gr_ctx_t ctx); @@ -728,6 +729,14 @@ int flint_vfprintf(FILE * fs, const char * ip, va_list vlist) res += out->len; ip += STRING_LENGTH("gr_poly}"); } + else if (IS_FLINT_TYPE(ip, "gr_ore_poly")) + { + const gr_ore_poly_struct * elem = va_arg(vlist, const gr_ore_poly_struct *); + gr_ctx_struct * ctx = va_arg(vlist, gr_ctx_struct *); + GR_MUST_SUCCEED(gr_ore_poly_write(out, elem, ctx)); + res += out->len; + ip += STRING_LENGTH("gr_ore_poly}"); + } else if (IS_FLINT_TYPE(ip, "gr_mat")) { const gr_mat_struct * elem = va_arg(vlist, const gr_mat_struct *); diff --git a/src/gr.h b/src/gr.h index 008dceb010..c572314a4a 100644 --- a/src/gr.h +++ b/src/gr.h @@ -714,6 +714,7 @@ typedef enum GR_CTX_FMPZ_MPOLY, GR_CTX_GR_MPOLY, GR_CTX_FMPZ_MPOLY_Q, GR_CTX_FMPZ_MOD_MPOLY_Q, + GR_CTX_GR_ORE_POLY, GR_CTX_GR_FRACTION, GR_CTX_GR_SERIES, GR_CTX_SERIES_MOD_GR_POLY, GR_CTX_GR_MAT, diff --git a/src/gr/ore_poly.c b/src/gr/ore_poly.c new file mode 100644 index 0000000000..ee7a6f9ac8 --- /dev/null +++ b/src/gr/ore_poly.c @@ -0,0 +1,18 @@ +/* + Copyright (C) 2025 Ricardo Buring + + 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 "gr_ore_poly.h" + +void +gr_ctx_init_gr_ore_poly(gr_ctx_t ctx, gr_ctx_t base_ring, slong base_var, const ore_algebra_t which_algebra) +{ + gr_ore_poly_ctx_init(ctx, base_ring, base_var, which_algebra); +} diff --git a/src/gr_ore_poly.h b/src/gr_ore_poly.h new file mode 100644 index 0000000000..522fdd7a55 --- /dev/null +++ b/src/gr_ore_poly.h @@ -0,0 +1,218 @@ +/* + Copyright (C) 2025 Ricardo Buring + + 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 . +*/ + +#ifndef GR_ORE_POLY_H +#define GR_ORE_POLY_H + +#ifdef GR_ORE_POLY_INLINES_C +#define GR_ORE_POLY_INLINE +#else +#define GR_ORE_POLY_INLINE static inline +#endif + +#include "gr.h" + +#ifdef __cplusplus +extern "C" { +#endif + +typedef enum { + ORE_ALGEBRA_DERIVATIVE, + ORE_ALGEBRA_EULER_DERIVATIVE, + ORE_ALGEBRA_FORWARD_SHIFT, + ORE_ALGEBRA_FORWARD_DIFFERENCE, +} +ore_algebra_t; + +#define ORE_POLY_NUM_ALGEBRAS 2 + +GR_ORE_POLY_INLINE +ore_algebra_t ore_algebra_randtest(flint_rand_t state) +{ + return (ore_algebra_t) n_randint(state, ORE_POLY_NUM_ALGEBRAS); +} + +/* Compatible with gr_poly_struct */ +typedef struct +{ + gr_ptr coeffs; + slong alloc; + slong length; +} +gr_ore_poly_struct; + +typedef gr_ore_poly_struct gr_ore_poly_t[1]; + +/* Compatible with polynomial_ctx_t */ +typedef struct +{ + gr_ctx_struct * base_ring; + slong degree_limit; + char * var; + ore_algebra_t which_algebra; + slong base_var; + /* TODO: Function pointers to \sigma, \delta. */ +} +_gr_ore_poly_ctx_struct; + +typedef gr_ctx_struct gr_ore_poly_ctx_struct; + +typedef gr_ore_poly_ctx_struct gr_ore_poly_ctx_t[1]; + +#define GR_ORE_POLY_CTX(ring_ctx) ((_gr_ore_poly_ctx_struct *)((ring_ctx))) +#define GR_ORE_POLY_ELEM_CTX(ring_ctx) (GR_ORE_POLY_CTX(ring_ctx)->base_ring) + +/* Context object */ + +void gr_ctx_init_gr_ore_poly(gr_ctx_t ctx, gr_ctx_t base_ring, slong base_var, const ore_algebra_t which_algebra); + +void gr_ore_poly_ctx_init(gr_ore_poly_ctx_t ctx, gr_ctx_t base_ring, slong base_var, const ore_algebra_t which_algebra); +void gr_ore_poly_ctx_clear(gr_ore_poly_ctx_t ctx); + +void gr_ore_poly_ctx_init_rand(gr_ore_poly_ctx_t ctx, flint_rand_t state, gr_ctx_t base_ring); + +WARN_UNUSED_RESULT int _gr_ore_poly_ctx_set_gen_name(gr_ctx_t ctx, const char * s); +WARN_UNUSED_RESULT int _gr_ore_poly_ctx_set_gen_names(gr_ctx_t ctx, const char ** s); +WARN_UNUSED_RESULT int gr_ore_poly_gens_recursive(gr_vec_t vec, gr_ore_poly_ctx_t ctx); + +int gr_ore_poly_ctx_write(gr_stream_t out, gr_ore_poly_ctx_t ctx); + +truth_t gr_ore_poly_ctx_is_ring(gr_ore_poly_ctx_t ctx); +truth_t gr_ore_poly_ctx_is_zero_ring(gr_ore_poly_ctx_t ctx); +truth_t gr_ore_poly_ctx_is_commutative_ring(gr_ore_poly_ctx_t ctx); +truth_t gr_ore_poly_ctx_is_integral_domain(gr_ore_poly_ctx_t ctx); +truth_t gr_ore_poly_ctx_is_threadsafe(gr_ore_poly_ctx_t ctx); + +/* Memory management */ + +void gr_ore_poly_init(gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx); +void gr_ore_poly_init2(gr_ore_poly_t poly, slong len, gr_ore_poly_ctx_t ctx); +void gr_ore_poly_clear(gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx); + +GR_ORE_POLY_INLINE gr_ptr +gr_ore_poly_coeff_ptr(gr_ore_poly_t poly, slong i, gr_ore_poly_ctx_t ctx) +{ + return GR_ENTRY(poly->coeffs, i, GR_ORE_POLY_ELEM_CTX(ctx)->sizeof_elem); +} + +GR_ORE_POLY_INLINE gr_srcptr +gr_ore_poly_coeff_srcptr(const gr_ore_poly_t poly, slong i, gr_ore_poly_ctx_t ctx) +{ + return GR_ENTRY(poly->coeffs, i, GR_ORE_POLY_ELEM_CTX(ctx)->sizeof_elem); +} + +GR_ORE_POLY_INLINE slong gr_ore_poly_length(const gr_ore_poly_t poly, gr_ore_poly_ctx_t FLINT_UNUSED(ctx)) +{ + return poly->length; +} + +GR_ORE_POLY_INLINE void +gr_ore_poly_swap(gr_ore_poly_t poly1, gr_ore_poly_t poly2, gr_ore_poly_ctx_t FLINT_UNUSED(ctx)) +{ + FLINT_SWAP(gr_ore_poly_struct, *poly1, *poly2); +} + +GR_ORE_POLY_INLINE void +gr_ore_poly_set_shallow(gr_ore_poly_t res, const gr_ore_poly_t x, const gr_ore_poly_ctx_t ctx) +{ + *res = *x; +} + +void gr_ore_poly_fit_length(gr_ore_poly_t poly, slong len, gr_ore_poly_ctx_t ctx); +void _gr_ore_poly_set_length(gr_ore_poly_t poly, slong len, gr_ore_poly_ctx_t ctx); +void _gr_ore_poly_normalise(gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx); + +/* Basic manipulation */ + +WARN_UNUSED_RESULT int gr_ore_poly_set(gr_ore_poly_t res, const gr_ore_poly_t src, gr_ore_poly_ctx_t ctx); + +WARN_UNUSED_RESULT int gr_ore_poly_truncate(gr_ore_poly_t poly, const gr_ore_poly_t src, slong newlen, gr_ore_poly_ctx_t ctx); + +GR_ORE_POLY_INLINE WARN_UNUSED_RESULT int +gr_ore_poly_zero(gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx) +{ + _gr_ore_poly_set_length(poly, 0, ctx); + return GR_SUCCESS; +} + +WARN_UNUSED_RESULT int gr_ore_poly_one(gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx); +WARN_UNUSED_RESULT int gr_ore_poly_neg_one(gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx); +WARN_UNUSED_RESULT int gr_ore_poly_gen(gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx); + +truth_t _gr_ore_poly_equal(gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, gr_ore_poly_ctx_t ctx); +truth_t gr_ore_poly_equal(const gr_ore_poly_t poly1, const gr_ore_poly_t poly2, gr_ore_poly_ctx_t ctx); + +truth_t gr_ore_poly_is_zero(const gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx); +truth_t gr_ore_poly_is_one(const gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx); +truth_t gr_ore_poly_is_gen(const gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx); + +/* Input and output */ + +int _gr_ore_poly_write(gr_stream_t out, gr_srcptr poly, slong len, gr_ore_poly_ctx_t ctx); +int gr_ore_poly_write(gr_stream_t out, const gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx); +int _gr_ore_poly_get_str(char ** res, gr_srcptr f, slong len, gr_ore_poly_ctx_t ctx); +int gr_ore_poly_get_str(char ** res, const gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx); +int _gr_ore_poly_set_str(gr_ptr res, const char * s, slong len, gr_ore_poly_ctx_t ctx); +int gr_ore_poly_set_str(gr_ore_poly_t res, const char * s, gr_ore_poly_ctx_t ctx); +int gr_ore_poly_print(const gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx); + +/* Random generation */ + +int gr_ore_poly_randtest(gr_ore_poly_t poly, flint_rand_t state, slong len, gr_ore_poly_ctx_t ctx); + +GR_ORE_POLY_INLINE WARN_UNUSED_RESULT int +_gr_ore_poly_randtest_default(gr_ore_poly_t res, flint_rand_t state, gr_ore_poly_ctx_t ctx) +{ + return gr_ore_poly_randtest(res, state, n_randint(state, 5), ctx); +} + +/* Constants */ + +WARN_UNUSED_RESULT int gr_ore_poly_set_si(gr_ore_poly_t poly, slong x, gr_ore_poly_ctx_t ctx); +WARN_UNUSED_RESULT int gr_ore_poly_set_ui(gr_ore_poly_t poly, ulong x, gr_ore_poly_ctx_t ctx); +WARN_UNUSED_RESULT int gr_ore_poly_set_fmpz(gr_ore_poly_t poly, const fmpz_t x, gr_ore_poly_ctx_t ctx); +WARN_UNUSED_RESULT int gr_ore_poly_set_fmpq(gr_ore_poly_t poly, const fmpq_t x, gr_ore_poly_ctx_t ctx); +WARN_UNUSED_RESULT int gr_ore_poly_set_other(gr_ore_poly_t poly, gr_srcptr x, gr_ctx_t x_ctx, 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); + +WARN_UNUSED_RESULT int _gr_ore_poly_add(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, gr_ore_poly_ctx_t ctx); +WARN_UNUSED_RESULT int gr_ore_poly_add(gr_ore_poly_t res, const gr_ore_poly_t poly1, const gr_ore_poly_t poly2, gr_ore_poly_ctx_t ctx); + +WARN_UNUSED_RESULT int _gr_ore_poly_sub(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, gr_ore_poly_ctx_t ctx); +WARN_UNUSED_RESULT int gr_ore_poly_sub(gr_ore_poly_t res, const gr_ore_poly_t poly1, const gr_ore_poly_t poly2, gr_ore_poly_ctx_t ctx); + +WARN_UNUSED_RESULT int gr_ore_poly_add_ui(gr_ore_poly_t res, const gr_ore_poly_t poly, ulong c, gr_ore_poly_ctx_t ctx); +WARN_UNUSED_RESULT int gr_ore_poly_add_si(gr_ore_poly_t res, const gr_ore_poly_t poly, slong c, gr_ore_poly_ctx_t ctx); +WARN_UNUSED_RESULT int gr_ore_poly_add_fmpz(gr_ore_poly_t res, const gr_ore_poly_t poly, const fmpz_t c, gr_ore_poly_ctx_t ctx); +WARN_UNUSED_RESULT int gr_ore_poly_add_fmpq(gr_ore_poly_t res, const gr_ore_poly_t poly, const fmpq_t c, gr_ore_poly_ctx_t ctx); +WARN_UNUSED_RESULT int gr_ore_poly_add_other(gr_ore_poly_t res, const gr_ore_poly_t poly, gr_srcptr x, gr_ctx_t x_ctx, gr_ore_poly_ctx_t ctx); + +WARN_UNUSED_RESULT int gr_ore_poly_sub_ui(gr_ore_poly_t res, const gr_ore_poly_t poly, ulong c, gr_ore_poly_ctx_t ctx); +WARN_UNUSED_RESULT int gr_ore_poly_sub_si(gr_ore_poly_t res, const gr_ore_poly_t poly, slong c, gr_ore_poly_ctx_t ctx); +WARN_UNUSED_RESULT int gr_ore_poly_sub_fmpz(gr_ore_poly_t res, const gr_ore_poly_t poly, const fmpz_t c, gr_ore_poly_ctx_t ctx); +WARN_UNUSED_RESULT int gr_ore_poly_sub_fmpq(gr_ore_poly_t res, const gr_ore_poly_t poly, const fmpq_t c, gr_ore_poly_ctx_t ctx); +WARN_UNUSED_RESULT int gr_ore_poly_sub_other(gr_ore_poly_t res, const gr_ore_poly_t poly, gr_srcptr x, gr_ctx_t x_ctx, gr_ore_poly_ctx_t ctx); + +WARN_UNUSED_RESULT int gr_ore_poly_mul_ui(gr_ore_poly_t res, const gr_ore_poly_t poly, ulong c, gr_ore_poly_ctx_t ctx); +WARN_UNUSED_RESULT int gr_ore_poly_mul_si(gr_ore_poly_t res, const gr_ore_poly_t poly, slong c, gr_ore_poly_ctx_t ctx); +WARN_UNUSED_RESULT int gr_ore_poly_mul_fmpz(gr_ore_poly_t res, const gr_ore_poly_t poly, const fmpz_t c, gr_ore_poly_ctx_t ctx); +WARN_UNUSED_RESULT int gr_ore_poly_mul_fmpq(gr_ore_poly_t res, const gr_ore_poly_t poly, const fmpq_t c, gr_ore_poly_ctx_t ctx); +WARN_UNUSED_RESULT int gr_ore_poly_other_mul(gr_ore_poly_t res, gr_srcptr x, gr_ctx_t x_ctx, const gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx); + + +#ifdef __cplusplus +} +#endif + +#endif diff --git a/src/gr_ore_poly/ctx.c b/src/gr_ore_poly/ctx.c new file mode 100644 index 0000000000..d476c6445c --- /dev/null +++ b/src/gr_ore_poly/ctx.c @@ -0,0 +1,261 @@ +/* + Copyright (C) 2025 Ricardo Buring + + 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 . +*/ + +/* Ore polynomials over generic rings */ + +#include +#include +#include "gr_vec.h" +#include "gr_ore_poly.h" +#include "gr_generic.h" + +static const char * default_var = "D"; + +int gr_ore_poly_ctx_write(gr_stream_t out, gr_ore_poly_ctx_t ctx) +{ + gr_stream_write(out, "Ring of Ore polynomials over "); + gr_ctx_write(out, GR_ORE_POLY_ELEM_CTX(ctx)); + return GR_SUCCESS; +} + +int _gr_ore_poly_ctx_set_gen_name(gr_ctx_t ctx, const char * s) +{ + slong len; + len = strlen(s); + + if (GR_ORE_POLY_CTX(ctx)->var == default_var) + GR_ORE_POLY_CTX(ctx)->var = NULL; + + GR_ORE_POLY_CTX(ctx)->var = flint_realloc(GR_ORE_POLY_CTX(ctx)->var, len + 1); + memcpy(GR_ORE_POLY_CTX(ctx)->var, s, len + 1); + return GR_SUCCESS; +} + +int _gr_ore_poly_ctx_set_gen_names(gr_ctx_t ctx, const char ** s) +{ + return _gr_ore_poly_ctx_set_gen_name(ctx, s[0]); +} + +void +gr_ore_poly_ctx_clear(gr_ore_poly_ctx_t ctx) +{ + if (GR_ORE_POLY_CTX(ctx)->var != default_var) + { + flint_free(GR_ORE_POLY_CTX(ctx)->var); + } +} + +truth_t +gr_ore_poly_ctx_is_ring(gr_ore_poly_ctx_t ctx) +{ + return gr_ctx_is_ring(GR_ORE_POLY_ELEM_CTX(ctx)); +} + +truth_t +gr_ore_poly_ctx_is_zero_ring(gr_ore_poly_ctx_t ctx) +{ + return gr_ctx_is_zero_ring(GR_ORE_POLY_ELEM_CTX(ctx)); +} + +truth_t +gr_ore_poly_ctx_is_commutative_ring(gr_ore_poly_ctx_t ctx) +{ + return T_UNKNOWN; +} + +truth_t +gr_ore_poly_ctx_is_integral_domain(gr_ore_poly_ctx_t ctx) +{ + return T_UNKNOWN; +} + +truth_t +gr_ore_poly_ctx_is_unique_factorization_domain(gr_ore_poly_ctx_t ctx) +{ + return T_UNKNOWN; +} + +truth_t +gr_ore_poly_ctx_is_threadsafe(gr_ore_poly_ctx_t ctx) +{ + return gr_ctx_is_threadsafe(GR_ORE_POLY_ELEM_CTX(ctx)); +} + +int +ore_poly_write(gr_stream_t out, gr_ore_poly_t poly, gr_ctx_t ctx) +{ + /* todo */ + if (poly->length == 0) + { + gr_stream_write(out, "0"); + return GR_SUCCESS; + } + + return gr_ore_poly_write(out, poly, ctx); +} + +int +gr_ore_poly_i(gr_ore_poly_t res, gr_ore_poly_ctx_t ctx) +{ + int status; + gr_ore_poly_fit_length(res, 1, ctx); + _gr_ore_poly_set_length(res, 1, ctx); + status = gr_i(res->coeffs, GR_ORE_POLY_ELEM_CTX(ctx)); + _gr_ore_poly_normalise(res, ctx); + return status; +} + +int +gr_ore_poly_gens_recursive(gr_vec_t vec, gr_ore_poly_ctx_t ctx) +{ + int status; + gr_vec_t vec1; + slong i, n; + + /* Get generators of the element ring */ + gr_vec_init(vec1, 0, GR_ORE_POLY_ELEM_CTX(ctx)); + status = gr_gens_recursive(vec1, GR_ORE_POLY_ELEM_CTX(ctx)); + n = vec1->length; + + gr_vec_set_length(vec, n + 1, ctx); + + /* Promote to Ore polynomials */ + for (i = 0; i < n; i++) + status |= gr_ore_poly_set_other(gr_vec_entry_ptr(vec, i, ctx), + gr_vec_entry_srcptr(vec1, i, GR_ORE_POLY_ELEM_CTX(ctx)), + GR_ORE_POLY_ELEM_CTX(ctx), + ctx); + + status |= gr_ore_poly_gen(gr_vec_entry_ptr(vec, n, ctx), ctx); + + gr_vec_clear(vec1, GR_ORE_POLY_ELEM_CTX(ctx)); + + return status; +} + +int _gr_ore_poly_methods_initialized = 0; + +gr_static_method_table _gr_ore_poly_methods; + +gr_method_tab_input _gr_ore_poly_methods_input[] = +{ + {GR_METHOD_CTX_WRITE, (gr_funcptr) gr_ore_poly_ctx_write}, + {GR_METHOD_CTX_CLEAR, (gr_funcptr) gr_ore_poly_ctx_clear}, + + {GR_METHOD_CTX_IS_RING, (gr_funcptr) gr_ore_poly_ctx_is_ring}, + {GR_METHOD_CTX_IS_ZERO_RING, (gr_funcptr) gr_ore_poly_ctx_is_zero_ring}, + {GR_METHOD_CTX_IS_COMMUTATIVE_RING, (gr_funcptr) gr_ore_poly_ctx_is_commutative_ring}, + {GR_METHOD_CTX_IS_INTEGRAL_DOMAIN, (gr_funcptr) gr_ore_poly_ctx_is_integral_domain}, + {GR_METHOD_CTX_IS_UNIQUE_FACTORIZATION_DOMAIN, (gr_funcptr) gr_ore_poly_ctx_is_unique_factorization_domain}, + {GR_METHOD_CTX_IS_FIELD, (gr_funcptr) gr_generic_ctx_predicate_false}, + {GR_METHOD_CTX_IS_THREADSAFE, (gr_funcptr) gr_ore_poly_ctx_is_threadsafe}, + {GR_METHOD_CTX_SET_GEN_NAME, (gr_funcptr) _gr_ore_poly_ctx_set_gen_name}, + {GR_METHOD_CTX_SET_GEN_NAMES, (gr_funcptr) _gr_ore_poly_ctx_set_gen_names}, + + {GR_METHOD_INIT, (gr_funcptr) gr_ore_poly_init}, + {GR_METHOD_CLEAR, (gr_funcptr) gr_ore_poly_clear}, + {GR_METHOD_SWAP, (gr_funcptr) gr_ore_poly_swap}, + {GR_METHOD_SET_SHALLOW, (gr_funcptr) gr_ore_poly_set_shallow}, + {GR_METHOD_RANDTEST, (gr_funcptr) _gr_ore_poly_randtest_default}, + {GR_METHOD_WRITE, (gr_funcptr) ore_poly_write}, + {GR_METHOD_ZERO, (gr_funcptr) gr_ore_poly_zero}, + {GR_METHOD_ONE, (gr_funcptr) gr_ore_poly_one}, + {GR_METHOD_NEG_ONE, (gr_funcptr) gr_ore_poly_neg_one}, + {GR_METHOD_I, (gr_funcptr) gr_ore_poly_i}, + + {GR_METHOD_GEN, (gr_funcptr) gr_ore_poly_gen}, + {GR_METHOD_GENS, (gr_funcptr) gr_generic_gens_single}, + {GR_METHOD_GENS_RECURSIVE, (gr_funcptr) gr_ore_poly_gens_recursive}, + +/* + {GR_METHOD_IS_ZERO, (gr_funcptr) gr_ore_poly_is_zero}, + {GR_METHOD_IS_ONE, (gr_funcptr) gr_ore_poly_is_one}, + {GR_METHOD_IS_NEG_ONE, (gr_funcptr) gr_ore_poly_is_neg_one}, +*/ + {GR_METHOD_EQUAL, (gr_funcptr) gr_ore_poly_equal}, + {GR_METHOD_SET, (gr_funcptr) gr_ore_poly_set}, + {GR_METHOD_SET_UI, (gr_funcptr) gr_ore_poly_set_ui}, + {GR_METHOD_SET_SI, (gr_funcptr) gr_ore_poly_set_si}, + {GR_METHOD_SET_FMPZ, (gr_funcptr) gr_ore_poly_set_fmpz}, + {GR_METHOD_SET_FMPQ, (gr_funcptr) gr_ore_poly_set_fmpq}, + {GR_METHOD_SET_OTHER, (gr_funcptr) gr_ore_poly_set_other}, +/* + {GR_METHOD_SET_INTERVAL_MID_RAD, (gr_funcptr) gr_ore_poly_set_interval_mid_rad}, +*/ + {GR_METHOD_SET_STR, (gr_funcptr) gr_ore_poly_set_str}, + {GR_METHOD_NEG, (gr_funcptr) gr_ore_poly_neg}, + {GR_METHOD_ADD_UI, (gr_funcptr) gr_ore_poly_add_ui}, + {GR_METHOD_ADD_SI, (gr_funcptr) gr_ore_poly_add_si}, + {GR_METHOD_ADD_FMPZ, (gr_funcptr) gr_ore_poly_add_fmpz}, + {GR_METHOD_ADD_FMPQ, (gr_funcptr) gr_ore_poly_add_fmpq}, + {GR_METHOD_ADD, (gr_funcptr) gr_ore_poly_add}, + {GR_METHOD_ADD_OTHER, (gr_funcptr) gr_ore_poly_add_other}, + {GR_METHOD_SUB_UI, (gr_funcptr) gr_ore_poly_sub_ui}, + {GR_METHOD_SUB_SI, (gr_funcptr) gr_ore_poly_sub_si}, + {GR_METHOD_SUB_FMPZ, (gr_funcptr) gr_ore_poly_sub_fmpz}, + {GR_METHOD_SUB_FMPQ, (gr_funcptr) gr_ore_poly_sub_fmpq}, + {GR_METHOD_SUB, (gr_funcptr) gr_ore_poly_sub}, + {GR_METHOD_SUB_OTHER, (gr_funcptr) gr_ore_poly_sub_other}, +/* + {GR_METHOD_MUL, (gr_funcptr) gr_ore_poly_mul}, + {GR_METHOD_MUL_OTHER, (gr_funcptr) gr_ore_poly_mul_other}, +*/ + {GR_METHOD_OTHER_MUL, (gr_funcptr) gr_ore_poly_other_mul}, + {GR_METHOD_MUL_UI, (gr_funcptr) gr_ore_poly_mul_ui}, + {GR_METHOD_MUL_SI, (gr_funcptr) gr_ore_poly_mul_si}, + {GR_METHOD_MUL_FMPZ, (gr_funcptr) gr_ore_poly_mul_fmpz}, + {GR_METHOD_MUL_FMPQ, (gr_funcptr) gr_ore_poly_mul_fmpq}, +/* + {GR_METHOD_POW_UI, (gr_funcptr) gr_ore_poly_pow_ui}, + {GR_METHOD_POW_SI, (gr_funcptr) gr_ore_poly_pow_si}, + {GR_METHOD_POW_FMPZ, (gr_funcptr) gr_ore_poly_pow_fmpz}, + {GR_METHOD_DIV, (gr_funcptr) gr_ore_poly_div}, + {GR_METHOD_INV, (gr_funcptr) gr_ore_poly_inv}, + + {GR_METHOD_EUCLIDEAN_DIV, (gr_funcptr) gr_ore_poly_euclidean_div}, + {GR_METHOD_EUCLIDEAN_REM, (gr_funcptr) gr_ore_poly_euclidean_rem}, + {GR_METHOD_EUCLIDEAN_DIVREM, (gr_funcptr) gr_ore_poly_euclidean_divrem}, + + {GR_METHOD_GCD, (gr_funcptr) gr_ore_poly_gcd}, + + {GR_METHOD_FACTOR, (gr_funcptr) gr_ore_poly_factor}, +*/ + + {0, (gr_funcptr) NULL}, +}; + +void +gr_ore_poly_ctx_init(gr_ore_poly_ctx_t ctx, gr_ctx_t base_ring, slong base_var, const ore_algebra_t which_algebra) +{ + ctx->which_ring = GR_CTX_GR_ORE_POLY; + ctx->sizeof_elem = sizeof(gr_ore_poly_struct); + ctx->size_limit = WORD_MAX; + + GR_ORE_POLY_CTX(ctx)->base_ring = (gr_ctx_struct *) base_ring; + GR_ORE_POLY_CTX(ctx)->degree_limit = WORD_MAX; + GR_ORE_POLY_CTX(ctx)->var = (char *) default_var; + GR_ORE_POLY_CTX(ctx)->which_algebra = which_algebra; + GR_ORE_POLY_CTX(ctx)->base_var = base_var; + + ctx->methods = _gr_ore_poly_methods; + + if (!_gr_ore_poly_methods_initialized) + { + gr_method_tab_init(_gr_ore_poly_methods, _gr_ore_poly_methods_input); + _gr_ore_poly_methods_initialized = 1; + } +} + +void +gr_ore_poly_ctx_init_rand(gr_ore_poly_ctx_t ctx, flint_rand_t state, gr_ctx_t base_ring) +{ + gr_ore_poly_ctx_init(ctx, base_ring, 0, ore_algebra_randtest(state)); +} diff --git a/src/gr_ore_poly/inlines.c b/src/gr_ore_poly/inlines.c new file mode 100644 index 0000000000..c336ce14bd --- /dev/null +++ b/src/gr_ore_poly/inlines.c @@ -0,0 +1,3 @@ +#define GR_ORE_POLY_INLINES_C +#include "gr_ore_poly.h" +#undef GR_ORE_POLY_INLINES_C diff --git a/src/gr_ore_poly/ring.c b/src/gr_ore_poly/ring.c new file mode 100644 index 0000000000..098d1697cd --- /dev/null +++ b/src/gr_ore_poly/ring.c @@ -0,0 +1,385 @@ +/* + Copyright (C) 2025 Ricardo Buring + + 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 +#include "fmpq.h" +#include "gr_poly.h" +#include "gr_ore_poly.h" + +/* Memory management */ + +void gr_ore_poly_init(gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx) +{ + poly->coeffs = NULL; + poly->length = 0; + poly->alloc = 0; +} + +void +gr_ore_poly_init2(gr_ore_poly_t poly, slong len, gr_ore_poly_ctx_t ctx) +{ + gr_poly_init((gr_poly_struct *) poly, GR_ORE_POLY_ELEM_CTX(ctx)); + gr_poly_fit_length((gr_poly_struct *) poly, len, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +void gr_ore_poly_clear(gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx) +{ + gr_poly_clear((gr_poly_struct *) poly, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +void +_gr_ore_poly_set_length(gr_ore_poly_t poly, slong len, gr_ore_poly_ctx_t ctx) +{ + _gr_poly_set_length((gr_poly_struct *) poly, len, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +void +gr_ore_poly_fit_length(gr_ore_poly_t poly, slong len, gr_ore_poly_ctx_t ctx) +{ + gr_poly_fit_length((gr_poly_struct *) poly, len, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +/* Basic manipulation */ + +void +_gr_ore_poly_normalise(gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx) +{ + _gr_poly_normalise((gr_poly_struct *) poly, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +int +gr_ore_poly_set(gr_ore_poly_t res, const gr_ore_poly_t src, gr_ore_poly_ctx_t ctx) +{ + return gr_poly_set((gr_poly_struct *) res, (const gr_poly_struct *) src, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +int +gr_ore_poly_truncate(gr_ore_poly_t poly, const gr_ore_poly_t src, slong newlen, gr_ore_poly_ctx_t ctx) +{ + return gr_poly_truncate((gr_poly_struct *) poly, (gr_poly_struct *) src, newlen, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +int +gr_ore_poly_one(gr_ore_poly_t res, gr_ore_poly_ctx_t ctx) +{ + return gr_poly_one((gr_poly_struct *) res, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +int +gr_ore_poly_neg_one(gr_ore_poly_t res, gr_ore_poly_ctx_t ctx) +{ + return gr_poly_neg_one((gr_poly_struct *) res, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +int +gr_ore_poly_gen(gr_ore_poly_t res, gr_ore_poly_ctx_t ctx) +{ + return gr_poly_gen((gr_poly_struct *) res, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +int +_gr_ore_poly_write(gr_stream_t out, gr_srcptr poly, slong n, gr_ore_poly_ctx_t ctx) +{ + return _gr_poly_write(out, poly, n, GR_ORE_POLY_CTX(ctx)->var, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +int +gr_ore_poly_write(gr_stream_t out, const gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx) +{ + return _gr_poly_write(out, poly->coeffs, poly->length, GR_ORE_POLY_CTX(ctx)->var, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +int +_gr_ore_poly_get_str(char ** res, gr_srcptr f, slong len, gr_ore_poly_ctx_t ctx) +{ + return _gr_poly_get_str(res, f, len, GR_ORE_POLY_CTX(ctx)->var, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +int +gr_ore_poly_get_str(char ** res, const gr_ore_poly_t f, gr_ore_poly_ctx_t ctx) +{ + return _gr_poly_get_str(res, f->coeffs, f->length, GR_ORE_POLY_CTX(ctx)->var, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +int +gr_ore_poly_print(const gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx) +{ + gr_stream_t out; + gr_stream_init_file(out, stdout); + return gr_ore_poly_write(out, poly, ctx); +} + +int +gr_ore_poly_set_str(gr_ore_poly_t res, const char * s, gr_ore_poly_ctx_t ctx) +{ + return gr_poly_set_str((gr_poly_struct *) res, s, GR_ORE_POLY_CTX(ctx)->var, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +int +_gr_ore_poly_set_str(gr_ptr res, const char * s, slong len, gr_ore_poly_ctx_t ctx) +{ + return _gr_poly_set_str(res, s, GR_ORE_POLY_CTX(ctx)->var, len, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +int +gr_ore_poly_randtest(gr_ore_poly_t poly, flint_rand_t state, slong len, gr_ore_poly_ctx_t ctx) +{ + return gr_poly_randtest((gr_poly_struct *) poly, state, len, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +truth_t +_gr_ore_poly_equal(gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, gr_ore_poly_ctx_t ctx) +{ + return _gr_poly_equal(poly1, len1, poly2, len2, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +truth_t +gr_ore_poly_equal(const gr_ore_poly_t poly1, const gr_ore_poly_t poly2, gr_ore_poly_ctx_t ctx) +{ + return gr_poly_equal((const gr_poly_struct *) poly1, (const gr_poly_struct *) poly2, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +truth_t +gr_ore_poly_is_zero(const gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx) +{ + return gr_poly_is_zero((const gr_poly_struct *) poly, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +truth_t +gr_ore_poly_is_one(const gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx) +{ + return gr_poly_is_one((const gr_poly_struct *) poly, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +truth_t +gr_ore_poly_is_gen(const gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx) +{ + return gr_poly_is_gen((const gr_poly_struct *) poly, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +int +gr_ore_poly_set_si(gr_ore_poly_t poly, slong x, gr_ore_poly_ctx_t ctx) +{ + return gr_poly_set_si((gr_poly_struct *) poly, x, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +int +gr_ore_poly_set_ui(gr_ore_poly_t poly, ulong x, gr_ore_poly_ctx_t ctx) +{ + return gr_poly_set_ui((gr_poly_struct *) poly, x, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +int +gr_ore_poly_set_fmpz(gr_ore_poly_t poly, const fmpz_t x, gr_ore_poly_ctx_t ctx) +{ + return gr_poly_set_fmpz((gr_poly_struct *) poly, x, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +int +gr_ore_poly_set_fmpq(gr_ore_poly_t poly, const fmpq_t x, gr_ore_poly_ctx_t ctx) +{ + return gr_poly_set_fmpq((gr_poly_struct *) poly, x, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +int +gr_ore_poly_set_other(gr_ore_poly_t poly, gr_srcptr x, gr_ctx_t x_ctx, gr_ore_poly_ctx_t ctx) +{ + if (x_ctx == ctx) { + return gr_poly_set((gr_poly_struct *) poly, x, GR_ORE_POLY_ELEM_CTX(ctx)); + } else if (x_ctx == GR_ORE_POLY_ELEM_CTX(ctx)) { + return gr_poly_set_scalar((gr_poly_struct *) poly, x, GR_ORE_POLY_ELEM_CTX(ctx)); + } else if (GR_ORE_POLY_CTX(ctx)->base_ring->which_ring == GR_CTX_GR_POLY && \ + POLYNOMIAL_CTX(GR_ORE_POLY_CTX(ctx)->base_ring)->base_ring == x_ctx) { + gr_poly_t tmp; + tmp->coeffs = (gr_ptr) x; + tmp->length = 1; + tmp->alloc = 1; + return gr_poly_set_scalar((gr_poly_struct *) poly, tmp, GR_ORE_POLY_ELEM_CTX(ctx)); + } else if (x_ctx->which_ring == GR_CTX_FMPZ) { + return gr_poly_set_fmpz((gr_poly_struct *) poly, x, GR_ORE_POLY_ELEM_CTX(ctx)); + } else if (x_ctx->which_ring == GR_CTX_FMPQ) { + return gr_poly_set_fmpq((gr_poly_struct *) poly, x, GR_ORE_POLY_ELEM_CTX(ctx)); + } + return GR_UNABLE; +} + +/* Arithmetic */ + +int +gr_ore_poly_neg(gr_ore_poly_t res, const gr_ore_poly_t src, gr_ore_poly_ctx_t ctx) +{ + return gr_poly_neg((gr_poly_struct *) res, (const gr_poly_struct *) src, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +int +_gr_ore_poly_add(gr_ptr res, gr_srcptr poly1, slong len1, + gr_srcptr poly2, slong len2, gr_ore_poly_ctx_t ctx) +{ + return _gr_poly_add(res, poly1, len1, poly2, len2, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +int +gr_ore_poly_add(gr_ore_poly_t res, const gr_ore_poly_t poly1, const gr_ore_poly_t poly2, gr_ore_poly_ctx_t ctx) +{ + return gr_poly_add((gr_poly_struct *) res, (const gr_poly_struct *) poly1, (const gr_poly_struct *) poly2, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +int +_gr_ore_poly_sub(gr_ptr res, gr_srcptr poly1, slong len1, + gr_srcptr poly2, slong len2, gr_ore_poly_ctx_t ctx) +{ + return _gr_poly_sub((gr_poly_struct *) res, poly1, len1, poly2, len2, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +int +gr_ore_poly_sub(gr_ore_poly_t res, const gr_ore_poly_t poly1, const gr_ore_poly_t poly2, gr_ore_poly_ctx_t ctx) +{ + return gr_poly_sub((gr_poly_struct *) res, (const gr_poly_struct *) poly1, (const gr_poly_struct *) poly2, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +int +gr_ore_poly_add_other(gr_ore_poly_t res, const gr_ore_poly_t poly, gr_srcptr x, gr_ctx_t x_ctx, gr_ore_poly_ctx_t ctx) +{ + if (x_ctx == ctx) { + return gr_poly_add((gr_poly_struct *) res, (const gr_poly_struct *) poly, x, GR_ORE_POLY_ELEM_CTX(ctx)); + } else if (x_ctx == GR_ORE_POLY_ELEM_CTX(ctx)) { + return gr_poly_add_scalar((gr_poly_struct *) res, (const gr_poly_struct *) poly, x, GR_ORE_POLY_ELEM_CTX(ctx)); + } else if (GR_ORE_POLY_CTX(ctx)->base_ring->which_ring == GR_CTX_GR_POLY && \ + POLYNOMIAL_CTX(GR_ORE_POLY_CTX(ctx)->base_ring)->base_ring == x_ctx) { + gr_poly_t tmp; + tmp->coeffs = (gr_ptr) x; + tmp->length = 1; + tmp->alloc = 1; + return gr_poly_add_scalar((gr_poly_struct *) poly, (gr_poly_struct *) poly, tmp, GR_ORE_POLY_ELEM_CTX(ctx)); + } else if (x_ctx->which_ring == GR_CTX_FMPZ) { + return gr_poly_add_fmpz((gr_poly_struct *) res, (const gr_poly_struct *) poly, x, GR_ORE_POLY_ELEM_CTX(ctx)); + } else if (x_ctx->which_ring == GR_CTX_FMPQ) { + return gr_poly_add_fmpq((gr_poly_struct *) res, (const gr_poly_struct *) poly, x, GR_ORE_POLY_ELEM_CTX(ctx)); + } + return GR_UNABLE; +} + +int +gr_ore_poly_add_ui(gr_ore_poly_t res, const gr_ore_poly_t poly, ulong c, gr_ore_poly_ctx_t ctx) +{ + return gr_poly_add_ui((gr_poly_struct *) res, (const gr_poly_struct *) poly, c, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +int +gr_ore_poly_add_si(gr_ore_poly_t res, const gr_ore_poly_t poly, slong c, gr_ore_poly_ctx_t ctx) +{ + return gr_poly_add_si((gr_poly_struct *) res, (const gr_poly_struct *) poly, c, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +int +gr_ore_poly_add_fmpz(gr_ore_poly_t res, const gr_ore_poly_t poly, const fmpz_t c, gr_ore_poly_ctx_t ctx) +{ + return gr_poly_add_fmpz((gr_poly_struct *) res, (const gr_poly_struct *) poly, c, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +int +gr_ore_poly_add_fmpq(gr_ore_poly_t res, const gr_ore_poly_t poly, const fmpq_t c, gr_ore_poly_ctx_t ctx) +{ + return gr_poly_add_fmpq((gr_poly_struct *) res, (const gr_poly_struct *) poly, c, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +int +gr_ore_poly_sub_other(gr_ore_poly_t res, const gr_ore_poly_t poly, gr_srcptr x, gr_ctx_t x_ctx, gr_ore_poly_ctx_t ctx) +{ + if (x_ctx == ctx) { + return gr_poly_sub((gr_poly_struct *) res, (const gr_poly_struct *) poly, x, GR_ORE_POLY_ELEM_CTX(ctx)); + } else if (x_ctx == GR_ORE_POLY_ELEM_CTX(ctx)) { + return gr_poly_sub_scalar((gr_poly_struct *) res, (gr_poly_struct *) poly, x, GR_ORE_POLY_ELEM_CTX(ctx)); + } else if (GR_ORE_POLY_CTX(ctx)->base_ring->which_ring == GR_CTX_GR_POLY && \ + POLYNOMIAL_CTX(GR_ORE_POLY_CTX(ctx)->base_ring)->base_ring == x_ctx) { + gr_poly_t tmp; + tmp->coeffs = (gr_ptr) x; + tmp->length = 1; + tmp->alloc = 1; + return gr_poly_sub_scalar((gr_poly_struct *) poly, (gr_poly_struct *) poly, tmp, GR_ORE_POLY_ELEM_CTX(ctx)); + } else if (x_ctx->which_ring == GR_CTX_FMPZ) { + return gr_poly_sub_fmpz((gr_poly_struct *) res, (const gr_poly_struct *) poly, x, GR_ORE_POLY_ELEM_CTX(ctx)); + } else if (x_ctx->which_ring == GR_CTX_FMPQ) { + return gr_poly_sub_fmpq((gr_poly_struct *) res, (const gr_poly_struct *) poly, x, GR_ORE_POLY_ELEM_CTX(ctx)); + } + return GR_UNABLE; +} + +int +gr_ore_poly_sub_ui(gr_ore_poly_t res, const gr_ore_poly_t poly, ulong c, gr_ore_poly_ctx_t ctx) +{ + return gr_poly_sub_ui((gr_poly_struct *) res, (const gr_poly_struct *) poly, c, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +int +gr_ore_poly_sub_si(gr_ore_poly_t res, const gr_ore_poly_t poly, slong c, gr_ore_poly_ctx_t ctx) +{ + return gr_poly_sub_si((gr_poly_struct *) res, (const gr_poly_struct *) poly, c, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +int +gr_ore_poly_sub_fmpz(gr_ore_poly_t res, const gr_ore_poly_t poly, const fmpz_t c, gr_ore_poly_ctx_t ctx) +{ + return gr_poly_sub_fmpz((gr_poly_struct *) res, (const gr_poly_struct *) poly, c, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +int +gr_ore_poly_sub_fmpq(gr_ore_poly_t res, const gr_ore_poly_t poly, const fmpq_t c, gr_ore_poly_ctx_t ctx) +{ + return gr_poly_sub_fmpq((gr_poly_struct *) res, (const gr_poly_struct *) poly, c, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +int +gr_ore_poly_other_mul(gr_ore_poly_t res, gr_srcptr x, gr_ctx_t x_ctx, const gr_ore_poly_t poly, gr_ore_poly_ctx_t ctx) +{ + if (x_ctx == ctx) { + /* TODO: Multiplication of Ore polynomials. */ + return GR_UNABLE; + } else if (x_ctx == GR_ORE_POLY_ELEM_CTX(ctx)) { + return gr_poly_scalar_mul((gr_poly_struct *) res, x, (const gr_poly_struct *) poly, GR_ORE_POLY_ELEM_CTX(ctx)); + } else if (GR_ORE_POLY_CTX(ctx)->base_ring->which_ring == GR_CTX_GR_POLY && \ + POLYNOMIAL_CTX(GR_ORE_POLY_CTX(ctx)->base_ring)->base_ring == x_ctx) { + gr_poly_t tmp; + tmp->coeffs = (gr_ptr) x; + tmp->length = 1; + tmp->alloc = 1; + return gr_poly_scalar_mul((gr_poly_struct *) poly, tmp, (const gr_poly_struct *) poly, GR_ORE_POLY_ELEM_CTX(ctx)); + } else if (x_ctx->which_ring == GR_CTX_FMPZ) { + return gr_poly_mul_fmpz((gr_poly_struct *) res, (const gr_poly_struct *) poly, x, GR_ORE_POLY_ELEM_CTX(ctx)); + } else if (x_ctx->which_ring == GR_CTX_FMPQ) { + return gr_poly_mul_fmpq((gr_poly_struct *) res, (const gr_poly_struct *) poly, x, GR_ORE_POLY_ELEM_CTX(ctx)); + } + return GR_UNABLE; +} + +int +gr_ore_poly_mul_ui(gr_ore_poly_t res, const gr_ore_poly_t poly, ulong c, gr_ore_poly_ctx_t ctx) +{ + return gr_poly_mul_ui((gr_poly_struct *) res, (const gr_poly_struct *) poly, c, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +int +gr_ore_poly_mul_si(gr_ore_poly_t res, const gr_ore_poly_t poly, slong c, gr_ore_poly_ctx_t ctx) +{ + return gr_poly_mul_si((gr_poly_struct *) res, (const gr_poly_struct *) poly, c, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +int +gr_ore_poly_mul_fmpz(gr_ore_poly_t res, const gr_ore_poly_t poly, const fmpz_t c, gr_ore_poly_ctx_t ctx) +{ + return gr_poly_mul_fmpz((gr_poly_struct *) res, (const gr_poly_struct *) poly, c, GR_ORE_POLY_ELEM_CTX(ctx)); +} + +int +gr_ore_poly_mul_fmpq(gr_ore_poly_t res, const gr_ore_poly_t poly, const fmpq_t c, gr_ore_poly_ctx_t ctx) +{ + return gr_poly_mul_fmpq((gr_poly_struct *) res, (const gr_poly_struct *) poly, c, GR_ORE_POLY_ELEM_CTX(ctx)); +} diff --git a/src/gr_ore_poly/test/main.c b/src/gr_ore_poly/test/main.c new file mode 100644 index 0000000000..a5ecd5293f --- /dev/null +++ b/src/gr_ore_poly/test/main.c @@ -0,0 +1,27 @@ +/* + Copyright (C) 2025 Ricardo Buring + + 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 functions *********************************************************/ + +#include "t-ring.c" +#include "t-set_str.c" + +/* Array of test functions ***************************************************/ + +test_struct tests[] = +{ + TEST_FUNCTION(gr_ore_poly_ring), + TEST_FUNCTION(gr_ore_poly_set_str), +}; + +/* main function *************************************************************/ + +TEST_MAIN(tests) diff --git a/src/gr_ore_poly/test/t-ring.c b/src/gr_ore_poly/test/t-ring.c new file mode 100644 index 0000000000..4f5666507f --- /dev/null +++ b/src/gr_ore_poly/test/t-ring.c @@ -0,0 +1,69 @@ +/* + Copyright (C) 2025 Ricardo Buring + + This file is part of FLINT. + + FLINT is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 3 of the License, or + (at your option) any later version. See . +*/ + +#include "test_helpers.h" +#include "gr_vec.h" +#include "gr_ore_poly.h" + +/* This is mostly a copy of src/gr_mpoly/test/t-ring.c */ + +FLINT_DLL extern gr_static_method_table _ca_methods; + +TEST_FUNCTION_START(gr_ore_poly_ring, state) +{ + slong iter; + + for (iter = 0; iter < 30 * flint_test_multiplier(); iter++) + { + gr_ctx_t ctx; + gr_ore_poly_ctx_t ore_ctx; + slong reps; + + gr_ctx_init_random(ctx, state); + /* TODO: ensure ctx has a generator */ + gr_ore_poly_ctx_init_rand(ore_ctx, state, ctx); + + if (gr_ctx_is_finite(ctx) == T_TRUE || + gr_ctx_has_real_prec(ctx) == T_TRUE) + { + reps = 10; + } + else if (ctx->methods == _ca_methods) /* hack: slow */ + { + reps = 1; + } + else + { + reps = 3; + } + + /* Hack: for string conversion tests, make sure we don't have + overlapping generator names. */ + gr_vec_t vec; + gr_vec_init(vec, 0, ctx); + if (gr_gens_recursive(vec, ctx) == GR_SUCCESS) + { + const char * vars[] = { "DD" }; + + GR_MUST_SUCCEED(gr_ctx_set_gen_names(ore_ctx, vars)); + + } + gr_vec_clear(vec, ctx); + + /* gr_ctx_println(ore_ctx); */ + gr_test_ring(ore_ctx, reps, 0 * GR_TEST_VERBOSE); + + gr_ore_poly_ctx_clear(ore_ctx); + gr_ctx_clear(ctx); + } + + TEST_FUNCTION_END(state); +} diff --git a/src/gr_ore_poly/test/t-set_str.c b/src/gr_ore_poly/test/t-set_str.c new file mode 100644 index 0000000000..d8bb3b080f --- /dev/null +++ b/src/gr_ore_poly/test/t-set_str.c @@ -0,0 +1,85 @@ +/* + Copyright (C) 2025 Ricardo Buring + + This file is part of FLINT. + + FLINT is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 3 of the License, or + (at your option) any later version. See . +*/ + +#include "test_helpers.h" +#include "fmpz.h" +#include "gr_vec.h" +#include "gr_ore_poly.h" + +/* This is mostly a copy of src/gr_poly/test/t-set_str.c */ + +TEST_FUNCTION_START(gr_ore_poly_set_str, state) +{ + slong iter; + + for (iter = 0; iter < 1000 * flint_test_multiplier(); iter++) + { + int status = GR_SUCCESS; + gr_ctx_t ctx; + gr_ore_poly_ctx_t ore_ctx; + gr_ore_poly_t f, g; + slong glen; + char * s; + + gr_ctx_init_random(ctx, state); + /* TODO: ensure ctx has a generator */ + gr_ore_poly_ctx_init_rand(ore_ctx, state, ctx); + gr_ore_poly_init(f, ore_ctx); + gr_ore_poly_init(g, ore_ctx); + + status |= gr_ore_poly_randtest(f, state, 1 + n_randint(state, 10), ore_ctx); + status |= gr_ore_poly_randtest(g, state, 1 + n_randint(state, 10), ore_ctx); + + if (n_randint(state, 2)) + { + status |= gr_ore_poly_get_str(&s, f, ore_ctx); + status |= gr_ore_poly_set_str(g, s, ore_ctx); + glen = -1; + + if ((ctx->which_ring == GR_CTX_FMPZ || ctx->which_ring == GR_CTX_FMPQ) && + status != GR_SUCCESS) + { + flint_printf("FAIL: unable over ZZ or QQ\n\n"); + flint_printf("f = %{gr_ore_poly}\n\n", f, ctx); + flint_printf("s = %s\n\n", s); + flint_abort(); + } + } + else + { + status |= _gr_ore_poly_get_str(&s, f->coeffs, f->length, ore_ctx); + glen = n_randint(state, 10); + gr_ore_poly_fit_length(g, glen, ore_ctx); + status |= _gr_ore_poly_set_str(g->coeffs, s, glen, ore_ctx); + _gr_ore_poly_set_length(g, glen, ore_ctx); + _gr_ore_poly_normalise(g, ore_ctx); + } + + if (status == GR_SUCCESS && gr_ore_poly_equal(f, g, ore_ctx) == T_FALSE) + { + flint_printf("FAIL: get_str, set_str roundtrip\n\n"); + flint_printf("f = %{gr_ore_poly}\n\n", f, ore_ctx); + flint_printf("s = %s\n\n", s); + flint_printf("glen = %wd\n\n", glen); + flint_printf("g = %{gr_ore_poly}\n\n", g, ore_ctx); + flint_abort(); + } + + flint_free(s); + + gr_ore_poly_clear(f, ore_ctx); + gr_ore_poly_clear(g, ore_ctx); + gr_ctx_clear(ore_ctx); + gr_ctx_clear(ctx); + } + + TEST_FUNCTION_END(state); +} From b1af5149c4271b04ecd12f75cf5bc13c0dad3749 Mon Sep 17 00:00:00 2001 From: Marc Mezzarobba Date: Wed, 17 Sep 2025 20:50:24 +0200 Subject: [PATCH 2/3] add gr_ctx_init_random_{poly,mpoly,series} --- doc/source/gr_domains.rst | 12 ++++++++++ src/gr.h | 3 +++ src/gr/init_random.c | 49 +++++++++++++++++++++++++++++++++++++++ 3 files changed, 64 insertions(+) diff --git a/doc/source/gr_domains.rst b/doc/source/gr_domains.rst index d1f3deda18..d6677506da 100644 --- a/doc/source/gr_domains.rst +++ b/doc/source/gr_domains.rst @@ -304,6 +304,10 @@ Polynomial rings over the given *base_ring*. Elements have type :type:`gr_poly_struct`. +.. function:: void gr_ctx_init_random_poly(gr_ctx_t ctx, flint_rand_t state) + + Initializes *ctx* to a random univariate polynomial ring. + .. function:: void gr_ctx_init_fmpz_mpoly(gr_ctx_t ctx, slong nvars, const ordering_t ord) Initializes *ctx* to a ring of sparsely represented multivariate @@ -318,6 +322,10 @@ Polynomial rings with monomial ordering *ord*. Elements have type :type:`gr_mpoly_struct`. +.. function:: void gr_ctx_init_random_mpoly(gr_ctx_t ctx, flint_rand_t state) + + Initializes *ctx* to a random multivariate polynomial ring. + Ore polynomials ------------------------------------------------------------------------------- @@ -335,6 +343,10 @@ Power series See :func:`gr_series_ctx_init` and :func:`gr_series_mod_ctx_init` in :ref:`gr-series`. +.. function:: void gr_ctx_init_random_series(gr_ctx_t ctx, flint_rand_t state) + + Initializes *ctx* to a random power series ring. + Fraction fields ------------------------------------------------------------------------------- diff --git a/src/gr.h b/src/gr.h index c572314a4a..2b3d73232c 100644 --- a/src/gr.h +++ b/src/gr.h @@ -1355,6 +1355,9 @@ truth_t gr_generic_ctx_predicate_false(gr_ctx_t ctx); /* Some base rings */ void gr_ctx_init_random(gr_ctx_t ctx, flint_rand_t state); +void gr_ctx_init_random_poly(gr_ctx_t ctx, flint_rand_t state); +void gr_ctx_init_random_mpoly(gr_ctx_t ctx, flint_rand_t state); +void gr_ctx_init_random_series(gr_ctx_t ctx, flint_rand_t state); void gr_ctx_init_fmpz(gr_ctx_t ctx); void gr_ctx_init_fmpq(gr_ctx_t ctx); diff --git a/src/gr/init_random.c b/src/gr/init_random.c index 73f7042e05..876158bb42 100644 --- a/src/gr/init_random.c +++ b/src/gr/init_random.c @@ -268,3 +268,52 @@ void gr_ctx_init_random(gr_ctx_t ctx, flint_rand_t state) gr_ctx_println(ctx); */ } + +void +gr_ctx_init_random_poly(gr_ctx_t ctx, flint_rand_t state) +{ + switch (n_randint(state, 8)) + { + case 0: + gr_ctx_init_fmpz_poly(ctx); + break; + case 1: + gr_ctx_init_fmpq_poly(ctx); + break; + default: + gr_ctx_init_gr_poly(ctx, _gr_random_base_ring(state)); + break; + } +} + +void +gr_ctx_init_random_mpoly(gr_ctx_t ctx, flint_rand_t state) +{ + ordering_t ordering = mpoly_ordering_randtest(state); + + switch (n_randint(state, 2)) + { + case 0: + gr_ctx_init_gr_mpoly(ctx, _gr_random_base_ring(state), n_randint(state, 3), ordering); + break; + case 1: + gr_ctx_init_fmpz_mpoly(ctx, n_randint(state, 3), mpoly_ordering_randtest(state)); + break; + } +} + +void +gr_ctx_init_random_series(gr_ctx_t ctx, flint_rand_t state) +{ + gr_ctx_struct * base_ring = _gr_random_base_ring(state); + + switch (n_randint(state, 2)) + { + case 0: + gr_series_ctx_init(ctx, base_ring, n_randint(state, 6)); + break; + case 1: + gr_series_mod_ctx_init(ctx, base_ring, n_randint(state, 6)); + break; + } +} From 885b42195b96da9ac823efe577eb7d6120acfaa9 Mon Sep 17 00:00:00 2001 From: Marc Mezzarobba Date: Wed, 17 Sep 2025 21:10:21 +0200 Subject: [PATCH 3/3] gr_ore_poly: sigma, delta --- doc/source/gr_ore_poly.rst | 85 +++++++++- src/gr_ore_poly.h | 90 +++++++++-- src/gr_ore_poly/ctx.c | 224 +++++++++++++++++++++++++- src/gr_ore_poly/ring.c | 1 - src/gr_ore_poly/sigma_delta.c | 232 +++++++++++++++++++++++++++ src/gr_ore_poly/test/main.c | 2 + src/gr_ore_poly/test/t-ring.c | 4 +- src/gr_ore_poly/test/t-set_str.c | 6 +- src/gr_ore_poly/test/t-sigma_delta.c | 156 ++++++++++++++++++ src/gr_poly/compose.c | 25 ++- 10 files changed, 785 insertions(+), 40 deletions(-) create mode 100644 src/gr_ore_poly/sigma_delta.c create mode 100644 src/gr_ore_poly/test/t-sigma_delta.c diff --git a/doc/source/gr_ore_poly.rst b/doc/source/gr_ore_poly.rst index 8db8cb6d9e..c94cd27e55 100644 --- a/doc/source/gr_ore_poly.rst +++ b/doc/source/gr_ore_poly.rst @@ -24,34 +24,71 @@ management and handles degenerate cases. Ore algebra types -------------------------------------------------------------------------------- + .. type:: ore_algebra_t Represents one of the following supported Ore algebra types: + .. macro:: ORE_ALGEBRA_CUSTOM + + Custom Ore polynomials. + + .. macro:: ORE_ALGEBRA_COMMUTATIVE + + Standard polynomials. + .. macro:: ORE_ALGEBRA_DERIVATIVE + Linear differential operators in the standard derivative. + The endomorphism `\sigma` is the identity, and the `\sigma`-derivation `\delta` is the derivative `\frac{d}{dx}` with respect to a generator `x` of the base ring. .. macro:: ORE_ALGEBRA_EULER_DERIVATIVE + Linear differential operators in the Euler derivative. + The endomorphism `\sigma` is the identity, and the `\sigma`-derivation `\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 + Linear difference operators in the standard forward shift. + 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` is the zero map. .. macro:: ORE_ALGEBRA_FORWARD_DIFFERENCE + Linear difference operator in the forward finite difference operator. + 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 + + Linear difference operators in the standard backward shift. + + .. macro:: ORE_ALGEBRA_BACKWARD_DIFFERENCE + + Linear difference operator in the backward finite difference operator. + + .. macro:: ORE_ALGEBRA_Q_SHIFT + + Linear q-difference operators. + + .. macro:: ORE_ALGEBRA_MAHLER + + Linear Mahler operators. + + .. macro:: ORE_ALGEBRA_FROBENIUS + + Ore polynomials over a field twisted by the Frobenius endomorphism. + .. function:: ore_algebra_t ore_algebra_randtest(flint_rand_t state) Return a random Ore algebra type. @@ -101,14 +138,40 @@ Context object methods generator of ``base_ring``; this will be the generator of index ``base_var``. -.. function:: void gr_ore_poly_ctx_clear(gr_ore_poly_ctx_t ctx) + This function can be used with all Ore algebra types for which no more + specific initialization function is listed below. - Clears the context object ``ctx``. +.. 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) + + Like :func:`gr_ore_poly_ctx_init` for predefined Ore polynomial types where + `\sigma` and `\delta` depend on parameters. + +.. function:: void * gr_ore_poly_ctx_data_ptr(gr_ore_poly_ctx_t ctx) + +.. function:: void gr_ore_poly_ctx_init_custom(gr_ore_poly_ctx_t ctx, gr_ctx_t base_ring, const gr_ore_poly_sigma_delta_t sigma_delta, void * ore_data) -.. function:: void gr_ore_poly_ctx_init_rand(gr_ore_poly_ctx_t ctx, flint_rand_t state, gr_ctx_t base_ring) + Initializes ``ctx`` to a ring of densely represented Ore polynomials over + the given ``base_ring``, with a custom Ore algebra structure specified by a + pointer ``sigma_delta`` to an implementation of + :func:`gr_ore_poly_sigma_delta`. + The ``ore_data`` argument is accessible to ``sigma_delta`` as + ``gr_ore_poly_ctx_data_ptr(ctx)``. + +.. function:: void gr_ore_poly_ctx_init_randtest(gr_ore_poly_ctx_t ctx, flint_rand_t state, gr_ctx_t base_ring) Initializes ``ctx`` with a random Ore algebra structure. +.. function:: void gr_ore_poly_ctx_init_randtest2(gr_ctx_t base_ring, gr_ore_poly_ctx_t ctx, flint_rand_t state) + + Initializes ``ctx`` with a random Ore algebra structure over a random base + ring. + +.. function:: void gr_ore_poly_ctx_clear(gr_ore_poly_ctx_t ctx) + + Clears the context object ``ctx``. + The following methods implement parts of the standard interface for ``gr`` context objects. @@ -189,6 +252,22 @@ Basic manipulation int gr_ore_poly_set_fmpq(gr_ore_poly_t poly, const fmpq_t c, gr_ore_poly_ctx_t ctx) int gr_ore_poly_set_other(gr_ore_poly_t poly, gr_srcptr x, gr_ctx_t x_ctx, gr_ore_poly_ctx_t ctx) +Action +------------------------------------------------------------------------------- + +.. function:: int gr_ore_poly_sigma(gr_ptr res, gr_srcptr a, gr_ore_poly_ctx_t ctx) + int gr_ore_poly_delta(gr_ptr res, gr_srcptr a, gr_ore_poly_ctx_t ctx) + int gr_ore_poly_sigma_delta(gr_ptr sigma, gr_ptr delta, gr_srcptr a, gr_ore_poly_ctx_t ctx) + + Compute *σ(a)*, *δ(a)*, or both, where *a* is an element of the base ring. + In the *sigma_delta* variant, the output variables *sigma* or *delta* can be + `NULL`. + +.. type:: gr_ore_poly_sigma_delta_t + + A pointer to a function with the same specification as + :func:`gr_ore_poly_sigma_delta`. + Arithmetic ------------------------------------------------------------------------------- diff --git a/src/gr_ore_poly.h b/src/gr_ore_poly.h index 522fdd7a55..039aa9a939 100644 --- a/src/gr_ore_poly.h +++ b/src/gr_ore_poly.h @@ -1,5 +1,6 @@ /* Copyright (C) 2025 Ricardo Buring + Copyright (C) 2025 Marc Mezzarobba This file is part of FLINT. @@ -24,21 +25,27 @@ extern "C" { #endif +/* Ore polynomial rings */ + typedef enum { + ORE_ALGEBRA_CUSTOM, + ORE_ALGEBRA_COMMUTATIVE, ORE_ALGEBRA_DERIVATIVE, ORE_ALGEBRA_EULER_DERIVATIVE, ORE_ALGEBRA_FORWARD_SHIFT, ORE_ALGEBRA_FORWARD_DIFFERENCE, + ORE_ALGEBRA_BACKWARD_SHIFT, + ORE_ALGEBRA_BACKWARD_DIFFERENCE, + ORE_ALGEBRA_Q_SHIFT, + ORE_ALGEBRA_MAHLER, + /* todo: q-derivative? general substitution?... */ + ORE_ALGEBRA_FROBENIUS, + + ORE_POLY_NUM_ALGEBRAS } ore_algebra_t; -#define ORE_POLY_NUM_ALGEBRAS 2 - -GR_ORE_POLY_INLINE -ore_algebra_t ore_algebra_randtest(flint_rand_t state) -{ - return (ore_algebra_t) n_randint(state, ORE_POLY_NUM_ALGEBRAS); -} +ore_algebra_t ore_algebra_randtest(flint_rand_t state); /* Compatible with gr_poly_struct */ typedef struct @@ -51,6 +58,23 @@ gr_ore_poly_struct; typedef gr_ore_poly_struct gr_ore_poly_t[1]; +typedef gr_ctx_struct gr_ore_poly_ctx_struct; + +typedef gr_ore_poly_ctx_struct gr_ore_poly_ctx_t[1]; + +typedef int (* gr_ore_poly_sigma_delta_t) (gr_ptr, gr_ptr, gr_srcptr, gr_ore_poly_ctx_struct *); + +typedef struct +{ + slong base_var; + gr_ptr sigma_x; + union { + gr_ptr q; + slong mahler_base; + }; +} +gr_ore_poly_ore_data_t; + /* Compatible with polynomial_ctx_t */ typedef struct { @@ -58,17 +82,20 @@ typedef struct slong degree_limit; char * var; ore_algebra_t which_algebra; - slong base_var; - /* TODO: Function pointers to \sigma, \delta. */ + gr_ore_poly_sigma_delta_t sigma_delta; + void * ore_data; } _gr_ore_poly_ctx_struct; -typedef gr_ctx_struct gr_ore_poly_ctx_struct; - -typedef gr_ore_poly_ctx_struct gr_ore_poly_ctx_t[1]; - #define GR_ORE_POLY_CTX(ring_ctx) ((_gr_ore_poly_ctx_struct *)((ring_ctx))) #define GR_ORE_POLY_ELEM_CTX(ring_ctx) (GR_ORE_POLY_CTX(ring_ctx)->base_ring) +#define GR_ORE_POLY_ORE_DATA(ring_ctx) ((gr_ore_poly_ore_data_t *)(GR_ORE_POLY_CTX(ring_ctx)->ore_data)) + +GR_ORE_POLY_INLINE void * +gr_ore_poly_ctx_data_ptr(gr_ore_poly_ctx_t ctx) +{ + return GR_ORE_POLY_CTX(ctx)->ore_data; +} /* Context object */ @@ -77,7 +104,8 @@ void gr_ctx_init_gr_ore_poly(gr_ctx_t ctx, gr_ctx_t base_ring, slong base_var, c void gr_ore_poly_ctx_init(gr_ore_poly_ctx_t ctx, gr_ctx_t base_ring, slong base_var, const ore_algebra_t which_algebra); void gr_ore_poly_ctx_clear(gr_ore_poly_ctx_t ctx); -void gr_ore_poly_ctx_init_rand(gr_ore_poly_ctx_t ctx, flint_rand_t state, gr_ctx_t base_ring); +void gr_ore_poly_ctx_init_randtest(gr_ore_poly_ctx_t ctx, flint_rand_t state, gr_ctx_t base_ring); +void gr_ore_poly_ctx_init_randtest2(gr_ctx_t base_ring, gr_ore_poly_ctx_t ctx, flint_rand_t state); WARN_UNUSED_RESULT int _gr_ore_poly_ctx_set_gen_name(gr_ctx_t ctx, const char * s); WARN_UNUSED_RESULT int _gr_ore_poly_ctx_set_gen_names(gr_ctx_t ctx, const char ** s); @@ -182,6 +210,40 @@ WARN_UNUSED_RESULT int gr_ore_poly_set_fmpz(gr_ore_poly_t poly, const fmpz_t x, WARN_UNUSED_RESULT int gr_ore_poly_set_fmpq(gr_ore_poly_t poly, const fmpq_t x, gr_ore_poly_ctx_t ctx); WARN_UNUSED_RESULT int gr_ore_poly_set_other(gr_ore_poly_t poly, gr_srcptr x, gr_ctx_t x_ctx, gr_ore_poly_ctx_t ctx); +/* Action on the base ring */ + +WARN_UNUSED_RESULT int sigma_delta_unable(gr_ptr sigma, gr_ptr delta, gr_srcptr a, gr_ore_poly_ctx_struct * ctx); +WARN_UNUSED_RESULT int sigma_delta_commutative(gr_ptr sigma, gr_ptr delta, gr_srcptr a, gr_ore_poly_ctx_struct * ctx); +WARN_UNUSED_RESULT int sigma_delta_derivative(gr_ptr sigma, gr_ptr delta, gr_srcptr a, gr_ore_poly_ctx_struct * ctx); +WARN_UNUSED_RESULT int sigma_delta_euler_derivative(gr_ptr sigma, gr_ptr delta, gr_srcptr a, gr_ore_poly_ctx_struct * ctx); +WARN_UNUSED_RESULT int sigma_delta_forward_shift(gr_ptr sigma, gr_ptr delta, gr_srcptr a, gr_ore_poly_ctx_struct * ctx); +WARN_UNUSED_RESULT int sigma_delta_backward_shift(gr_ptr sigma, gr_ptr delta, gr_srcptr a, gr_ore_poly_ctx_struct * ctx); +WARN_UNUSED_RESULT int sigma_delta_forward_difference(gr_ptr sigma, gr_ptr delta, gr_srcptr a, gr_ore_poly_ctx_struct * ctx); +WARN_UNUSED_RESULT int sigma_delta_backward_difference(gr_ptr sigma, gr_ptr delta, gr_srcptr a, gr_ore_poly_ctx_struct * ctx); +WARN_UNUSED_RESULT int sigma_delta_q_shift(gr_ptr sigma, gr_ptr delta, gr_srcptr a, gr_ore_poly_ctx_struct * ctx); +WARN_UNUSED_RESULT int sigma_delta_mahler(gr_ptr sigma, gr_ptr delta, gr_srcptr a, gr_ore_poly_ctx_struct * ctx); +WARN_UNUSED_RESULT int sigma_delta_frobenius(gr_ptr sigma, gr_ptr delta, gr_srcptr a, gr_ore_poly_ctx_struct * ctx); + +GR_ORE_POLY_INLINE WARN_UNUSED_RESULT int +gr_ore_poly_sigma_delta(gr_ptr sigma, gr_ptr delta, gr_srcptr a, gr_ore_poly_ctx_t ctx) +{ + return GR_ORE_POLY_CTX(ctx)->sigma_delta(sigma, delta, a, ctx); +} + +GR_ORE_POLY_INLINE WARN_UNUSED_RESULT int +gr_ore_poly_sigma(gr_ptr res, gr_srcptr a, gr_ore_poly_ctx_t ctx) +{ + return GR_ORE_POLY_CTX(ctx)->sigma_delta(res, NULL, a, ctx); +} + +GR_ORE_POLY_INLINE WARN_UNUSED_RESULT int +gr_ore_poly_delta(gr_ptr res, gr_srcptr a, gr_ore_poly_ctx_t ctx) +{ + return GR_ORE_POLY_CTX(ctx)->sigma_delta(NULL, res, a, ctx); +} + +extern const gr_ore_poly_sigma_delta_t _gr_ore_poly_default_sigma_delta[]; + /* 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/ctx.c b/src/gr_ore_poly/ctx.c index d476c6445c..f14c269f1c 100644 --- a/src/gr_ore_poly/ctx.c +++ b/src/gr_ore_poly/ctx.c @@ -13,15 +13,60 @@ #include #include +#include "fmpz_mpoly_q.h" +#include "fmpz_mod_mpoly_q.h" #include "gr_vec.h" #include "gr_ore_poly.h" +#include "gr_poly.h" #include "gr_generic.h" +#include "long_extras.h" +#include "mpoly.h" +#include "ulong_extras.h" static const char * default_var = "D"; int gr_ore_poly_ctx_write(gr_stream_t out, gr_ore_poly_ctx_t ctx) { - gr_stream_write(out, "Ring of Ore polynomials over "); + switch (GR_ORE_POLY_CTX(ctx)->which_algebra) + { + case ORE_ALGEBRA_CUSTOM: + gr_stream_write(out, "Custom Ore polynomials"); + break; + case ORE_ALGEBRA_COMMUTATIVE: + gr_stream_write(out, "Commutative Ore polynomials"); + break; + case ORE_ALGEBRA_DERIVATIVE: + gr_stream_write(out, "Differential operators"); + break; + case ORE_ALGEBRA_EULER_DERIVATIVE: + gr_stream_write(out, "Differential operators (Euler derivative)"); + break; + case ORE_ALGEBRA_FORWARD_SHIFT: + gr_stream_write(out, "Difference operators (forward shift)"); + break; + case ORE_ALGEBRA_FORWARD_DIFFERENCE: + gr_stream_write(out, "Difference operators (forward differences)"); + break; + case ORE_ALGEBRA_BACKWARD_SHIFT: + gr_stream_write(out, "Difference operators (backward shift)"); + break; + case ORE_ALGEBRA_BACKWARD_DIFFERENCE: + gr_stream_write(out, "Difference operators (backward differences)"); + break; + case ORE_ALGEBRA_Q_SHIFT: + gr_stream_write(out, "q-difference operators (q-shift)"); + break; + case ORE_ALGEBRA_MAHLER: + gr_stream_write(out, "Mahler operators"); + break; + case ORE_ALGEBRA_FROBENIUS: + gr_stream_write(out, "Ore-Frobenius polynomials"); + break; + default: + gr_stream_write(out, "Ore polynomials"); + break; + } + gr_stream_write(out, " over "); gr_ctx_write(out, GR_ORE_POLY_ELEM_CTX(ctx)); return GR_SUCCESS; } @@ -47,10 +92,23 @@ int _gr_ore_poly_ctx_set_gen_names(gr_ctx_t ctx, const char ** s) void gr_ore_poly_ctx_clear(gr_ore_poly_ctx_t ctx) { - if (GR_ORE_POLY_CTX(ctx)->var != default_var) + switch (GR_ORE_POLY_CTX(ctx)->which_algebra) { - flint_free(GR_ORE_POLY_CTX(ctx)->var); + case ORE_ALGEBRA_Q_SHIFT: + gr_heap_clear(GR_ORE_POLY_ORE_DATA(ctx)->q, + GR_ORE_POLY_ELEM_CTX(ctx)); + gr_heap_clear(GR_ORE_POLY_ORE_DATA(ctx)->sigma_x, + GR_ORE_POLY_ELEM_CTX(ctx)); + break; + default: + ; } + + if (GR_ORE_POLY_CTX(ctx)->which_algebra != ORE_ALGEBRA_CUSTOM) + flint_free(GR_ORE_POLY_CTX(ctx)->ore_data); + + if (GR_ORE_POLY_CTX(ctx)->var != default_var) + flint_free(GR_ORE_POLY_CTX(ctx)->var); } truth_t @@ -233,7 +291,7 @@ gr_method_tab_input _gr_ore_poly_methods_input[] = }; void -gr_ore_poly_ctx_init(gr_ore_poly_ctx_t ctx, gr_ctx_t base_ring, slong base_var, const ore_algebra_t which_algebra) +_gr_ore_poly_ctx_init(gr_ore_poly_ctx_t ctx, gr_ctx_t base_ring, const ore_algebra_t which_algebra) { ctx->which_ring = GR_CTX_GR_ORE_POLY; ctx->sizeof_elem = sizeof(gr_ore_poly_struct); @@ -243,7 +301,7 @@ gr_ore_poly_ctx_init(gr_ore_poly_ctx_t ctx, gr_ctx_t base_ring, slong base_var, GR_ORE_POLY_CTX(ctx)->degree_limit = WORD_MAX; GR_ORE_POLY_CTX(ctx)->var = (char *) default_var; GR_ORE_POLY_CTX(ctx)->which_algebra = which_algebra; - GR_ORE_POLY_CTX(ctx)->base_var = base_var; + GR_ORE_POLY_CTX(ctx)->sigma_delta = _gr_ore_poly_default_sigma_delta[which_algebra]; ctx->methods = _gr_ore_poly_methods; @@ -255,7 +313,159 @@ gr_ore_poly_ctx_init(gr_ore_poly_ctx_t ctx, gr_ctx_t base_ring, slong base_var, } void -gr_ore_poly_ctx_init_rand(gr_ore_poly_ctx_t ctx, flint_rand_t state, gr_ctx_t base_ring) +gr_ore_poly_ctx_init(gr_ore_poly_ctx_t ctx, gr_ctx_t base_ring, slong base_var, const ore_algebra_t which_algebra) +{ + _gr_ore_poly_ctx_init(ctx, base_ring, which_algebra); + + /* todo: FLINT_ASSERT(base_var < gr_ctx_ngens(base_ring)) */ + + GR_ORE_POLY_CTX(ctx)->ore_data = flint_malloc(sizeof(gr_ore_poly_ore_data_t)); + GR_ORE_POLY_ORE_DATA(ctx)->base_var = base_var; +} + +void +gr_ore_poly_ctx_init_custom(gr_ore_poly_ctx_t ctx, gr_ctx_t base_ring, + const gr_ore_poly_sigma_delta_t sigma_delta, + void * ore_data) +{ + _gr_ore_poly_ctx_init(ctx, base_ring, ORE_ALGEBRA_CUSTOM); + GR_ORE_POLY_CTX(ctx)->sigma_delta = sigma_delta; + GR_ORE_POLY_CTX(ctx)->ore_data = ore_data; +} + +/* although we currently require base_ring to be a univariate polynomial ring + * and in this case q will typically be a constant, for multivariate polynomial + * rings one typically wants q to be an element of the polynomial ring, not its + * base ring */ +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) +{ + if (base_ring->which_ring != GR_CTX_GR_POLY) + return GR_UNABLE; + + int status = GR_SUCCESS; + + gr_ore_poly_ctx_init(ctx, base_ring, base_var, ORE_ALGEBRA_Q_SHIFT); + GR_ORE_POLY_ORE_DATA(ctx)->q = gr_heap_init(base_ring); + gr_ptr sigma_x = gr_heap_init(base_ring); + GR_ORE_POLY_ORE_DATA(ctx)->sigma_x = sigma_x; + + status |= gr_set(GR_ORE_POLY_ORE_DATA(ctx)->q, q, base_ring); + status |= gr_gen(sigma_x, base_ring); + status |= gr_mul(sigma_x, q, sigma_x, base_ring); + + return status; +} + +int +gr_ore_poly_ctx_init_mahler(gr_ore_poly_ctx_t ctx, gr_ctx_t base_ring, + slong base_var, long mahler_base) +{ + if (mahler_base == 0) + return GR_DOMAIN; + + gr_ore_poly_ctx_init(ctx, base_ring, base_var, ORE_ALGEBRA_MAHLER); + GR_ORE_POLY_ORE_DATA(ctx)->mahler_base = mahler_base; + + return GR_SUCCESS; +} + +ore_algebra_t +ore_algebra_randtest(flint_rand_t state) +{ + switch (n_randint(state, 3)) + { + case 0: + return ORE_ALGEBRA_DERIVATIVE; + case 1: + return ORE_ALGEBRA_FORWARD_SHIFT; + default: + return (ore_algebra_t) n_randint(state, ORE_POLY_NUM_ALGEBRAS); + } +} + +void +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; + + int status = GR_SUCCESS; + + if (alg == ORE_ALGEBRA_CUSTOM) + { + /* sigma_delta_commutative does not use ore_data */ + gr_ore_poly_ctx_init_custom(ctx, base_ring, sigma_delta_commutative, NULL); + } + else if (alg == ORE_ALGEBRA_Q_SHIFT) + { + if (base_ring->which_ring == GR_CTX_GR_POLY) + { + gr_ptr q; + GR_TMP_INIT(q, base_ring); + status |= gr_randtest_not_zero(q, state, base_ring); + if (status == GR_SUCCESS) + status |= gr_ore_poly_ctx_init_q_shift(ctx, base_ring, base_var, q); + GR_TMP_CLEAR(q, base_ring); + } + else + { + status = GR_UNABLE; + } + } + else if (alg == ORE_ALGEBRA_MAHLER) + { + slong b = 2 + n_randint(state, 3); + status |= gr_ore_poly_ctx_init_mahler(ctx, base_ring, base_var, b); + } + else + { + gr_ore_poly_ctx_init(ctx, base_ring, base_var, alg); + } + + if (status != GR_SUCCESS) + gr_ore_poly_ctx_init_custom(ctx, base_ring, sigma_delta_unable, NULL); +} + +void +gr_ore_poly_ctx_init_randtest2(gr_ctx_t base_ring, gr_ore_poly_ctx_t ctx, flint_rand_t state) { - gr_ore_poly_ctx_init(ctx, base_ring, 0, ore_algebra_randtest(state)); + fmpz_t mod; + + switch (n_randint(state, 10)) + { + case 0: + gr_ctx_init_random_mpoly(base_ring, state); + break; + case 1: + gr_ctx_init_random_series(base_ring, state); + break; + case 2: + gr_ctx_init_fmpz_mpoly_q(base_ring, n_randint(state, 3), + mpoly_ordering_randtest(state)); + break; + case 3: + fmpz_init(mod); + fmpz_set_ui(mod, n_randtest_prime(state, 1)); + gr_ctx_init_fmpz_mod_mpoly_q(base_ring, n_randint(state, 3), + mpoly_ordering_randtest(state), + mod); + fmpz_clear(mod); + break; + case 4: + fmpz_init(mod); + fmpz_set_ui(mod, n_randtest_prime(state, 1)); + gr_ctx_init_fq(base_ring, mod, 1 + n_randint(state, 3), NULL); + gr_ore_poly_ctx_init(ctx, base_ring, 0, ORE_ALGEBRA_FROBENIUS); + fmpz_clear(mod); + return; + default: + gr_ctx_init_random_poly(base_ring, state); + break; + } + + gr_ore_poly_ctx_init_randtest(ctx, state, base_ring); } diff --git a/src/gr_ore_poly/ring.c b/src/gr_ore_poly/ring.c index 098d1697cd..1af0f51611 100644 --- a/src/gr_ore_poly/ring.c +++ b/src/gr_ore_poly/ring.c @@ -10,7 +10,6 @@ */ #include -#include "fmpq.h" #include "gr_poly.h" #include "gr_ore_poly.h" diff --git a/src/gr_ore_poly/sigma_delta.c b/src/gr_ore_poly/sigma_delta.c new file mode 100644 index 0000000000..0a6b93a14d --- /dev/null +++ b/src/gr_ore_poly/sigma_delta.c @@ -0,0 +1,232 @@ +#include "gr_poly.h" +#include "gr_ore_poly.h" + +/* gr_poly/compose.c */ +int _gr_poly_inflate(gr_ptr poly, slong len, slong n, gr_ctx_t ctx); + +int +sigma_delta_unable(gr_ptr sigma, gr_ptr delta, + gr_srcptr a, + gr_ore_poly_ctx_struct * ctx) +{ + return GR_UNABLE; +} + +int +sigma_delta_commutative(gr_ptr sigma, gr_ptr delta, + gr_srcptr a, + gr_ore_poly_ctx_struct * ctx) +{ + gr_ctx_ptr cctx = GR_ORE_POLY_ELEM_CTX(ctx); + int status = GR_SUCCESS; + + if (sigma != NULL && sigma != a) + status |= gr_set(sigma, a, cctx); + + if (delta != NULL) + status |= gr_zero(delta, cctx); + + return status; +} + +int +sigma_delta_derivative(gr_ptr sigma, gr_ptr delta, + gr_srcptr a, + gr_ore_poly_ctx_struct * ctx) +{ + gr_ctx_struct * cctx = GR_ORE_POLY_ELEM_CTX(ctx); + /* todo: gr_derivative taking a generator index? */ + if (cctx->which_ring != GR_CTX_GR_POLY) + return GR_UNABLE; + + int status = GR_SUCCESS; + + if (sigma != NULL && sigma != a) + status |= gr_set(sigma, a, cctx); + + if (delta != NULL) + status |= gr_poly_derivative(delta, a, POLYNOMIAL_ELEM_CTX(cctx)); + + return status; +} + +int +sigma_delta_euler_derivative(gr_ptr sigma, gr_ptr delta, + gr_srcptr a, + gr_ore_poly_ctx_struct * ctx) +{ + int status = GR_SUCCESS; + gr_ctx_struct * cctx = GR_ORE_POLY_ELEM_CTX(ctx); + if (cctx->which_ring != GR_CTX_GR_POLY) + return GR_UNABLE; + + /* - for polynomial base rings, could call _gr_poly_derivative and avoid the + * post-shift + * - gr_euler_derivative for arbitrary rings??? */ + + status |= sigma_delta_derivative(sigma, delta, a, ctx); + + if (status == GR_SUCCESS && delta != NULL) + status |= gr_poly_shift_left(delta, delta, 1, POLYNOMIAL_ELEM_CTX(cctx)); + + return status; +} + +static int +sigma_delta_shift_si(gr_ptr sigma, gr_ptr delta, + gr_srcptr a, slong shift, slong difference, + gr_ore_poly_ctx_struct * ctx) +{ + gr_ctx_struct * cctx = GR_ORE_POLY_ELEM_CTX(ctx); + if (cctx->which_ring != GR_CTX_GR_POLY) + return GR_UNABLE; + gr_ctx_struct * sctx = POLYNOMIAL_ELEM_CTX(cctx); + + int status = GR_SUCCESS; + + gr_ptr _shift; + gr_poly_t shifted; + GR_TMP_INIT(_shift, sctx); + gr_poly_init(shifted, sctx); + + status |= gr_set_si(_shift, shift, sctx); + status |= gr_poly_taylor_shift(shifted, a, _shift, sctx); + + if (delta != NULL) + { + if (difference == 0) + status |= gr_zero(delta, cctx); + else if (difference == 1) + status |= gr_poly_sub(delta, shifted, a, sctx); + else if (difference == -1) + status |= gr_poly_sub(delta, a, shifted, sctx); + } + + if (sigma != NULL) + gr_swap(sigma, shifted, cctx); + + gr_poly_clear(shifted, sctx); + GR_TMP_CLEAR(_shift, sctx); + + return status; +} + +int +sigma_delta_forward_shift(gr_ptr sigma, gr_ptr delta, + gr_srcptr a, + gr_ore_poly_ctx_struct * ctx) +{ + return sigma_delta_shift_si(sigma, delta, a, +1, 0, ctx); +} + +int +sigma_delta_backward_shift(gr_ptr sigma, gr_ptr delta, + gr_srcptr a, + gr_ore_poly_ctx_struct * ctx) +{ + return sigma_delta_shift_si(sigma, delta, a, -1, 0, ctx); +} + +int +sigma_delta_forward_difference(gr_ptr sigma, gr_ptr delta, + gr_srcptr a, + gr_ore_poly_ctx_struct * ctx) +{ + return sigma_delta_shift_si(sigma, delta, a, +1, +1, ctx); +} + +int +sigma_delta_backward_difference(gr_ptr sigma, gr_ptr delta, + gr_srcptr a, + gr_ore_poly_ctx_struct * ctx) +{ + return sigma_delta_shift_si(sigma, delta, a, -1, -1, ctx); +} + +int +sigma_delta_compose(gr_ptr sigma, gr_ptr delta, + gr_srcptr a, + gr_ore_poly_ctx_struct * ctx) +{ + gr_ctx_struct * cctx = GR_ORE_POLY_ELEM_CTX(ctx); + if (cctx->which_ring != GR_CTX_GR_POLY) + return GR_UNABLE; + gr_ctx_struct * sctx = POLYNOMIAL_ELEM_CTX(cctx); + + int status = GR_SUCCESS; + + if (sigma != NULL) + status |= gr_poly_compose(sigma, a, GR_ORE_POLY_ORE_DATA(ctx)->sigma_x, sctx); + + if (delta != NULL) + status |= gr_zero(delta, cctx); + + return status; +} + +int +sigma_delta_mahler(gr_ptr sigma, gr_ptr delta, + gr_srcptr a, + gr_ore_poly_ctx_struct * ctx) +{ + gr_ctx_struct * cctx = GR_ORE_POLY_ELEM_CTX(ctx); + if (cctx->which_ring != GR_CTX_GR_POLY) + return GR_UNABLE; + gr_ctx_struct * sctx = POLYNOMIAL_ELEM_CTX(cctx); + + int status = GR_SUCCESS; + + if (sigma != NULL) + { + gr_poly_struct * _sigma = sigma; + slong b = GR_ORE_POLY_ORE_DATA(ctx)->mahler_base; + status |= gr_set(sigma, a, cctx); + if (_sigma->length > 1) + { + slong newlen = (_sigma->length - 1) * b + 1; + gr_poly_fit_length(sigma, newlen, sctx); + status |= _gr_poly_inflate(_sigma->coeffs, _sigma->length, b, sctx); + gr_poly_fit_length(sigma, newlen, sctx); + _gr_poly_set_length(sigma, newlen, sctx); + _gr_poly_normalise(sigma, sctx); + } + } + + if (delta != NULL) + status |= gr_zero(delta, cctx); + + return status; +} + +int +sigma_delta_frobenius(gr_ptr sigma, gr_ptr delta, + gr_srcptr a, + gr_ore_poly_ctx_struct * ctx) +{ + gr_ctx_ptr cctx = GR_ORE_POLY_ELEM_CTX(ctx); + int status = GR_SUCCESS; + + if (sigma != NULL) + status |= gr_fq_frobenius(sigma, a, 1, cctx); + + if (delta != NULL) + status |= gr_zero(delta, cctx); + + return status; +} + +const gr_ore_poly_sigma_delta_t _gr_ore_poly_default_sigma_delta[] = +{ + [ORE_ALGEBRA_CUSTOM] = sigma_delta_unable, + [ORE_ALGEBRA_COMMUTATIVE] = sigma_delta_commutative, + [ORE_ALGEBRA_DERIVATIVE] = sigma_delta_derivative, + [ORE_ALGEBRA_EULER_DERIVATIVE] = sigma_delta_euler_derivative, + [ORE_ALGEBRA_FORWARD_SHIFT] = sigma_delta_forward_shift, + [ORE_ALGEBRA_FORWARD_DIFFERENCE] = sigma_delta_forward_difference, + [ORE_ALGEBRA_BACKWARD_SHIFT] = sigma_delta_backward_shift, + [ORE_ALGEBRA_BACKWARD_DIFFERENCE] = sigma_delta_backward_difference, + [ORE_ALGEBRA_Q_SHIFT] = sigma_delta_compose, + [ORE_ALGEBRA_MAHLER] = sigma_delta_mahler, + [ORE_ALGEBRA_FROBENIUS] = sigma_delta_frobenius, + [ORE_POLY_NUM_ALGEBRAS] = NULL, +}; diff --git a/src/gr_ore_poly/test/main.c b/src/gr_ore_poly/test/main.c index a5ecd5293f..bb27782d88 100644 --- a/src/gr_ore_poly/test/main.c +++ b/src/gr_ore_poly/test/main.c @@ -13,6 +13,7 @@ #include "t-ring.c" #include "t-set_str.c" +#include "t-sigma_delta.c" /* Array of test functions ***************************************************/ @@ -20,6 +21,7 @@ test_struct tests[] = { TEST_FUNCTION(gr_ore_poly_ring), TEST_FUNCTION(gr_ore_poly_set_str), + TEST_FUNCTION(gr_ore_poly_sigma_delta), }; /* main function *************************************************************/ diff --git a/src/gr_ore_poly/test/t-ring.c b/src/gr_ore_poly/test/t-ring.c index 4f5666507f..82768e04f9 100644 --- a/src/gr_ore_poly/test/t-ring.c +++ b/src/gr_ore_poly/test/t-ring.c @@ -27,9 +27,7 @@ TEST_FUNCTION_START(gr_ore_poly_ring, state) gr_ore_poly_ctx_t ore_ctx; slong reps; - gr_ctx_init_random(ctx, state); - /* TODO: ensure ctx has a generator */ - gr_ore_poly_ctx_init_rand(ore_ctx, state, ctx); + gr_ore_poly_ctx_init_randtest2(ctx, ore_ctx, state); if (gr_ctx_is_finite(ctx) == T_TRUE || gr_ctx_has_real_prec(ctx) == T_TRUE) diff --git a/src/gr_ore_poly/test/t-set_str.c b/src/gr_ore_poly/test/t-set_str.c index d8bb3b080f..4a6418b550 100644 --- a/src/gr_ore_poly/test/t-set_str.c +++ b/src/gr_ore_poly/test/t-set_str.c @@ -10,8 +10,6 @@ */ #include "test_helpers.h" -#include "fmpz.h" -#include "gr_vec.h" #include "gr_ore_poly.h" /* This is mostly a copy of src/gr_poly/test/t-set_str.c */ @@ -29,9 +27,7 @@ TEST_FUNCTION_START(gr_ore_poly_set_str, state) slong glen; char * s; - gr_ctx_init_random(ctx, state); - /* TODO: ensure ctx has a generator */ - gr_ore_poly_ctx_init_rand(ore_ctx, state, ctx); + gr_ore_poly_ctx_init_randtest2(ctx, ore_ctx, state); gr_ore_poly_init(f, ore_ctx); gr_ore_poly_init(g, ore_ctx); diff --git a/src/gr_ore_poly/test/t-sigma_delta.c b/src/gr_ore_poly/test/t-sigma_delta.c new file mode 100644 index 0000000000..7ce7d343e9 --- /dev/null +++ b/src/gr_ore_poly/test/t-sigma_delta.c @@ -0,0 +1,156 @@ +/* + Copyright (C) 2025 Marc Mezzarobba + + This file is part of FLINT. + + FLINT is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 3 of the License, or + (at your option) any later version. See . +*/ + +#include "test_helpers.h" +#include "gr_ore_poly.h" + +FLINT_DLL extern gr_static_method_table _ca_methods; + +TEST_FUNCTION_START(gr_ore_poly_sigma_delta, state) +{ + for (slong iter = 0; iter < 1000 * flint_test_multiplier(); iter++) + { + gr_ctx_t cctx, ctx; + + gr_ore_poly_ctx_init_randtest2(cctx, ctx, state); + + gr_ptr a = gr_heap_init(cctx); + gr_ptr b = gr_heap_init(cctx); + gr_ptr c = gr_heap_init(cctx); + gr_ptr sa = gr_heap_init(cctx); + gr_ptr sb = gr_heap_init(cctx); + gr_ptr sc = gr_heap_init(cctx); + gr_ptr sc1 = gr_heap_init(cctx); + gr_ptr da = gr_heap_init(cctx); + gr_ptr db = gr_heap_init(cctx); + gr_ptr dc = gr_heap_init(cctx); + gr_ptr dc1 = gr_heap_init(cctx); + + int status = GR_SUCCESS; + + status |= gr_randtest(a, state, cctx); + status |= gr_randtest(b, state, cctx); + + int call = n_randint(state, 4); + switch (call) + { + case 0: + status |= gr_ore_poly_sigma_delta(sa, da, a, ctx); + break; + case 1: + status |= gr_ore_poly_sigma(sa, a, ctx); + status |= gr_ore_poly_delta(da, a, ctx); + break; + case 2: + status |= gr_set(sa, a, cctx); + status |= gr_ore_poly_sigma_delta(sa, da, sa, ctx); + break; + case 3: + status |= gr_set(da, a, cctx); + status |= gr_ore_poly_sigma_delta(sa, da, da, ctx); + break; + default: + FLINT_ASSERT(0); + } + + status |= gr_ore_poly_sigma_delta(sb, db, b, ctx); + + if (status != GR_SUCCESS) + { + if (cctx->which_ring == GR_CTX_GR_POLY) + { + switch (POLYNOMIAL_CTX(cctx)->base_ring->which_ring) + { + case GR_CTX_FMPZ: + case GR_CTX_FMPQ: + break; + default: + goto cleanup; + } + switch (GR_ORE_POLY_CTX(ctx)->which_algebra) + { + case ORE_ALGEBRA_COMMUTATIVE: + case ORE_ALGEBRA_DERIVATIVE: + case ORE_ALGEBRA_EULER_DERIVATIVE: + case ORE_ALGEBRA_FORWARD_SHIFT: + case ORE_ALGEBRA_FORWARD_DIFFERENCE: + case ORE_ALGEBRA_BACKWARD_SHIFT: + case ORE_ALGEBRA_BACKWARD_DIFFERENCE: + case ORE_ALGEBRA_Q_SHIFT: + case ORE_ALGEBRA_MAHLER: + flint_printf("FAIL: unexpected failure\n"); + flint_abort(); + default: + goto cleanup; + } + } + goto cleanup; + } + + status = gr_add(c, a, b, cctx); + status |= gr_ore_poly_sigma_delta(sc, dc, c, ctx); + + if (status == GR_SUCCESS) + { + status = gr_add(sc1, sa, sb, cctx); + if (status == GR_SUCCESS && gr_equal(sc, sc1, cctx) == T_FALSE) + { + flint_printf("FAIL: σ(a + b) = σ(a) + σ(b)\n"); + flint_abort(); + } + + status = gr_add(dc1, da, db, cctx); + if (status == GR_SUCCESS && gr_equal(dc, dc1, cctx) == T_FALSE) + { + flint_printf("FAIL: δ(a + b) = δ(a) + δ(b)\n"); + flint_abort(); + } + } + + status = gr_mul(c, a, b, cctx); + status |= gr_ore_poly_sigma_delta(sc, dc, c, ctx); + + if (status == GR_SUCCESS) + { + status |= gr_mul(sc1, sa, sb, cctx); + if (status == GR_SUCCESS && gr_equal(sc, sc1, cctx) == T_FALSE) + { + flint_printf("FAIL: σ(a·b) = σ(a)·σ(b)\n"); + flint_abort(); + } + + status |= gr_mul(dc1, da, b, cctx); + status |= gr_addmul(dc1, sa, db, cctx); + if (status == GR_SUCCESS && gr_equal(dc, dc1, cctx) == T_FALSE) + { + flint_printf("FAIL: δ(a·b) = δ(a)·b + σ(a)·δ(b)\n"); + flint_abort(); + } + } + +cleanup: + gr_heap_clear(dc1, cctx); + gr_heap_clear(dc, cctx); + gr_heap_clear(db, cctx); + gr_heap_clear(da, cctx); + gr_heap_clear(sc1, cctx); + gr_heap_clear(sc, cctx); + gr_heap_clear(sb, cctx); + gr_heap_clear(sa, cctx); + gr_heap_clear(c, cctx); + gr_heap_clear(b, cctx); + gr_heap_clear(a, cctx); + gr_ore_poly_ctx_clear(ctx); + gr_ctx_clear(cctx); + } + + TEST_FUNCTION_END(state); +} diff --git a/src/gr_poly/compose.c b/src/gr_poly/compose.c index 111d57f17c..60cc300942 100644 --- a/src/gr_poly/compose.c +++ b/src/gr_poly/compose.c @@ -14,16 +14,31 @@ #include "gr_vec.h" #include "gr_poly.h" +int +_gr_poly_inflate(gr_ptr poly, slong len, slong n, gr_ctx_t ctx) +{ + slong i, sz = ctx->sizeof_elem; + int status = GR_SUCCESS; + + for (i = len - 1; i >= 1 && n > 1; i--) + { + gr_swap(GR_ENTRY(poly, i * n, sz), GR_ENTRY(poly, i, sz), ctx); + status |= _gr_vec_zero(GR_ENTRY(poly, (i - 1) * n + 1, sz), n - 1, ctx); + } + + return status; +} + /* compose by poly2 = a*x^n + c, no aliasing; n >= 1 */ int _gr_poly_compose_axnc(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr c, gr_srcptr a, slong n, gr_ctx_t ctx) { slong i, sz = ctx->sizeof_elem; - int status; + int status = GR_SUCCESS; /* shift by c (c = 0 case will be fast) */ - status = _gr_poly_taylor_shift(res, poly1, len1, c, ctx); + status |= _gr_poly_taylor_shift(res, poly1, len1, c, ctx); /* multiply by powers of a */ if (gr_is_one(a, ctx) != T_TRUE) @@ -56,11 +71,7 @@ _gr_poly_compose_axnc(gr_ptr res, gr_srcptr poly1, slong len1, } /* stretch */ - for (i = len1 - 1; i >= 1 && n > 1; i--) - { - gr_swap(GR_ENTRY(res, i * n, sz), GR_ENTRY(res, i, sz), ctx); - status |= _gr_vec_zero(GR_ENTRY(res, (i - 1) * n + 1, sz), n - 1, ctx); - } + status |= _gr_poly_inflate(res, len1, n, ctx); return status; }