Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add some mpoly context util functions #228

Merged
merged 16 commits into from
Sep 22, 2024
Merged

Conversation

Jake-Moss
Copy link
Contributor

@Jake-Moss Jake-Moss commented Sep 16, 2024

This adds methods to

  • remove generators from a context, either via index of variable name, and
  • coerce a mpoly from one context to another, either with an inferred mapping or a provided one.

As well as some doc string updates and support for negative indecies in the variable_to_index method.

Missing ATM:

  • ability to add generators, should be easy,
  • test cases, and
  • the fmpz_mod_mpoly_compose_fmpz_mod_mpoly_gen function in Flint, PR here Add missing fmpz_mod_mpoly_compose_ functions flint#2068
    I've disabled support for it ATM and prevented the linker from complaining on import with a macro.

@oscarbenjamin
Copy link
Collaborator

  • coerce a mpoly from one context to another, either with an inferred mapping or a provided one.

We should check what generic rings does before adding any sort of coercion like this. I didn't add coercion in gh-218 but we could easily add some explicit coercion there.

@oscarbenjamin
Copy link
Collaborator

we could easily add some explicit coercion there.

There are two functions related to this gr_set_other and gr_ctx_cmp_coercion.

@Jake-Moss
Copy link
Contributor Author

There are two functions related to this gr_set_other and gr_ctx_cmp_coercion.

I've taken a quick look at these and their mpoly variants.

The gr_ctx_cmp_coercion function says "Returns 1 if coercing elements into ctx1 is more meaningful, and returns -1 otherwise.", it seems to measure this by the ordering of the gr_which_structure enum. So if ctx1->which_ring is defined later than ctx2->which_ring then it is "more meaningful".

The coerce_to_context method in this PR only works between mpoly contexts of the same type. In that cause we'd fall all the way through in this function and return 1.

// flint/src/gr/cmp_coercion.c
int gr_ctx_cmp_coercion(gr_ctx_t ctx1, gr_ctx_t ctx2)
{
    if (ctx1->which_ring < ctx2->which_ring)
        return -1;
    if (ctx1->which_ring > ctx2->which_ring)
        return 1;

    if (ctx1->which_ring == GR_CTX_GR_POLY)
    {
        return gr_ctx_cmp_coercion(POLYNOMIAL_ELEM_CTX(ctx1), POLYNOMIAL_ELEM_CTX(ctx2));
    }

    if (ctx1->which_ring == GR_CTX_GR_MAT)
    {
        return gr_ctx_cmp_coercion(MATRIX_CTX(ctx1)->base_ring, MATRIX_CTX(ctx2)->base_ring);
    }

    return 1;
}

// flint/src/gr.h
typedef enum
{
    GR_CTX_FMPZ, GR_CTX_FMPQ, GR_CTX_FMPZI,
    GR_CTX_FMPZ_MOD, GR_CTX_NMOD, GR_CTX_NMOD8, GR_CTX_NMOD32, GR_CTX_MPN_MOD,
    GR_CTX_FQ, GR_CTX_FQ_NMOD, GR_CTX_FQ_ZECH,
    GR_CTX_NF,
    GR_CTX_REAL_ALGEBRAIC_QQBAR, GR_CTX_COMPLEX_ALGEBRAIC_QQBAR,
    GR_CTX_REAL_ALGEBRAIC_CA, GR_CTX_COMPLEX_ALGEBRAIC_CA,
    GR_CTX_RR_CA, GR_CTX_CC_CA,
    GR_CTX_COMPLEX_EXTENDED_CA,
    GR_CTX_RR_ARB, GR_CTX_CC_ACB,
    GR_CTX_REAL_FLOAT_ARF, GR_CTX_COMPLEX_FLOAT_ACF,
    GR_CTX_NFLOAT, GR_CTX_NFLOAT_COMPLEX,
    GR_CTX_FMPZ_POLY, GR_CTX_FMPQ_POLY, GR_CTX_GR_POLY,
    GR_CTX_FMPZ_MPOLY, GR_CTX_GR_MPOLY,
    GR_CTX_FMPZ_MPOLY_Q,
    GR_CTX_GR_SERIES, GR_CTX_SERIES_MOD_GR_POLY,
    GR_CTX_GR_MAT,
    GR_CTX_GR_VEC,
    GR_CTX_PSL2Z, GR_CTX_DIRICHLET_GROUP, GR_CTX_PERM,
    GR_CTX_FEXPR,
    GR_CTX_UNKNOWN_DOMAIN,
    GR_CTX_WHICH_STRUCTURE_TAB_SIZE
}
gr_which_structure;

Now looking at the gr_set_other method, the fmpz_mpoly gr defines this to be

// flint/src/gr/fmpz_mpoly.c
gr_method_tab_input _gr_fmpz_mpoly_methods_input[] =
{
    // ...
    {GR_METHOD_SET_OTHER,   (gr_funcptr) _gr_fmpz_mpoly_set_other},
    // ...
};

which directly sets the polynomial if the two contexts match in nvars and ordering

int
_gr_fmpz_mpoly_set_other(fmpz_mpoly_t res, gr_srcptr x, gr_ctx_t x_ctx, gr_ctx_t ctx)
{
    if (x_ctx->which_ring == GR_CTX_FMPZ_MPOLY)
    {
        /* fmpz_mpoly_set_fmpz_poly */
        if (MPOLYNOMIAL_MCTX(ctx)->minfo->nvars == MPOLYNOMIAL_MCTX(x_ctx)->minfo->nvars &&
            MPOLYNOMIAL_MCTX(ctx)->minfo->ord == MPOLYNOMIAL_MCTX(x_ctx)->minfo->ord)
        {
            fmpz_mpoly_set(res, x, MPOLYNOMIAL_MCTX(ctx));
            return GR_SUCCESS;
        }
    }

    return gr_generic_set_other(res, x, x_ctx, ctx);
}

In any other case it falls back to

// flint/src/gr_generic/generic.c
int gr_generic_set_other(gr_ptr res, gr_srcptr x, gr_ctx_t xctx, gr_ctx_t ctx)
{
    if (xctx == ctx)
    {
        return gr_set(res, x, ctx);
    }
    else if (xctx->which_ring == GR_CTX_FMPZ)
    {
        return gr_set_fmpz(res, x, ctx);
    }
    else if (xctx->which_ring == GR_CTX_FMPQ)
    {
        return gr_set_fmpq(res, x, ctx);
    }
    else if (xctx->which_ring == GR_CTX_FEXPR)
    {
        gr_vec_t vec;
        fexpr_vec_t fvec;
        int status;
        fexpr_vec_init(fvec, 0);
        gr_vec_init(vec, 0, ctx);
        status = gr_set_fexpr(res, fvec, vec, x, ctx);
        gr_vec_clear(vec, ctx);
        fexpr_vec_clear(fvec);
        return status;
    }
    else
    {
        return GR_UNABLE;
    }
}

which attempts some easy settings and what seems to be a conversion to a symbolic expression, otherwise it returns GR_UNABLE meaning the conversion is unimplemented.

GR_CTX_GR_MPOLY doesn't define an explicit GR_METHOD_SET_OTHER so I believe it falls back to the same generic method above.

@Jake-Moss
Copy link
Contributor Author

It seems like the generic rings are only able to change mpoly context directly when the contexts match exactly. Maybe theres some method to go to a symbolic expression then back to another context, I'm not sure.

The coerce_to_context here is able to convert a mpoly between contexts of the same python-flint type, i.e. fmpz_mpoly_ctx to fmpz_mpoly_ctx via the compose generator functions. It's able to handle a change in nvars by mapping any absent ones to 0, and a change in ordering as the *_mpoly_compose_mat method used internally performs a sort, then combines like terms to normalise the polynomial before returning.

@Jake-Moss Jake-Moss marked this pull request as ready for review September 17, 2024 07:36
@oscarbenjamin
Copy link
Collaborator

Sorry, I missed that this was marked as ready for review.

Comment on lines 944 to 947
def drop_unused_gens(self):
"""
Remove unused generators from this polynomial. Returns a potentially new
context, and potentially new polynomial.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In general I don't like the idea of generating a new context on the fly like this especially based on the value of the polynomial.

I think it should be up to the caller to have the context that they want to use. There could be a context method like ctx.drop_gens() that returns a new context and then the caller should use that for coercing to that specific context.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh, I see there is a ctx.drop_gens method.

What is this particular mpoly.drop_unused_gens method needed for?

I think I would prefer something like:

ctx2 = ctx.drop_gens(p.unused_gens())
p2 = ctx2(p)

In general if you did want to drop some generators you would want to check the generators that are used or not from several polynomials. Perhaps like:

ctx2, [p1_2, p2_2, ...] = ctx.drop_unused_gens([p1, p2, ...])

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is this particular mpoly.drop_unused_gens method needed for?

Just a convenience method from before that skipped the index -> gen name -> index process. unused_gens is definitely more useful and easier enough to compose.

In general if you did want to drop some generators you would want to check the generators that are used or not from several polynomials. Perhaps like:

ctx2, [p1_2, p2_2, ...] = ctx.drop_unused_gens([p1, p2, ...])

Certainly can add this if desired

Comment on lines 968 to 975
def coerce_to_context(self, other_ctx, mapping: dict[str | int, str | int] = None):
"""
Coerce this polynomial to a different context.

This is equivalent to composing this polynomial with the generators of another
context. By default the mapping between contexts is inferred based on the name
of the generators. Generators with names that are not found within the other
context are mapped to 0. The mapping can be explicitly provided.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Mapping the unused generators to zero makes this not really what I would call a "coercion". This is perhaps better described as a "projection".

For a coercion I would expect an error if the element cannot be represented in the target context.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's fair enough. I've renamed it to project_to_context.

For a coercion I would expect an error if the element cannot be represented in the target context.

It could certainly check if that was the case and raise a warning/exception

Comment on lines +509 to +512
>>> mapping # doctest: +SKIP
{3: 1, 4: 0}
>>> list(sorted(mapping.items())) # Set ordering is not stable
[(3, 1), (4, 0)]
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Order of dicts is stable in since Python 3.6.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, however AFAIK the set operations on dictionary view objects are not. I found the order here differed between versions as well. Happy to just have the skipped doc test, the sorting + list is a little clunky for a doc string

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh, I see. Thinking about it though maybe the simple way to represent a mapping from integer to integer is just a list of integers. In this case it could be mapping -> [4, 3].

Comment on lines +976 to +979
>>> f.project_to_context(ctx2)
3*a + 4*b
>>> f.project_to_context(ctx2, mapping={"x": "a", "b": 0})
5*a
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Setting "b": 0 here seems ambiguous to me. Is it setting "b" to the variable of index 0 or to the value 0?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fair. It means the generator at index 0 here. I meant to show that it also accepts generator indices. Perhaps I should make that clearer

@Jake-Moss
Copy link
Contributor Author

Sorry, I missed that this was marked as ready for review.

No worries at all

@oscarbenjamin
Copy link
Collaborator

Okay, let's get this in.

@oscarbenjamin oscarbenjamin merged commit d76d40b into flintlib:main Sep 22, 2024
40 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants