Begin adding fmpz_mod_poly #87

GiacomoPope · 2023-09-20T15:04:39Z

Opening a PR just to have the ability for group discussion while working on this new type. Only work so far is an import of the level0 files, and a simple print test to ensure everything compiles happily.

GiacomoPope · 2023-09-20T15:27:20Z

Main thing I'm not so sure about is the new context extension, whether this should be built on top of the fmpz_mod_ctx or instead be a class which takes fmpz_mod_ctx as an argument instead. This second option would allow the creation of fmpz_mod elements by passing the context stored in fmpz_mod_poly_ctx as an arg...

Something like

cdef class fmpx_mod_poly_ctx:
    def __init__(self, mod_ctx):
        self.val = mod_ctx

Then the way of working would be:

F = fmpz_mod_ctx(163)
a = F(3) # fmpz_mod(3, 163)

R = fmpz_mod_poly_ctx(F)
f = R([1,2,3]) # fmpz_mod_poly([1,2,3], 11)

Of course, we could allow either, so we could have:

F = fmpz_mod_ctx(163)
R = fmpz_mod_poly_ctx(F)
assert R == fmpz_mod_poly_ctx(163)

Where when 163 can be cast to fmpz type, we first create the fmpz_mod_ctx and then use this as the value for fmpz_mod_poly_ctx.

oscarbenjamin · 2023-09-20T16:18:59Z

If we were working at a higher level I might suggest something like R = ZZ_n(163).poly_ring() or something.

Otherwise I prefer this version:

F = fmpz_mod_ctx(163)
R = fmpz_mod_poly_ctx(F)

The reason I like it is because then R can keep a reference to F and then any method from fmpz_mod_poly that returns fmpz_mod can always ensure to use the exact same context object F. That way any calculation that mixes R and F can also work and it can be fast to check that an element of R matches with an element of F etc in R(2)*F([1,2]).

If we had cached contexts then it wouldn't matter because fmpz_mod_poly_ctx(163) could retrieve F from the cache.

src/flint/types/fmpz_mod_poly.pyx

GiacomoPope · 2023-09-20T16:47:43Z

Made the change above about the __getitem__() and also reworked the context. I have also changed __hash__ for mod_ctx and mod_poly_ctx as before it was only the modulus being hashed and I figured these two objects should hash to different values.

oscarbenjamin · 2023-09-20T23:44:49Z

I figured these two objects should hash to different values.

It is okay for them to hash to the same value. The hash invariant requires only that equal objects have the same hash. It is not required that unequal objects have different hash so if you want all hashes can just be 1. In practice I doubt that it matters much as I can't imagine someone putting these two objects into a dict/set together at least not in a performance sensitive context.

GiacomoPope · 2023-09-21T14:23:26Z

There's a whole bunch of functions available in flint which might not make it into the first PR for this type, but I might add something like

def berlekamp_massey(self):
    raise NotImplementedError("Function available in Flint but not yet implemeneted in python-flint")

To entice others to work on it if they're interested in having the feature.

GiacomoPope · 2023-09-21T16:19:16Z

So I think I'm missing something about how factor works for Z/NZ.

In flint, I call:

>>> from flint import *
>>> F = fmpz_mod_ctx(163)
>>> R = fmpz_mod_poly_ctx(F)
>>> f = R([1,2,3])
>>> f.factor()
[(x + 115, 1), (x + 103, 1)]

sage: F = GF(163)
sage: R.<x> = PolynomialRing(F)
sage: f = R([1,2,3])
sage: f.factor()
(3) * (x + 103) * (x + 115)

We see that the constant term is missing.

Now, in the fmpz_poly_factor_struct, there's this fmpz_t value which is called c and is the constant term, but for fmpz_mod_poly we have:

    ctypedef struct fmpz_mod_poly_factor_struct:
        fmpz_mod_poly_struct * poly
        slong *exp
        slong num
        slong alloc
    ctypedef fmpz_mod_poly_factor_struct fmpz_mod_poly_factor_t[1]

and there's no storage for this constant term. Is the intention that this is only ever called on monic polynomials?

Another option I guess would be to first collect the leading coefficient with *fmpz_mod_poly_lead to get a pointer to the leading coefficient, then factor the input polynomial (which is treated as a monic polynomial) and then return a tuple (leading_coeff, monic_factors)

GiacomoPope · 2023-09-21T16:47:59Z

Added a few methods to do what I suggest above:

>>> from flint import *
>>> F = fmpz_mod_ctx(163)
>>> R = fmpz_mod_poly_ctx(F)
>>> f = R([1,2,3])
>>> f.leading_coefficient()
fmpz_mod(3, 163)
>>> f.factor()
(fmpz_mod(3, 163), [(x + 115, 1), (x + 103, 1)])
>>> f.monic()
x^2 + 55*x + 109
>>> f.monic().factor()
(fmpz_mod(1, 163), [(x + 115, 1), (x + 103, 1)])

oscarbenjamin · 2023-09-21T17:33:22Z

Is the intention that this is only ever called on monic polynomials?

The doc says:

Factorises a non-constant polynomial f into monic irreducible factors

So I guess it always makes the polynomial monic.

The docs don't say if the modulus is assumed to be prime but I imagine that many of the factor related functions would only work for prime modulus in which case the polynomial could always be made monic.

Looks like the modulus should be prime:

In [2]: import flint

In [3]: F = flint.fmpz_mod_ctx(12)
   ...: R1 = flint.fmpz_mod_poly_ctx(F)

In [4]: p = R1([1, 2, 3, 4])

In [5]: p.factor()
Flint exception (Impossible inverse):
    Exception in fmpz_mod_inv: Cannot invert.
Aborted (core dumped)

I guess that failure is from trying to make it monic.

Here it is monic already but the modulue is not prime:

In [1]: import flint

In [2]: F = flint.fmpz_mod_ctx(12)
   ...: R1 = flint.fmpz_mod_poly_ctx(F)

In [3]: p = R1([4, 3, 2, 1])

In [4]: p.factor()
GR_MUST_SUCCEED failure: src/fmpz_mod_poly/gcd.cAborted (core dumped)

src/flint/types/fmpz_mod.pyx

src/flint/types/fmpz_mod_poly.pyx

GiacomoPope · 2023-09-21T17:45:32Z

What's your opinion on adding a is_prime bool to fmpz_mod_ctx, which is set on initalisation and checked for the cases when the modulus is needed to be prime? We could set this with a check_primality flag, which people can set to off if for some reason they want to include some enormous modulus?

oscarbenjamin · 2023-09-21T18:13:47Z

We probably should have prime checking by default if the alternative is a segfault. If it can be done only once per context then I imagine that is fine.

src/flint/types/fmpz_mod_poly.pyx

GiacomoPope · 2023-09-22T10:38:32Z

Pushed some more work from last night.

fmpz_mod_ctz now checks for primality on initialisation, saving us from segfaults
Many more new methods, with docstrings for most (all)?
Arithmetic is not tested yet, it's boring and I will do it soon when I am feeling more patient

GiacomoPope · 2023-09-22T13:28:28Z

Added support to get random elements from both fmpz_mod_ctx and fmpz_mod_poly_ctx, other methods and addressed the above issue with __cinit__

src/flint/types/fmpz_mod_poly.pyx

oscarbenjamin · 2023-09-25T20:18:11Z

Let me know if at some point you would just like to get what is here merged and move on. Not everything needs to be added in one pull request.

I would like anything that is added to be tested though. We should aim for 100% line coverage in this library because it is mostly wrapper code. The bugs are unlikely to be in the detailed algorithms and are mostly just going to be that calling a function with the wrong arguments causes a segfault or something. Knowing that a function here ever succeeds without a segfault is probably enough to have some reasonable confidence that it works for most generic inputs. It is most important just to see that a very simple test passes on each platform rather than checking particularly complicated or difficult inputs.

GiacomoPope · 2023-09-25T20:46:50Z

I'm going to clean up what I have an add a few extra functions and then start working on coverage. I had my wedding yesterday and so have been not-coding, but will get this ready for merging this week :)

oscarbenjamin · 2023-09-29T17:09:23Z

src/flint/test/test.py

-        assert raises(lambda: 1 / P([1, 1]), TypeError)
-        assert raises(lambda: P([1, 2, 1]) / P([1, 1]), TypeError)
-        assert raises(lambda: P([1, 2, 1]) / P([1, 2]), TypeError)
+        # TODO:
+        # I think this should be a ValueError and not a Type Error?
+        # assert raises(lambda: 1 / P([1, 1]), TypeError)
+        # assert raises(lambda: P([1, 2, 1]) / P([1, 1]), TypeError)
+        # assert raises(lambda: P([1, 2, 1]) / P([1, 2]), TypeError)


The way this is defined for the other poly types is just that it is a TypeError. Note that P([1, 2, 1]) / P([1, 1]) can be an exact division but it gives a TypeError: you need to use // or divmod instead. I think that this was a deliberate choice by @fredrik-johansson so that the operations are always well defined.

I did add scalar division where the ground domain is a field though (nmod_poly and fmpq_poly):

In [14]: flint.fmpq_poly([1, 1]) / 7 Out[14]: 1/7*x + 1/7

In SymPy's algebra system there is another operation which is call exquo that performs an exact division when possible or otherwise raises an ExactQuotientFailed error. The exquo operation is very useful but it is not implemented as a / b in SymPy at least. You have to use it explicitly:

In [16]: ZZ.exquo(ZZ(6), ZZ(3)) Out[16]: mpz(2) In [17]: ZZ.exquo(ZZ(6), ZZ(5)) --------------------------------------------------------------------------- ExactQuotientFailed

Note also that fmpz refuses an exact division:

In [18]: flint.fmpz(5)/5 --------------------------------------------------------------------------- TypeError Traceback (most recent call last) Cell In[18], line 1 ----> 1 flint.fmpz(5)/5 TypeError: unsupported operand type(s) for /: 'flint.types.fmpz.fmpz' and 'int'

Oh ok. Thank you for the detailed reply. Seems like the best thing then is to remove the support for division in this way and potentially have a specialised formula.

For scalar division, this is a TODO atm but could easily be implemented to allow the work to match other polynomials.

I don't know whether we would potentially want to change a / b to be something like exquo. There is a lot of inconsistency in defining a / b for non-fields between different environments even just within the Python ecosystem. Maybe it is better just to avoid that ambiguity by not allowing it at all.

oscarbenjamin · 2023-09-29T17:16:56Z

Currently flint returns None when no sqrt is found, I think this should be a ValueError

I don't have a strong preference here.

There is no integral method for fmpz_mod_poly available from FLINT, so I have to skip this for fmpz_mod_poly

It is easy to add this in Python though:

In [30]: p = flint.fmpq_poly([5, 5, 5, 5, 5])

In [31]: p2 = flint.fmpq_poly([0] + [c/n for n, c in enumerate(p.coeffs(), 1)])

In [32]: p.integral() == p2
Out[32]: True

GiacomoPope · 2023-09-29T17:55:21Z

Ok, made no sqrt be a value error, removed the exact division from _div_ and made it a sep. function, fixed the pow todo and implemented integral pythonically.

Main thing now is get back to 100% coverage and clean up the tests file, and also figure out how to solve the __init__ issue which I have introduced.

Edit:

I also called exact division exact_division, but if you like quoex, we can do this. I don't really mind. More converned about the missing init at the moment :/

oscarbenjamin · 2023-09-29T22:31:19Z

src/flint/types/fmpz_mod_poly.pyx

+        else:
+            return self.evalutate(input)
+
+    def evalutate(self, input):


oscarbenjamin · 2023-09-29T22:37:29Z

More converned about the missing init at the moment :/

That's the thing that concerns me as well. Otherwise I would say let's just merge and continue to improve later but we can't mess around with getting the memory management right: it needs to be correct now and always to the best of our knowledge or otherwise we invite a world of pain in future.

GiacomoPope · 2023-09-29T23:15:58Z

So I've tried something, but I don't know if I like it...

    def __cinit__(self, val, ctx):
        if not typecheck(ctx, fmpz_mod_poly_ctx):
            raise TypeError
        self.ctx = ctx
        fmpz_mod_poly_init(self.val, self.ctx.mod.val)

    def __dealloc__(self):
        fmpz_mod_poly_clear(self.val, self.ctx.mod.val)

    def __init__(self, val, ctx):
        if not typecheck(ctx, fmpz_mod_poly_ctx):
            raise TypeError
        self.ctx = ctx
        check = self.ctx.set_any_as_fmpz_mod_poly(self.val, val)
        if check is NotImplemented:
            raise TypeError 
        self.initialized = True

If we then call

res = fmpz_mod_poly.__new__(fmpz_mod_poly, None, context)

Instead of simply:

res = fmpz_mod_poly.__new__(fmpz_mod_poly)

Then __cinit__() has the ctx it needs to initialise the value. However, I really don't like having to pass None for the value argument when initialising :/

This seems to compile fine though and all the tests pass

oscarbenjamin · 2023-09-30T11:34:27Z

So I've tried something, but I don't know if I like it...

If __cinit__ raises TypeError then __dealloc__ will be called and will try to call fmpz_mod_poly_clear without any previous call to fmpz_mod_poly_init.

I think the way that you had the code before was good but just we need to ensure that an fmpz_mod_poly_init function is always used when __new__ is used. It can be like:

def __cinit__(self):
    self.initialized = False

def __dealloc__(self):
    if self.initialized:
        fmpz_mod_poly_clear(self.val, self.ctx.mod.val)

def __init__(self, val, ctx):
    if not typecheck(ctx, fmpz_mod_poly_ctx):
        raise TypeError
    self.ctx = ctx
    fmpz_mod_poly_init(obj.val, ctx)
    self.initialized = True
    check = self.ctx.set_any_as_fmpz_mod_poly(self.val, val)
    if check is NotImplemented:
        raise TypeError

Then when you use __new__ you would do:

obj = fmpz_mod_poly.__new__(fmpz_mod_poly)
fmpz_mod_poly_init(obj.val, ctx)

In other words the code using this becomes responsible for calling the init function if it chooses to use __new__.

The advantage of letting the caller call init is that they can choose which init function to call e.g. in this case there is also fmpz_mod_poly_init2. Since there are only two init functions for fmpz_mod_poly though I would be inclined just to make static methods that handle calling them as well as __new__ so then the caller does:

obj = fmpz_mod_poly.new(ctx) # already initialised
obj = fmpz_mod_poly.new_alloc(ctx, n) # already initalised

Alternatively we can have __cinit__ set this up like:

def __cinit__(self, val, ctx):
    self.initialized = False
    if not typecheck(ctx, fmpz_mod_poly_ctx):
        raise TypeError
    self.ctx = ctx
    fmpz_mod_poly_init(self.val, ctx)
    self.initialized = True

def __dealloc__(self):
    if self.initialized:
        fmpz_mod_poly_clear(self.val, self.ctx)

Which can be used like your version but it handles the initialized/dealloc correctly if __cinit__ does not initialise self.val.

oscarbenjamin · 2023-09-30T11:42:16Z

A simpler version without the initialised flag could be:

def __cinit__(self, val, ctx):
    # self.ctx is already set to None here I think
    if not typecheck(ctx, fmpz_mod_poly_ctx):
        raise TypeError
    fmpz_mod_poly_init(self.val, ctx)
    self.ctx = ctx

def __dealloc__(self):
    if self.ctx is not None:
        fmpz_mod_poly_init(self.val, self.ctx)

oscarbenjamin · 2023-09-30T12:00:22Z

src/flint/types/fmpz_mod_poly.pyx

+    def __cinit__(self, val, ctx):
+        if not typecheck(ctx, fmpz_mod_poly_ctx):
+            raise TypeError
+        self.ctx = ctx
+        fmpz_mod_poly_init(self.val, self.ctx.mod.val)
+
+    def __dealloc__(self):
+        if self.ctx is not None:
+            fmpz_mod_poly_clear(self.val, self.ctx.mod.val)
+
+    def __init__(self, val, ctx):
+        check = self.ctx.set_any_as_fmpz_mod_poly(self.val, val)
+        if check is NotImplemented:
+            raise TypeError


Okay, I think this is correct now.

oscarbenjamin · 2023-09-30T12:09:34Z

src/flint/types/fmpz_mod_poly.pyx

+
+        G = fmpz_mod_poly.__new__(fmpz_mod_poly, None, self.ctx)
+        S = fmpz_mod_poly.__new__(fmpz_mod_poly, None, self.ctx)
+        T = fmpz_mod_poly.__new__(fmpz_mod_poly, None, self.ctx)


I would make a staticmethod on ctx for this:

@staticmethod cdef fmpz_mod_poly new(fmpz_mod_poly_ctx self): return fmpz_mod_poly.__new__(fmpz_mod_poly, None, self)

Then it is something like:

cdef fmpz_mod_poly ctx cdef fmpz_mod_poly G, S, T ctx = self.ctx G = ctx.new() S = ctx.new() T = ctx.new()

I would check how the generated C code looks to see exactly what is happening though.

Sorry, not a staticmethod, just an ordinary method. I was going to suggest a staticmethod on fmpz_mod_poly like

G = fmpz_mod_poly.new(ctx)

but then I realised that it makes more sense as an instance method on the context:

G = ctx.new()

oscarbenjamin · 2023-09-30T13:32:39Z

I think this PR is probably good to go but I need to find time to read through it in detail.

GiacomoPope · 2023-09-30T14:55:45Z

I think I’m not quite at 100% coverage again so I’ll address this and the above comments as well. Hopefully by this evening

GiacomoPope · 2023-09-30T19:22:21Z

Added a new poly method, which i think makes the code easier to read and recovered 100% coverage.

src/flint/types/fmpz_mod_poly.pyx

oscarbenjamin · 2023-10-01T00:14:20Z

src/flint/types/fmpz_mod_poly.pyx

+            res[i] = (u, exp)
+        return self.leading_coefficient(), res
+
+    def roots(self, multiplicities=True):


I generally prefer separate functions rather than a boolean flag

So I think the default should be with multiplicities, so you think a root_without_multiplicities() function would be better?

Coming up with good names is hard :)

Maybe root_set?

In this situation I always ant to check what other systems do...

Root set is definitely better than my suggestion. Haha. For SageMath, they have a bool flag.

For SageMath, they have a bool flag.

Yeah, so does SymPy:

In [1]: roots((x**2 - 1)**2) Out[1]: {-1: 2, 1: 2} In [2]: roots((x**2 - 1)**2, multiple=True) Out[2]: [-1, -1, 1, 1]

I think it is better though if a function has a clear return type and so if you want to return different types then you should have different function names. Also it is clearer just to write roots_multiple(...) rather than e.g. roots(..., multiple=True).

So if I thought not multiplicities was default then I'd agree but I think it's fine for a user to essentially never know this is an option as worst case they filter out the multiplicities themselves in python? The bool flag here is really only for people who want it but it's not needed.

Agree about consistent return types tho. Not sure what's best

I'm not sure what's best either really.

It is usually better just to have separate functions if there is a need to provide different options. It can also be good just to refuse to provide any options that are potentially redundant.

For python-flint right now I would aim to avoid (postpone?) providing anything that is not strictly necessary because it is a lot of work just to provide the strictly necessary things when you consider that all types (fmpz_poly etc) should ideally support the same usage.

Maybe we could just remove the option to not compute multiplicities then? I'm not at my laptop now but I can edit it out tomorrow.

We can leave it for now.

src/flint/types/fmpz_mod_poly.pyx

oscarbenjamin · 2023-10-01T00:25:47Z

I think this looks good. I'll merge once CI completes.

oscarbenjamin · 2023-10-01T01:09:56Z

Okay looks good. Thanks!

Begin adding fmpz_mod_poly

fea7def

oscarbenjamin reviewed Sep 20, 2023

View reviewed changes

src/flint/types/fmpz_mod_poly.pyx Outdated Show resolved Hide resolved

Modify methods for contexts

6b1611a

Add conversion and comparison functions

6797150

Start work on factor; bug

e6b28d4

Monic, Leading Coefficent and Factor

4e64490

oscarbenjamin reviewed Sep 21, 2023

View reviewed changes

src/flint/types/fmpz_mod.pyx Outdated Show resolved Hide resolved

oscarbenjamin reviewed Sep 21, 2023

View reviewed changes

src/flint/types/fmpz_mod.pyx Outdated Show resolved Hide resolved

oscarbenjamin reviewed Sep 21, 2023

View reviewed changes

src/flint/types/fmpz_mod_poly.pyx Show resolved Hide resolved

More progress, but badly unittested currently

36724e3

oscarbenjamin reviewed Sep 21, 2023

View reviewed changes

src/flint/types/fmpz_mod_poly.pyx Outdated Show resolved Hide resolved

GiacomoPope added 2 commits September 22, 2023 11:36

Add additional bool to check if modulus is prime on init

0a84867

Continue adding functions and docstrings

6122937

Address review, add random elements and further methods and docstrings

4827066

oscarbenjamin reviewed Sep 22, 2023

View reviewed changes

src/flint/types/fmpz_mod_poly.pyx Outdated Show resolved Hide resolved

Add more functions to fmpz_mod_poly, also add is_unit to fmpz_mod

00b87f5

oscarbenjamin reviewed Sep 29, 2023

View reviewed changes

Address TODO in tests and modify how _div_ works

8fd39d1

oscarbenjamin reviewed Sep 29, 2023

View reviewed changes

src/flint/types/fmpz_mod_poly.pyx Outdated

else:

return self.evalutate(input)

def evalutate(self, input):

Copy link

Collaborator

oscarbenjamin Sep 29, 2023

Choose a reason for hiding this comment

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

evaluate

Hacky fix for memory?

3614a42

Simplify init function

9dfd197

better handling of dealloc

1945859

oscarbenjamin reviewed Sep 30, 2023

View reviewed changes

GiacomoPope added 2 commits September 30, 2023 19:35

Add new poly method to context

57b8346

Return to 100% coverage

6a3b447

Forgot to push test commit

52dee64

oscarbenjamin reviewed Sep 30, 2023

View reviewed changes

src/flint/types/fmpz_mod_poly.pyx Outdated Show resolved Hide resolved

oscarbenjamin reviewed Oct 1, 2023

View reviewed changes

src/flint/types/fmpz_mod_poly.pyx Outdated Show resolved Hide resolved

Fix docstring types and initalise int type

b2543b4

oscarbenjamin merged commit 98c2883 into flintlib:master Oct 1, 2023

oscarbenjamin mentioned this pull request Nov 15, 2023

Add fmpz_mod_mat #106

Merged

oscarbenjamin mentioned this pull request Apr 8, 2024

Support for multivariate power series #131

Closed

oscarbenjamin mentioned this pull request Aug 14, 2024

Add fq_default_poly #182

Merged

Begin adding fmpz_mod_poly #87

Begin adding fmpz_mod_poly #87

Conversation

GiacomoPope commented Sep 20, 2023

GiacomoPope commented Sep 20, 2023 • edited Loading

oscarbenjamin commented Sep 20, 2023

GiacomoPope commented Sep 20, 2023

oscarbenjamin commented Sep 20, 2023

GiacomoPope commented Sep 21, 2023

GiacomoPope commented Sep 21, 2023

GiacomoPope commented Sep 21, 2023

oscarbenjamin commented Sep 21, 2023

GiacomoPope commented Sep 21, 2023

oscarbenjamin commented Sep 21, 2023

GiacomoPope commented Sep 22, 2023

GiacomoPope commented Sep 22, 2023

oscarbenjamin commented Sep 25, 2023

GiacomoPope commented Sep 25, 2023

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

oscarbenjamin commented Sep 29, 2023

GiacomoPope commented Sep 29, 2023 • edited Loading

Choose a reason for hiding this comment

oscarbenjamin commented Sep 29, 2023

GiacomoPope commented Sep 29, 2023 • edited Loading

oscarbenjamin commented Sep 30, 2023

oscarbenjamin commented Sep 30, 2023

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

oscarbenjamin commented Sep 30, 2023

GiacomoPope commented Sep 30, 2023 via email • edited Loading

GiacomoPope commented Sep 30, 2023

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

oscarbenjamin commented Oct 1, 2023

oscarbenjamin commented Oct 1, 2023

GiacomoPope commented Sep 20, 2023 •

edited

Loading

GiacomoPope commented Sep 29, 2023 •

edited

Loading

GiacomoPope commented Sep 29, 2023 •

edited

Loading

GiacomoPope commented Sep 30, 2023 via email •

edited

Loading