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

Begin adding fmpz_mod_poly #87

Merged
merged 36 commits into from Oct 1, 2023

Conversation

GiacomoPope
Copy link
Contributor

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
Copy link
Contributor Author

GiacomoPope commented Sep 20, 2023

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
Copy link
Collaborator

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.

@GiacomoPope
Copy link
Contributor Author

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
Copy link
Collaborator

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
Copy link
Contributor Author

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
Copy link
Contributor Author

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
Copy link
Contributor Author

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
Copy link
Collaborator

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)

@GiacomoPope
Copy link
Contributor Author

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
Copy link
Collaborator

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.

@GiacomoPope
Copy link
Contributor Author

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
Copy link
Contributor Author

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__

@oscarbenjamin
Copy link
Collaborator

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
Copy link
Contributor Author

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 :)

@GiacomoPope
Copy link
Contributor Author

I have just pushed a commit which attempts to put in the new type into the generic test you wrote.

Some todo:

  • When a true div does not work, you expect a TypeError, but I believe this should be a ValueError
  • Currently flint returns None when no sqrt is found, I think this should be a ValueError
  • There is no integral method for fmpz_mod_poly available from FLINT, so I have to skip this for fmpz_mod_poly

I suppose the thing to then do is remove the specific tests for fmpz_mod_poly where there is duplication to minimise the test file, but I wanted to fix these todo first

Comment on lines 1975 to 2363
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)
Copy link
Collaborator

Choose a reason for hiding this comment

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

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'

Copy link
Contributor Author

Choose a reason for hiding this comment

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

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.

Copy link
Collaborator

Choose a reason for hiding this comment

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

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
Copy link
Collaborator

  • 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
Copy link
Contributor Author

GiacomoPope commented Sep 29, 2023

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 :/

else:
return self.evalutate(input)

def evalutate(self, input):
Copy link
Collaborator

Choose a reason for hiding this comment

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

evaluate

@oscarbenjamin
Copy link
Collaborator

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
Copy link
Contributor Author

GiacomoPope commented Sep 29, 2023

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
Copy link
Collaborator

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
Copy link
Collaborator

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)

Comment on lines +322 to +335
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
Copy link
Collaborator

Choose a reason for hiding this comment

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

Okay, I think this is correct now.


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)
Copy link
Collaborator

Choose a reason for hiding this comment

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

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.

Copy link
Collaborator

Choose a reason for hiding this comment

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

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
Copy link
Collaborator

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

@GiacomoPope
Copy link
Contributor Author

GiacomoPope commented Sep 30, 2023 via email

@GiacomoPope
Copy link
Contributor Author

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

res[i] = (u, exp)
return self.leading_coefficient(), res

def roots(self, multiplicities=True):
Copy link
Collaborator

Choose a reason for hiding this comment

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

I generally prefer separate functions rather than a boolean flag

Copy link
Contributor Author

Choose a reason for hiding this comment

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

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

Copy link
Collaborator

Choose a reason for hiding this comment

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

Coming up with good names is hard :)

Maybe root_set?

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

Copy link
Contributor Author

Choose a reason for hiding this comment

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

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

Copy link
Collaborator

Choose a reason for hiding this comment

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

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).

Copy link
Contributor Author

Choose a reason for hiding this comment

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

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

Copy link
Collaborator

Choose a reason for hiding this comment

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

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.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

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.

Copy link
Collaborator

Choose a reason for hiding this comment

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

We can leave it for now.

@oscarbenjamin
Copy link
Collaborator

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

@oscarbenjamin
Copy link
Collaborator

Okay looks good. Thanks!

@oscarbenjamin oscarbenjamin merged commit 98c2883 into flintlib:master Oct 1, 2023
16 checks passed
@oscarbenjamin oscarbenjamin mentioned this pull request Nov 15, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

2 participants