Fmpz mod mpoly and nmod mpoly

[Jake-Moss]

This PR aims to add support for the fmpz_mod_mpoly and nmod_mpoly flint objects. I'm using the recently added fmpz_mpoly as a basis, there may be remants of the large copy-paste that have escaped my changes.

I haven't given the docs any more attention that a quick search-replace as of yet so while they (and all tests locally) may be passing I'll have to give them a once over before I start on nmod_mpoly.

[Jake-Moss]

For now I've crammed the fmpz_mod_mpoly into the mpoly tests with both a prime and non-prime modulus. Though it's not quite as clean as Id like. The division functions assume the modulus is prime, so I've set it to raise an exception on any division function with a non-prime modulus, even if it would be within the domain. The same goes for sqrt and gcd. While the flint docs don't state they require a prime modulus for gcd, it was raising an impossible inversion error so I added the exception to it as well.

[oscarbenjamin]

I've set it to raise an exception on any division function with a non-prime modulus, even if it be within the domain. The same goes for sqrt and gcd.

That seems right to me.

Since nmod_mpoly has a context (unlike other nmod_* types; see gh-124) we could also check whether the modulus is prime when the context is created.

[Jake-Moss]

Since nmod_mpoly has a context (unlike other nmod_* types; see gh-124) we could also check whether the modulus is prime when the context is created.

I think that's also reasonable, currently I'm just checking the first time, storing the result and just returning that instead for later calls.

[oscarbenjamin]

I think that's also reasonable, currently I'm just checking the first time, storing the result and just returning that instead for later calls.

Either approach is reasonable I think. Not doing it when the context is created saves the one time cost of the check if it is never needed. In exchange we get a slight bit of overhead each time the attribute needs to be checked. I don't know whether that overhead is even measurable though.

What would be measurable I think is actually computing the prime check every time division happens. Particularly in the case of nmod rather than nmod_mpoly I would expect that to bring a measurable slowdown. This is why we need a context for nmod so that we can store whether or not the modulus is prime and only do the prime check once.

[Jake-Moss]

This is why we need a context for nmod so that we can store whether or not the modulus is prime and only do the prime check once.

That sounds reasonable, is that something we'd like in this PR?

I've brought the nmod_mpoly up to the same spec as fmpz_mod_mpoly. I've refrained from cramming them into the other nmod_poly and fmpz_mod_poly tests respectively just because of large API differences. Happy to have a crack at it but would either require a little more than just adding a for loop + some lambdas.

Whoops I looks like I missed that numpy isn't required here, was just using it for 32bit/64bit dependent memoryviews in place of a malloc of ulongs.

[oscarbenjamin]

This is why we need a context for nmod so that we can store whether or not the modulus is prime and only do the prime check once.

That sounds reasonable, is that something we'd like in this PR?

That's up to you. If it makes things easier in this PR then go ahead. Otherwise we can leave that to a separate PR.

More generally I think that we want all types to have associated contexts so that there is a uniform interface like ZZ.poly([1,2]) or GF(3).poly([1,2]) and we can express operations like ZZ.poly([1,2]).set_domain(QQ). That means we want contexts for fmpz, fmpq, arb, ...

We would also need contexts for all types for generic rings as well. I haven't looked enough into the generic rings interfaces to see if that should affect the overall design of contexts for the non-generic types.

For arb we need a contexts that can hold non-global precision. Likewise for the number of terms in a series. It isn't possible for SymPy, mpmath etc to make use of things like arf if the precision is global (e.g. SymPy and mpmath would both clobber each other). Even methods like __repr__ depend on global variables so there needs to be some sort of printer context or something.

If you want to add contexts for nmod then I would suggest just keeping it simple so that it resembles fmpz_mod. I mention these other things here just to give a sense for where I would want to go in future.

[oscarbenjamin]

The test failure is

Traceback (most recent call last):
  File "/opt/hostedtoolcache/Python/3.12.4/x64/lib/python3.12/site-packages/flint/test/__main__.py", line 36, in run_tests
    test()
  File "/opt/hostedtoolcache/Python/3.12.4/x64/lib/python3.12/site-packages/flint/test/test_all.py", line 2888, in test_mpolys
    assert P(ctx=ctx).is_zero() is True
           ^^^^^^^^^^^^^^^^^^
AttributeError: 'flint.types.fmpz_mod_mpoly.fmpz_mod_mpoly' object has no attribute 'is_zero'

I think that is a test added in gh-168 so you may need to merge master to see the same failure locally.

[Jake-Moss]

I think that is a test added in gh-168 so you may need to merge master to see the same failure locally.

Yep that was it thanks! Seems like CI was doing something and merging the PR before testing, definitely through me off for a minute.

If you want to add contexts for nmod then I would suggest just keeping it simple so that it resembles fmpz_mod. I mention these other things here just to give a sense for where I would want to go in future.

I think something like that might be best suited in it's own PR, this one only touches on nmod_mpoly and fmpz_mod_mpoly

[Jake-Moss]

More generally I think that we want all types to have associated contexts so that there is a uniform interface like ZZ.poly([1,2]) or GF(3).poly([1,2]) and we can express operations like ZZ.poly([1,2]).set_domain(QQ). That means we want contexts for fmpz, fmpq, arb, ...

This seems like it could be done as system that builds upon the existing contexts. It would have to be Python types rather than Cython types as cdef classs don't support multiple inheritance if something OO driven is desired.

We would also need contexts for all types for generic rings as well. I haven't looked enough into the generic rings interfaces to see if that should affect the overall design of contexts for the non-generic types.

Just from a quick skim they seem like quite the beast.

For arb we need a contexts that can hold non-global precision. Likewise for the number of terms in a series. It isn't possible for SymPy, mpmath etc to make use of things like arf if the precision is global (e.g. SymPy and mpmath would both clobber each other). Even methods like __repr__ depend on global variables so there needs to be some sort of printer context or something.

Would this be replacing the global FlintContext? Instead moving it towards something instance specific? Or at least swappable?

[oscarbenjamin]

Seems like CI was doing something and merging the PR before testing, definitely through me off for a minute.

That is usually how CI works at least for Travis and GitHub Actions. We were using Cirrus as well before which did not seem to work like that.

When you push to a PR (or close and open it) the CI system will make a temporary branch somewhere effectively like:

git checkout master
git checkout -b tmp_branch
git merge pr_branch # merge the PR into temporary copy of master
<run CI checks>

Sometimes you need to close and open a PR to get it to merge with latest master for CI (asking GA to "rerun" a job does not trigger a new merge).

There are also a bunch of complicated access rules in GitHub so if the PR modifies the GA workflow files then it might still run with the master version of the workflows depending on the event type that triggers the workflow.

[Jake-Moss]

Looks like another possible memory error. I'll hunt that down tomorrow.

[oscarbenjamin]

Would this be replacing the global FlintContext? Instead moving it towards something instance specific?

These are details that need to be worked out. For arf for example I would like to use it for SymPy's RR domain: sympy/sympy#26861. It also makes sense for mpmath to use arf at least as its internal representation for mpf. That would mean that both SymPy and mpmath start manipulating the global precision and would mess each other up. Likewise the global precision is not thread-safe even if SymPy is the only user of it.

The thread-safe part could be managed by using a thread-local variable in python-flint. I don't know what the overheads would be for that in the context of arf though. I still think SymPy and mpmath would clobber each other though because it is likely that SymPy wants to modify the precision exactly when it is doing arithmetic with mpf or calling mpmath functions. Using a thread-local but otherwise global precision would still require very careful usage of context managers like with ctx.extraprec() everywhere and that just seems fragile.

The cleanest solution is to use separate contexts with context-local precision etc settings. Then you can do e.g.:

ctx = arb_ctx(53)
result = ctx.sin(1)

Then SymPy and mpmath can each make their own contexts and manipulate the precision separately. Of course if SymPy makes this context a global variable and modifies its precision then we're back to not thread-safe etc.

In the case of something like the RR domain though we can attach a separate context to each domain object and we never need to modify its precision. Ideally creating this arb_ctx would be a very cheap operation so that it can be say the first step in a call to something like .evalf(). Then everything is naturally thread-safe and robust since we've just avoided having any global variables to worry about.

The problem with having separate arb contexts everywhere is that then you need to do e.g. ctx.add(x, y) rather than x + y. For many purposes that is fine but it's not really nice for end users and wouldn't work for SymPy's RR domain where the elements need to work in generic code that uses elementary arithmetic. Now then what you need is to have each arb, arf etc hold a reference to the arb_ctx. Then in arithmetic like x + y the x.__add__ method can check the context for the other object and if they have the same context then its precision can be used for the operation and the returned result can also hold a reference to the same context.

Then we could ensure that any arb created as e.g. flint.arb(2) just ends up using the global context. That way users who just want to set the precision globally can do so but it doesn't mess up SymPy or mpmath who would never use the global context directly.

Now the problem is what are you going to do when you have x + y and they have different contexts with different precision or even different rounding modes? Let's see what mpmath does:

In [1]: import mpmath

In [2]: ctx1 = mpmath.MPContext()

In [3]: ctx2 = mpmath.MPContext()

In [4]: ctx2.dps = 50

In [5]: ctx1.mpf(1) / ctx2.mpf(3)
Out[5]: mpf('0.33333333333333331')

In [6]: ctx2.mpf(1) / ctx1.mpf(3)
Out[6]: mpf('0.33333333333333333333333333333333333333333333333333311')

It basically just takes the context from the left-hand operand. I think that it would be better to say something like:

• If there are mixed context operands then take the one with higher precision (since dropping precision is irreversible)?
• Definitely error if they have different rounding modes though.
• Maybe actually just make all mixed context operations an error and require an explicit conversion or the use of context methods like x1+ctx1(x2) or ctx1.add(x1,x2).

There are likely two very different groups of users here:

• End users who are probably happy using a global context. These users probably want implicit conversions.
• Library users like SymPy, mpmath etc who are probably happy using explicit conversions. From SymPy's perspective it is nicer for debugging if mixed context operations give an immediate error so that we can fix it and add the explicit conversion in the right place.

Perhaps a way to satisfy both groups is to say that implicit coercions between contexts are allowed only when mixing the global context with a non-global context (the global context can coerce all others). Then SymPy can never use the global context but end users who do use it can take the output of SymPy functions and mix it into their own types.

Similar considerations apply to all of the other operations like the number of terms in a series how exactly everything gets printed. The problem with printing is you might have something like:

result = str([obj1, obj2])

Now the list's __str__ method is going to call obj1.__repr__(). Adding obj1.repr() doesn't help us here. This is basically the same problem as with x + y. One of two things is needed:

• Either we have a printer object that takes care of printing lists and things as well so you do printer.str([obj1, obj2]). Now __repr__ and __str__ are not used at all.
• Or otherwise each object needs to somehow hold a reference to a context that holds its printing settings but that adds to the memory footprint of every object.

Again the ideal solution here is probably different for different groups of users. SymPy for example already has its own complicated printing machinery so being able to call printer.str(...) or obj.repr() works fine. Many end users will find that inconvenient though and may want something more like the current global variables.

[oscarbenjamin]

There should be some handling equality checking against scalars:

In [26]: flint.fmpz_poly([1]) == 1
Out[26]: True

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

In [28]: ctx.from_dict({(0,0):1})
Out[28]: 1

In [29]: ctx.from_dict({(0,0):1}) == 1
Out[29]: False

[oscarbenjamin]

I've asked how to handle this on the Cython mailing list:
https://groups.google.com/g/cython-users/c/sw1D4DdJI0U

The suggested solution from Cython maintainer is to do:

cdef extern from *:
    """
    #define SIZEOF_ULONG sizeof(ulong)
    """
    int SIZEOF_ULONG

That is marginally better than using FLINT_BITS // 8.

It is also apparently considered to be a Cython bug that the code for sizeof(ulong) is not generated correctly. Apparently it should be okay to have an extern typedef like:

cdef extern from "flint/flint.h":
   ctypedef unsigned long ulong
   ctypedef long slong

Cython should be expected to respect the signed-ness of the types but otherwise not make any assumptions about the actual size or whether it really is typedef'd to unsigned long in the C code.

It doesn't look like we actually use sizeof(ulong) or sizeof(slong) anywhere else so hopefully there are no other problems lurking as a result of this.

[oscarbenjamin]

Tests are failing but only under Python 3.12. Specifically this test:

            assert p.subs({"x1": S(0), 0: S(0)}) == ctx.from_dict({(0, 0): 1})
>           assert p.compose(p.subs({"x1": 0}), ctx.from_dict({(0, 1): 1})) == mpoly({
                (2, 2): 36,
                (1, 2): 24,
                (1, 0): 9,
                (0, 2): 4,
                (0, 1): 2,
                (0, 0): 4
            })
E           assert 64*x0^4*x1^2 ...^2 + 2*x1 + 10 == 36*x0^2*x1^2 ...1^2 + 2*x1 + 4

The output is definitely quite different from that expected in the test:

(Pdb) type(p)
<class 'flint.types.fmpz_mod_mpoly.fmpz_mod_mpoly'>
(Pdb) p p.compose(p.subs({"x1": 0}), ctx.from_dict({(0, 1): 1}))
64*x0^4*x1^2 + 96*x0^3*x1^2 + 31*x0^2*x1^2 + 12*x0^2 + 72*x0*x1^2 + 9*x0 + 36*x1^2 + 2*x1 + 10

It looks like subs is the problem:

(Pdb) p p
4*x0^2*x1^2 + 3*x0 + 2*x1 + 1
(Pdb) p p.subs({"x1":0})
4*x0^2 + 3*x0 + 3

If we set x1 -> 0 then we should just get 3*x0 + 1 here. It seems to be setting x1 -> 1 instead.

The fmpz_mod_mpoly.subs method is here:

python-flint/src/flint/types/fmpz_mod_mpoly.pyx

Lines 894 to 925 in 043ca9d

	def subs(self, dict_args) -> fmpz_mod_mpoly:
	"""
	Partial evaluate this polynomial with select constants. Keys must be generator names or generator indices,
	all values must be fmpz.
	>>> from flint import Ordering
	>>> ctx = fmpz_mod_mpoly_ctx.get_context(2, Ordering.lex, 11, 'x')
	>>> f = ctx.from_dict({(0, 0): 1, (1, 0): 2, (0, 1): 3, (1, 1): 4})
	>>> f.subs({"x1": 0})
	2*x0 + 1
	"""
	cdef:
	fmpz_mod_mpoly res
	slong i
	args = tuple((self.ctx.variable_to_index(k), v) for k, v in dict_args.items())
	# Partial application with args in Z. We evaluate the polynomial one variable at a time
	res = create_fmpz_mod_mpoly(self.ctx)
	fmpz_mod_mpoly_set(res.val, self.val, self.ctx.val)
	for i, arg in args:
	if typecheck(arg, fmpz_mod):
	fmpz_mod_mpoly_evaluate_one_fmpz(res.val, res.val, i, (<fmpz_mod>arg).val, self.ctx.val)
	else:
	arg_fmpz = any_as_fmpz(arg)
	if arg_fmpz is NotImplemented:
	raise TypeError(f"cannot coerce {type(arg)} to fmpz")
	fmpz_mod_mpoly_evaluate_one_fmpz(res.val, res.val, i, (<fmpz>arg).val, self.ctx.val)
	return res

[oscarbenjamin]

If factor would raise an exception for non-prime modulus then that should be tested here so that the test shows the behaviour in all cases.

Also can these not be written with g in the output?

        characteristic_zero = not (P is flint.fmpz_mod_mpoly or P is flint.nmod_mpoly)

        if characteristic_zero:
            # Primitive factors over Z for fmpz_mpoly and fmpq_mpoly
            assert (f * g).factor() == (S(2), [(g / 2, 1), (f, 1)])
        elif is_field:
            # Monic polynomials over Z/pZ for nmod_mpoly and fmpz_mod_mpoly
            assert (f * g).factor() == (S(8), [(g / 2, 1), (f / 4, 1)])
        else:
            # Factorisation not allowed over Z/nZ for n not prime.
            # Flint would abort so we raise an exception instead:
            assert raises(...)

One thing that I hadn't noticed until now is that fmpq_mpoly.factor behaves differently compared to fmpq_poly.factor. With fmpq_poly factors are always monic:

In [29]: f = fmpq_poly([1, 4])

In [30]: g = fmpq_poly([2, 2])

In [31]: f
Out[31]: 4*x + 1

In [32]: g
Out[32]: 2*x + 2

In [33]: (f*g).factor()
Out[33]: (8, [(x + 1/4, 1), (x + 1, 1)])

However fmpq_mpoly gives factors primitive over Z:

In [35]: ctx = fmpq_mpoly_ctx(1, Ordering.lex, ["x"])

In [36]: f = fmpq_mpoly('4*x + 1', ctx)

In [37]: g = fmpq_mpoly('2*x + 2', ctx)

In [38]: (f*g).factor()
Out[38]: (2, [(x + 1, 1), (4*x + 1, 1)])

I'm wondering whether it makes sense to try to make this consistent somehow. One sort of consistency is that whenever we have a field all factors can be monic which would suggest that fmpq_mpoly should also return monic factors.

On the other hand though it is computationally better for Q to use primitive factors and absorb any denominator into the constant factor as fmpq_mpoly does so we could change fmpq_poly to return primitive factors instead:

python-flint/src/flint/types/fmpq_poly.pyx

Lines 396 to 416 in ba1a7fe

	def factor(self):
	"""
	Factors *self* into irreducible polynomials. Returns (*c*, *factors*)
	where *c* is the leading coefficient and *factors* is a list of
	(*poly*, *exp*) pairs with all *poly* monic.
	>>> fmpq_poly.legendre_p(5).factor()
	(63/8, [(x, 1), (x^4 + (-10/9)*x^2 + 5/21, 1)])
	>>> (fmpq_poly([1,-1],10) ** 5 * fmpq_poly([1,2,3],7)).factor()
	(-3/700000, [(x^2 + 2/3*x + 1/3, 1), (x + (-1), 5)])
	"""
	c, fac = self.numer().factor()
	c = fmpq(c)
	for i in range(len(fac)):
	base, exp = fac[i]
	lead = base[base.degree()]
	base = fmpq_poly(base, lead)
	c *= lead ** exp
	fac[i] = (base, exp)
	return c / self.denom(), fac

In SymPy it is fmpq_poly that is inconsistent with SymPy's behaviour and needs to be special-cased to remove precisely the denominators that python-flint inserts:
https://github.com/sympy/sympy/blob/31e957f47a636df9b7e62cd3951ec0afcc266369/sympy/polys/polyclasses.py#L2285-L2304

[Jake-Moss]

I'm wondering whether it makes sense to try to make this consistent somehow. One sort of consistency is that whenever we have a field all factors can be monic which would suggest that fmpq_mpoly should also return monic factors.

On the other hand though it is computationally better for Q to use primitive factors and absorb any denominator into the constant factor as fmpq_mpoly does so we could change fmpq_poly to return primitive factors instead:

In my mind I would expect python-flint to be largely transparent with the results it returns from flint, if flint is returning primitive factors then python-flint should as well, if not to just to be consistent. I've seen your #189 PR and will take a look now

[oscarbenjamin]

I had a couple of comments about the tests but generally I think this looks good and could be merged.

Is there anything outstanding to do?

[Jake-Moss]

Is there anything outstanding to do?

I don't believe so, nothing comes to mind. I will however open another PR to remove the if isinstance ... elif typecheck ladders from fmpz_mpoly and fmpq_mpoly as well once this is merged.

[oscarbenjamin]

Okay, this looks good. Let's get it in!

[oscarbenjamin]

This is nice because it is mote of the mpoly types from Flint now.

The only ones missing I think are gr_mpoly (all of gr is missing) and fq_nmod_mpoly (which makes less sense since we use fq_default rather than fq_nmod).

Certainly these are all of the mpoly types that SymPy would currently use until we add gr.