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 mpoly docs, functions for discriminant, resultant, term_content, and deflation, and run mod_mpoly doctests. #216

Merged
merged 28 commits into from
Sep 9, 2024

Conversation

Jake-Moss
Copy link
Contributor

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

A small collection of changes I've just pushed. This PR:

  • Adds mpoly docs, addresses Interface for making an mpoly context #195 (comment)
  • Adds mod_mpolys to the doc test run list. Previously they were not being run, I've fixed that up and fixed the docs tests themselves. There doesn't appear to be any other modules missing from the list.
    Ideally these would be discovered and run as part of pytest. There's a pytest plugin for Cython modules, however it relies on in-tree builds. Cython conveniently provides the __test__ dictionary in modules that contain function names + line no. to doc string for doc strings that contain >>>. I think is is all that is required so a pytest plugin could work with meson-python's out of tree builds.
  • Adds functions for discriminant, resultant, term_content, and deflation for all mpoly types.
  • Adds a primitive_part function for fmpz_mpoly.

I've added those functions as I intended to use them with SymPy's DMP_Flint, however these are the univariate functions adapted for multivariate polynomial types, not the multivariate variants that SymPy seems to expect. Regardless they're here now.

@oscarbenjamin
Copy link
Collaborator

  • Previously they were not being run, I've fixed that up and fixed the docs tests themselves. There doesn't appear to be any other modules missing from the list.

We might want to have a test for this a bit like test_all_tests that checks the tests are all listed in test_all.py.

  • Ideally these would be discovered and run as part of pytest. There's a pytest plugin for Cython modules, however it relies on in-tree builds. Cython conveniently provides the __test__ dictionary in modules that contain function names + line no. to doc string for doc strings that contain >>>. I think is is all that is required so a pytest plugin could work with meson-python's out of tree builds.

I'm not quite sure I understand what you mean by this. Firstly I believe that the __test__ thing is what causes doctest failures to be reported more than once sometimes. It doesn't happen consistently though so not sure how to reproduce right now...

Is there a simple way to adapt the current code so that instead of using a hard-coded list of modules it can walk the tree looking for __test__?

Secondly you seem to be saying that there is a plugin but it doesn't work. Or are you saying that we should make a new plugin or send a fix to the existing one?

I agree that it would be better if pytest could do this and if we could arrange it so that spin test would run the doctests by default. Writing a plugin is hard work though. I'm trying to upstream the coverage one.

It would still be nice even if pytest is used if there is a way to run all the doctests like python -m flint.test after install. I like the fact that users can run the tests:

pip install python-flint
python -m flint.test

@oscarbenjamin
Copy link
Collaborator

  • Adds functions for discriminant, resultant, term_content, and deflation for all mpoly types.

This is nice. I wonder for deflation if there should also be deflate and inflate. At least inflate is needed to be able to reverse deflation although I guess compose can do it:

In [6]: from flint import *

In [7]: x = fmpz_poly([0, 1])

In [8]: p = x**4 + x**2 + 1

In [9]: p.deflation()
Out[9]: (x^2 + x + 1, 2)

In [10]: p2, n = p.deflation()

In [11]: p2(x**n)
Out[11]: x^4 + x^2 + 1

Using compose is more awkward for mpoly though.

Looking at this across other types it seems to be somewhat unfortunately designed that we have:

  • p.inflate(n) - replace x -> x^n.
  • p.deflate(n) - replace x^n -> x.
  • p.deflation() - find n and replace x^n -> x.

So there is no method to get n without computing the deflated polynomial. The reason we might want n without the polynomial is that you might have more than one polynomial that you want to transform analogously e.g.:

n1 = p1.get_deflation()
n2 = p2.get_deflation()
n = gcd(n1, n2)
p1p = p1.deflate(n)
p2p = p2.deflate(n)

@oscarbenjamin
Copy link
Collaborator

I've added those functions as I intended to use them with SymPy's DMP_Flint, however these are the univariate functions adapted for multivariate polynomial types, not the multivariate variants that SymPy seems to expect.

I'm not sure I understand here. Note that there are many different function in SymPy as well as two different implementations of multivariate polynomials and the fact that any kind of polynomial can have multivariate polynomial coefficients...

We have e.g.

In [2]: Poly(2*x**2*y + 4*y).primitive()
Out[2]: (2, Poly(x**2*y + 2*y, x, y, domain='ZZ'))

In [3]: Poly(2*x**2*y + 4*y, x).primitive()
Out[3]: (2y, Poly(x**2 + 2, x, domain='ZZ[y]'))

In the second case the coefficients would by python-flint types if PolyElement could wrap python-flint's mpoly types.

In the first case the Poly would internally hold a DMP_FLint which would internally hold an fmpz_mpoly. The primitive_part method here only computes the primitive part but not also the content. Of course you can get the content by dividing the polynomials but it would have been better to compute the content and divide the polynomial by the (scalar) content which is what fmpz_poly_primitive_part does except that it discards the content.

Suppose that we had a way to get the coefficients of an fmpz_mpoly as an fmpz_vec and also that fmpz_vec had a gcd method. Then in python-flint (or in SymPy) you could do:

def content(p):
    content = p.coeffs_vec().gcd()

def content_primitive(p):
    primitive_part = p / p.content() # scalar div
    return content, primitive

Apart from making the leading coefficient positive that is exactly what fmpz_poly_primitive_part does. Otherwise we need:

def content_primitive(p):
    primitive_part = p.primitive_part()
    content = p / primitive_part # poly div
    content = content.leading_coefficient()
    if content < 0:
        content = -content
        primitive_part = -primitive_part
    return content, primitive

That's fine but needing a polynomial division is unfortunate.

Note that SymPy would attach the sign to the polynomial rather than the content (not sure if that's actually a good idea but that is what it does):

In [5]: primitive(-x)
Out[5]: (1, -x)

@oscarbenjamin
Copy link
Collaborator

For resultant and discriminant we have these in SymPy:

In [6]: Poly(x + y).resultant(Poly(x - y))
Out[6]: Poly(-2*y, y, domain='ZZ')

In this case we need to go from a ring with 2 generators to a ring with 1 generator. That is not what Flint's resultant etc functions do so we would need some way of converting the fmpz_mpoly from R[x,y] -> R[y].

@oscarbenjamin
Copy link
Collaborator

I think that the changes here look good and could be merged if you are done with this.

@oscarbenjamin
Copy link
Collaborator

Actually I don't think that deflation works quite right because of the shifts:

In [30]: f = x**3 * y + x * y**4 + x * y

In [31]: f
Out[31]: x^3*y + x*y^4 + x*y

In [32]: f_deflated, strides = f.deflation()

In [33]: f_deflated
Out[33]: x + y + 1

In [34]: strides
Out[34]: fmpz_vec(['2', '3'], 2)

In [36]: f_deflated.compose(x**2, y**3)
Out[36]: x^2 + y^3 + 1

In [37]: f_deflated.compose(x**2, y**3) * x*y
Out[37]: x^3*y + x*y^4 + x*y

The output is affects by the shifts but they aren't returned. Either this should not use the shifts or it needs to return them.

@Jake-Moss
Copy link
Contributor Author

I'm not quite sure I understand what you mean by this. Firstly I believe that the __test__ thing is what causes doctest failures to be reported more than once sometimes. It doesn't happen consistently though so not sure how to reproduce right now...

Secondly you seem to be saying that there is a plugin but it doesn't work. Or are you saying that we should make a new plugin or send a fix to the existing one?

I agree that it would be better if pytest could do this and if we could arrange it so that spin test would run the doctests by default.

Sorry I should have clearer, there is a pytest plugin by the name of pytest-cython that provides doctest support for Cython code. I briefly tried and I couldn't get it to work with an out-of-tree build, the documentation says

It is assumed that the C extension modules have been build in place before running py.test and there is a matching Cython .pyx file

I've briefly looked over the code and it's a small single file plugin (<200LOC), from there I saw mention of the __test__ attribute with this comment

       # If somehow the equivalent test does not already exist, we
       # keep the __test__ test (maybe it is something else not
       # generated by autotestdict)

https://github.com/lgpage/pytest-cython/blob/bdd8266f3a9f58e7efd931ab8dbec42bc9e481f7/src/pytest_cython/plugin.py#L137-L140

It seems like the plugin could pick up the "missing" doctests from the modules alone but they are being skipped because

# only run test if matching .so and .pyx files exist

https://github.com/lgpage/pytest-cython/blob/bdd8266f3a9f58e7efd931ab8dbec42bc9e481f7/src/pytest_cython/plugin.py#L45-L46

To me this sounds like the existing plugin could be modified to support running doctests from the modules alone, like how the existing flint.test does.

Writing a plugin is hard work though. I'm trying to upstream the coverage one.

It certainly seems like it, thank you for that. I gave the pytest docs for it a skim but had to move on to other things. I might give patching pytest-cython a crack this weekend but it's low on my list of things to do.

Is there a simple way to adapt the current code so that instead of using a hard-coded list of modules it can walk the tree looking for __test__?

I haven't experimented with this but I believe so.

It would still be nice even if pytest is used if there is a way to run all the doctests like python -m flint.test after install. I like the fact that users can run the tests:

That is nice just as a sanity check.

@Jake-Moss
Copy link
Contributor Author

I'm not sure I understand here. Note that there are many different function in SymPy as well as two different implementations of multivariate polynomials and the fact that any kind of polynomial can have multivariate polynomial coefficients...

Yeah I'm certainly realising the complexity of it. I'm first allowing a DMP_Flint to hold a fmpz_mpoly for the ZZ and QQ domains

For resultant and discriminant we have these in SymPy:

In [6]: Poly(x + y).resultant(Poly(x - y))
Out[6]: Poly(-2*y, y, domain='ZZ')

In this case we need to go from a ring with 2 generators to a ring with 1 generator. That is not what Flint's resultant etc functions do so we would need some way of converting the fmpz_mpoly from R[x,y] -> R[y].

These are the two functions I was referring to, in SymPy DMP.resultants and DMP.subresultants seem to operate without respect to a particular generator, while PolyElement.subresultants and mpoly.resultants do.

Suppose that we had a way to get the coefficients of an fmpz_mpoly as an fmpz_vec and also that fmpz_vec had a gcd method. Then in python-flint (or in SymPy) you could do:

def content(p):
    content = p.coeffs_vec().gcd()

Is this not what mpoly.term_content computes?

For resultant and discriminant we have these in SymPy:

In [6]: Poly(x + y).resultant(Poly(x - y))
Out[6]: Poly(-2*y, y, domain='ZZ')

In this case we need to go from a ring with 2 generators to a ring with 1 generator. That is not what Flint's resultant etc functions do so we would need some way of converting the fmpz_mpoly from R[x,y] -> R[y].

In [1]: from flint import *

In [5]: ctx = fmpz_mpoly_ctx.get_context(2, Ordering.lex, 'x')

In [6]: x, y = ctx.gens()\

In [7]: (x + y).resultant(x - y, 'x0')
Out[7]: -2*x1

In [8]: (x + y).resultant(x - y, 'x1')
Out[8]: 2*x0

This seems to align with what flint computes with respect to the first generator. Not sure if that is a coincidence. I've been having a little difficulty deciphering the DMP docs and implementations.

@oscarbenjamin
Copy link
Collaborator

I don't know if you've read this:
https://docs.sympy.org/latest/modules/polys/domainsintro.html#dmp-representation

The DMP representation is a recursive representation of multivariate polynomials where a polynomial in say QQ[x,y,z] is actually represented like an element of QQ[y,z][x] or more explicitly (QQ[y,z])[x]. It is a univariate polynomial whose coefficients are polynomials in QQ[y,z]. Likewise the elements of QQ[y,z] are represented as elements of QQ[z][y] and so on.

What this means is that a DMP always has a leading variable so if you want to compute the resultant of two multivariate polynomials then it is always known which variable you want to use: the leading one or "main variable". Technically the Sylvester resultant is a concept that is defined for univariate polynomials and the resultant itself should be an element of the ring of coefficients e.g. resultant of two polynomials in D[x] is an element of D not of D[x]. It just happens to be most useful though to compute resultants of multivariate polynomials that are viewed as univariate by choosing a main variable:
https://en.wikipedia.org/wiki/Resultant

The DMP recursive representation makes it natural that an operation like resultant produces a polynomial in a different ring e.g. projecting from QQ[x,y,z] to QQ[y,z] (or from QQ[y,z][x] to QQ[y,z]):

In [11]: p1 = Poly(x + y)

In [12]: p2 = Poly(x - y)

In [13]: p3 = p1.resultant(p2)

In [14]: p1
Out[14]: Poly(x + y, x, y, domain='ZZ')

In [15]: p2
Out[15]: Poly(x - y, x, y, domain='ZZ')

In [16]: p3
Out[16]: Poly(-2*y, y, domain='ZZ')

Notice that p3 does not have x as a generator. The internal representation is different:

In [19]: p1.rep.to_list()
Out[19]: [[1], [1, 0]]

In [20]: p2.rep.to_list()
Out[20]: [[1], [-1, 0]]

In [21]: p3.rep.to_list()
Out[21]: [-2, 0]

This is different from what FLINT's resultant function does because while it returns a polynomial free of x it is still represented as an element of the ring ZZ[x,y,z] i.e. it still has the same context. To implement what SymPy's resultant does there needs to be a way to convert a polynomial that is free of x from ZZ[x,y,z] to ZZ[y,z].

DMP itself does not have any notion of what its generators are:

In [22]: p1.rep
Out[22]: DMP_Python([[1], [1, 0]], ZZ)

The idea that the variables are "x" and "y" is maintained at a higher-level in the Poly type. What this means is that we need a general way to compute resultant of a polynomial wrt main variable and project down to the representation with the main variable removed from the context. We don't care what the variable names are and we always know which one is the main variable.

Note that I am not saying that what you have as a resultant function here is wrong for python-flint but just that it needs to be adapted somehow in SymPy and possibly something extra needs to be provided in python-flint to make that possible.

What I'm not sure about is how to hold the contexts for this sort of thing in memory. You could imagine a stack of contexts where the context for 3 variables holds a reference to the context for 2 variables etc. (I mean that this stack would be in SymPy, not in python-flint, although it could potentially make sense in python-flint as well).

@oscarbenjamin
Copy link
Collaborator

Suppose that we had a way to get the coefficients of an fmpz_mpoly as an fmpz_vec and also that fmpz_vec had a gcd method. Then in python-flint (or in SymPy) you could do:

def content(p):
    content = p.coeffs_vec().gcd()

Is this not what mpoly.term_content computes?

Almost but not quite:

In [26]: (4*x + 2*x*y).term_content()
Out[26]: 2*x

It is the gcd of the terms rather than of the coefficients and so it also includes the gcd of the monomials. SymPy would give:

In [31]: primitive(4*x + 2*x*y)
Out[31]: (2, x*y + 2*x)

In [32]: Poly(4*x + 2*x*y).primitive()
Out[32]: (2, Poly(x*y + 2*x, x, y, domain='ZZ'))

I suppose in python-flint it can be:

In [34]: (4*x + 2*x*y).term_content().leading_coefficient()
Out[34]: 2

It is not ideal to compute the gcd of the monomials when it is not needed though. Note that computing the gcd of many things like gcd(a, b, c, d, ...) is often not as costly as you might imagine because usually the gcd turns out to be 1 fairly early and you can shortcut checking all of the arguments. If you take a random polynomial and multiply it by x though then it will be very likely that the gcd of the coefficients is 1 but the gcd of the terms will not be and so it will not be possible to shortcut in term_content but it would be in some analogous function for computing the gcd of the coefficients.

Internally FLINT computes the equivalent of p.coeffs_vec().gcd():
https://github.com/flintlib/flint/blob/331171713cfb0c1b3e9dc4dbabb9849cbb3f7eb7/src/fmpz_mpoly/primitive_part.c#L33-L35
That operation is not exposed through a public function but we could make one in python-flint that does the same since it is straight-forward. Ideally we don't even need to make an fmpz_vec explicitly.

There are a few different natural notions of what "content" means when we are working with polynomials. For univariate polynomials there are two cases:

  • We might mean gcd of coefficients (an element of ZZ)
  • We might mean gcd of terms (a 1-term element of ZZ[x,y,z]).

Usually when someone says "content of a polynomial" they mean the gcd of coefficients:
https://en.wikipedia.org/wiki/Primitive_part_and_content

In a multivariate context we also need to define our interpretation of what is a coefficient i.e. whether some generators are viewed as being part of the coefficients or part of the monomials for the purpose of defining the content.

  • We might choose a main variable and regard all other variables as being part of the coefficients rather than the terms e.g. we think of ZZ[x,y,z] as ZZ[y,z][x] and so the content is a polynomial in ZZ[y,z]
  • We can also choose some subset of the variables to be generators and regard the remainder as being part of the coefficients e.g. we think of ZZ[x,y,z,t] as ZZ[z,t][x,y] and then the content is a polynomial in ZZ[z,t].

The Flint functions related to this are:

  • fmpz_mpoly_term_content. This treats all variables as generators but computes gcd of terms rather than coefficients. This is not usually what is meant by "content" and there is no analogue of this in SymPy.
  • fmpz_mpoly_content_vars. This allows choosing a subset of variables to be generators and computing the content as a polynomial in the remaining generators. There is no direct analogue of this operation in SymPy because it would be done by converting to a representation that already distinguishes which variables are generators.
  • fmpz_mpoly_primitive_part. This treats all variables as generators and uses the usual definition of content as just the gcd of the coefficients. This is directly analogous to what SymPy's content and primitive functions/methods do except that it does not return the content and also it handles sign normalisation differently.

There are several ways to adapt these functions so that they do what SymPy's equivalent functions would do:

  • c = p.term_content().leading_coefficient() gives the usual content and then prim = p / c gives the primitive part. This is less efficient than calculating c directly.
  • c = p.content_vars([0,1,2,...,n-1]) also gives c but I think it is computed in a less efficient way by a more general algorithm that is designed to do something more complicated.
  • prim = p.primitive_part() gives the usual primitive part and then c = p / prim gives the content but using a polynomial division here is less efficient.

None of these is quite optimal for just getting the normal idea of the "content" which is just the gcd of the integer coefficients. Generally term_content is probably more useful in many situations but still there are situations where we just want to factor out the gcd of the coefficients as an integer. I think that we just need to add a function that does this in python-flint using _fmpz_vec_content much like FLINT's primitive_part does.

What I'm not sure about is what names to give to methods that do these things. Probably there should be a method that is just called content for getting the gcd of the coefficients although we should check what existing names are used for any related poly or mpoly methods. In SymPy the analogous methods/functions are called:

  • c = content(p) get the gcd of coefficients as an integer.
  • (c, prim) = primitive(p) returns both content and primitive part such that p == c*prim.

It makes sense to have these two functions because it isn't possible to compute the primitive part without first computing the content so you might as well have a function called primitive that returns both.

Note also that we should not necessarily always follow exactly Flint's conventions for names of functions and methods etc. Flint itself is not even fully consistent across all types as in they may or may not have exactly the same methods which may or may not be defined to mean exactly the same things. Rather we should define methods that make sense consistently within python-flint from the perspective of someone using python-flint from Python. In some cases that means we need to piece together a few Flint functions to do something useful rather then just exposing precisely Flint's documented functions for each type.

We should also consider that in general it is useful if all types have the same methods e.g. if all polynomial types have primitive and content methods or even if those can be consistent across both univariate and multivariate polynomials. If we are going to have term_content for multivariate polynomials then we probably want it for univariate polynomials as well. It is easy to compute for univariate polynomials so I guess FLINT does not provide a function for it but it's nicer for users if the interface is as uniform as possible. Likewise if we have content and primitive for multivariate polynomials then we also want it for univariate polynomials and it should be defined consistently in both cases.

Note that in SymPy content and primitive are defined not just for polynomials with integer coefficients:

In [5]: primitive(2*x/3 + 2)
Out[5]: (2/3, x + 3)

In [6]: primitive(2*x + 3, modulus=5)
Out[6]: (1, 2x - 2)

The definition is potentially ambiguous over a field but it can still be useful to define some notion of "factor out a scalar".

@oscarbenjamin
Copy link
Collaborator

Note also that we should not necessarily always follow exactly Flint's conventions for names of functions and methods etc.

To summarise:

  • term_content is a useful function with a reasonable name. It makes sense for all polynomial types so let's add that.
  • We also want content though which gives the actual coefficient content and this should be for both fmpz_poly and fmpz_mpoly at least.
  • primitive_part seems poorly designed to me. It should return the content instead of discarding it. I don't think we should copy that interface. It would be better to have a function primitive or content_primitive or something that returns both the content and the primitive part. We could also have a primitive_part method but it does not add much value compared to content_primitive.
  • deflation is good but FLINT has defined it somewhat inconsistently for univariate and multivariate polynomials. I don't think that we should copy that inconsistency. We could have two versions like deflation and shift_deflation or something and then deflation for multivariate polynomials should be analogous to what it does with univariate polynomials.

@oscarbenjamin
Copy link
Collaborator

  • We could have two versions like deflation and shift_deflation or something

Note that this is also sort of unnecessary because you can use term_content for mostly the same effect:

In [1]: from flint import *

In [2]: ctx = fmpz_mpoly_ctx(2, Ordering.lex, ['x', 'y'])

In [3]: x, y = ctx.gens()

In [4]: f = x**3 * y + x * y**4 + x * y

In [5]: f
Out[5]: x^3*y + x*y^4 + x*y

In [6]: tc = f.term_content()

In [7]: tc
Out[7]: x*y

In [8]: fd, s = (f / tc).deflation()

In [9]: fd
Out[9]: x + y + 1

In [10]: s
Out[10]: fmpz_vec(['2', '3'], 2)

In [11]: tc * fd.compose(x**2, y**3)
Out[11]: x^3*y + x*y^4 + x*y

In [12]: tc * fd.compose(x**2, y**3) == f
Out[12]: True

Using term_content like this is not exactly equivalent because it also computes the gcd of the coefficients but in may cases that would be desired anyway. There could be a function (not sure of a good name) to compute the deflation like this:

m, fd, s = f.shift_deflation()
assert m*fd.inflate(s) == f

Another point I realise now is that it would be fine to have __call__ do the equivalent of what compose does. When we discussed this before the __call__ function was doing lots of different things but I think now that subs is a separate function there is probably no reason that you shouldn't be able to do e.g.:

In [13]: fd.compose(x**2, y**3)
Out[13]: x^2 + y^3 + 1

In [14]: fd(x**2, y**3) # maybe this could work
---------------------------------------------------------------------------
TypeError

That does work for univariate polynomials already:

In [15]: x = fmpz_poly([0, 1])

In [16]: p = x**2 + 1

In [17]: p(p)
Out[17]: x^4 + 2*x^2 + 2

@Jake-Moss
Copy link
Contributor Author

deflation is good but FLINT has defined it somewhat inconsistently for univariate and multivariate polynomials. I don't think that we should copy that inconsistency. We could have two versions like deflation and shift_deflation or something and then deflation for multivariate polynomials should be analogous to what it does with univariate polynomials.

Regarding this, what should these functions be? If I'm understanding correctly to make the interface consistent with the univariate polynomials the shift vector should optionally be zeros. So the deflation function could optionally accept the two vectors, stride and shift as arguments then if provided, use them, otherwise call fmpz_mpoly_deflation to compute the stride and shift, then call fmpz_mpoly_deflate with just the stride, ignoring the shift vector?
To make the above interface useful we'd then also need to have another method to compute the stride and shift so they can be passed as arguments. On top of deflation we could then provide an inflation method that takes the stride as a required argument and optionally the shift to undo the deflation.

@Jake-Moss
Copy link
Contributor Author

The DMP recursive representation makes it natural that an operation like resultant produces a polynomial in a different ring e.g. projecting from QQ[x,y,z] to QQ[y,z] (or from QQ[y,z][x] to QQ[y,z]):

That makes sense thank you.

This is different from what FLINT's resultant function does because while it returns a polynomial free of x it is still represented as an element of the ring ZZ[x,y,z] i.e. it still has the same context. To implement what SymPy's resultant does there needs to be a way to convert a polynomial that is free of x from ZZ[x,y,z] to ZZ[y,z].

There should still be a context coercion function laying around somewhere. Here it is

    # TODO: Rethink context conversions, particularly the proposed methods in #132
    def coerce_to_context(self, ctx):
        cdef:
            fmpz_mod_mpoly res
            slong *C
            slong i

        if not typecheck(ctx, fmpz_mod_mpoly_ctx):
            raise ValueError("provided context is not a fmpz_mod_mpoly_ctx")

        if self.ctx is ctx:
            return self

        C = <slong *> libc.stdlib.malloc(self.ctx.val.minfo.nvars * sizeof(slong *))
        if C is NULL:
            raise MemoryError("malloc returned a null pointer")
        res = create_fmpz_mod_mpoly(self.ctx)

        vars = {x: i for i, x in enumerate(ctx.py_names)}
        for i, var in enumerate(self.ctx.py_names):
            C[i] = <slong>vars[var]

        fmpz_mod_mpoly_compose_fmpz_mod_mpoly_gen(res.val, self.val, C, self.ctx.val, (<fmpz_mod_mpoly_ctx>ctx).val)

        libc.stdlib.free(C)
        return res

It works by composing the mpolys with the generators as the functions to convert between contexts. It's rather rough in this form but could be made to work. I'm not sure if Flint provides a better method than this.

What I'm not sure about is how to hold the contexts for this sort of thing in memory. You could imagine a stack of contexts where the context for 3 variables holds a reference to the context for 2 variables etc. (I mean that this stack would be in SymPy, not in python-flint, although it could potentially make sense in python-flint as well).

I think context interop and coercion is something is missing from python-flint (and Flint). A tree of contexts and their variables could work but it would have to be lazy or risk blowing up for lots of variables. Either way to change domain we need a new context, I think creating one on the fly and coercing within SymPy is reasonable for now, just storing them in python-flints context cache. It would have to be an explicit operation that's only used either when required or user-side to prevent slow down.

@Jake-Moss
Copy link
Contributor Author

Thank you for all your comments here, I'll try to get to the rest tomorrow evening.

@oscarbenjamin
Copy link
Collaborator

A tree of contexts and their variables could work but it would have to be lazy or risk blowing up for lots of variables

At least for DMP we don't need a tree. Remember that DMP does not care what the names of the variables are and only ever wants to compute the resultant wrt to the main variable. So it would be a singly linked list of contexts for different numbers of variables like the 100 variable context holds a reference to the 99 variable context and so on. I suppose that a DMP_Flint will need to have a reference to a *_mpoly_ctx. We could wrap the _mpoly_ctx in a Python class that could lazily generate and then cache a context for one fewer variable like:

class Flint_mpoly_ctx:
    def __init__(self, ...):
        self.fctx = fmpz_mpoly_ctx(...)
        self.fctx_drop = None

   @property
   def drop_ctx(self):
       if self.fctx_drop is None:
            ...
      return self.fctx_drop

Alternatively we could also just @lru_cache the function that creates one of these contexts for DMP.

I think context interop and coercion is something is missing from python-flint (and Flint).

It does look like fmpz_mod_mpoly_compose_fmpz_mod_mpoly_gen is the most direct function for this though. Internally it looks like it makes a matrix to perform the transformation but I would have expected some more direct way.

I think this is possibly something that the univar structs can be used for although I'm not sure:
https://flintlib.org/doc/fmpz_mpoly.html#univariate-functions

@oscarbenjamin
Copy link
Collaborator

So the deflation function could optionally accept the two vectors, stride and shift as arguments then if provided, use them, otherwise call fmpz_mpoly_deflation to compute the stride and shift, then call fmpz_mpoly_deflate with just the stride, ignoring the shift vector?

We already have the methods inflate, deflate and deflation for different univariate polynomial types although not all types have them. Obviously the multivariate ones need to return a list of integers for the deflation whereas the univariate ones return an int so there has to be that difference. Otherwise though we should try to keep them as consistent as possible.

I suggest the following:

class foo_poly:

    def inflate(self, n: int) -> foo_poly:
        """q  s.t.  q(x) = p(x^n)"""

    def deflate(self, n: int) -> foo_poly:
        """q  s.t.  p(x) = q(x^n)"""

    def deflation(self) -> (foo_poly, int):
        """(q, n)  s.t.  p(x) = q(x^n)  for maximal n"""

    def deflation_monom(self) -> (foo_poly, int, foo_poly):
        """(q, n, m)  s.t.  p(x) = m * q(x^n)  for maximal n and monomial m"""

    def deflation_index(self) -> (foo_poly, int, int):
        """(n, i)  s.t.  p(x) = x^i * q(x^n)  for maximal n and i"""


class foo_mpoly:

    def inflate(self, N: list[int]) -> foo_mpoly:
        """q  s.t.  q(X) = p(X^N)"""

    def deflate(self, N: list[int]) -> foo_mpoly:
        """q  s.t.  p(X) = q(X^N)"""

    def deflation(self) -> (foo_mpoly, list[int]):
        """(q, N)  s.t.  p(X) = q(X^N)  for maximal N"""

    def deflation_monom(self) -> (foo_mpoly, list[int], foo_mpoly):
        """(q, N, m)  s.t.  p(X) = m * q(X^N)  for maximal N and monomial m"""

    def deflation_index(self) -> (foo_mpoly, list[int], list[int]):
        """(N, I)  s.t.  p(X) = X^I * q(X^N)  for maximal N and I"""

Then the invariants for both cases are like:

assert p.inflate(n).deflate(n) == p

q, n = p.deflation()
assert q.inflate(n) == p

q, n, m = p.deflation_monom()
assert m * q.inflate(n) == p

n, i = p.deflation_index()
q = (p/x**i).deflate(n)
assert q.inflate(n) * x**i == p

The last case for deflation_index is a bit awkward with mpoly because you would need to somehow turn the list of indices i into a monomial but that is fine. Possibly we can do without deflation_index altogether or at least for now it is not necessary.

This keeps consistency with the different types in python-flint. The equivalent methods in SymPy correspond to inflate and deflation (called deflate in SymPy).

@oscarbenjamin
Copy link
Collaborator

It wouldn't be too bad to require pytest unconditionally. We could also add it to the requirements as an extra like pip install python-flint[test].

One thing though is that when I am debugging a crash (core dump etc) that happens during the test suite the output from pytest is particularly unhelpful. Typically I have to add print(1), print(2) etc throughout the code and then run with python -m flint.test so that I can actually see print output. Maybe there is a way to make that better with pytest.

@oscarbenjamin
Copy link
Collaborator

It wouldn't be too bad to require pytest unconditionally.

Although we need to ensure that it works after install with

pytest --pyargs flint

and also we should ensure that the currently documented approach still works:

python -m flint.test

@oscarbenjamin
Copy link
Collaborator

Only downside is there are a lot of dots, I'm sure we can live with that though.

I think that the dots are a good thing. You can also pass -v to turn each into something that prints the filename. It is hardly needed for python-flint but with pytest-xdist the different tests can be run in parallel with pytest -nauto.

@Jake-Moss
Copy link
Contributor Author

I think that the dots are a good thing. You can also pass -v to turn each into something that prints the filename. It is hardly needed for python-flint but with pytest-xdist the different tests can be run in parallel with pytest -nauto.

Previously they were identifiable in pytests verbose mode, I've just fixed that

@Jake-Moss
Copy link
Contributor Author

Only downside is there are a lot of dots, I'm sure we can live with that though.

I think that the dots are a good thing. You can also pass -v to turn each into something that prints the filename. It is hardly needed for python-flint but with pytest-xdist the different tests can be run in parallel with pytest -nauto.

These both work for me (when run under spin)

@oscarbenjamin
Copy link
Collaborator

I left a few comments. The only critical ones are about adding hand-written definitions in flintlib/functions. That can't work after gh-217. In gh-219 I have added a README in flintlib that explains this.

Otherwise this looks good. I think that the inflate etc functions are doing the right thing now and it is good that they are consistent across types.

@Jake-Moss
Copy link
Contributor Author

I believe all those should be resolved now. I've also added a skip directive to the arb.hypgeom_2f1 doctest, this brings my test runtime down from ~21s to ~2.5. Happy for that commit to be dropped if you'd prefer to have that run every time

Comment on lines +734 to +743
def content(self):
"""
Return the GCD of the coefficients of ``self``.

>>> fmpz_poly([3, 6, 0]).content()
3
"""
cdef fmpz res = fmpz()
_fmpz_vec_content(res.val, self.val.coeffs, self.val.length)
return res
Copy link
Collaborator

Choose a reason for hiding this comment

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

We may need to think a bit about the sign convention for primitive and content. Typically the primitive part is defined to be a canonical representative in some sense which may mean that the content should absorb the sign although there are different ways of defining this.

@oscarbenjamin
Copy link
Collaborator

Okay, this looks good.

Also the test runner changes are nice.

Let's get this in.

Thanks!

@oscarbenjamin oscarbenjamin merged commit 6cdc1bc into flintlib:main Sep 9, 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