Modules over Nemo.ZZ fail due to divrem/div definitions

[wbhart]

AbstractAlgebra internally defines div/divrem differently to Julia for Julia integer types (to agree with mod not rem, so that there's a minimal set of representatives mod n).

This means that the AbstractAlgebra tests for modules over ZZ will pass, but when they are run for modules over ZZ using Nemo, they will fail, as Nemo does not internally define div/divrem in the same way AbstractAlgebra does, but continues to use the Julia definitions.

julia> using Nemo

Welcome to Nemo version 0.17.5

Nemo comes with absolutely no warranty whatsoever

julia> import AbstractAlgebra

julia> function rand_module(R::AbstractAlgebra.Ring, vals...)
          rk = rand(0:5)
          M = FreeModule(R, rk)
          levels = rand(0:3)
          for i = 1:levels
             if ngens(M) == 0
                break
             end
             G = [rand(M, vals...) for i in 1:rand(1:ngens(M))]
                   S, f = sub(M, G)
             if rand(1:2) == 1
                M, f = quo(M, S)
             else
                M = S
             end
          end
          return M
       end
rand_module (generic function with 1 method)

julia> for R in [ZZ, QQ]
             for iter = 1:100
                M = rand_module(R, -10:10)

                I, f = snf(M)

                K, g = kernel(f)

                         @test length(invariant_factors(K)) == 0

                m = rand(I, -10:10)

                @test m == inv(f)(f(m))
                      end
          end
Test Failed at REPL[17]:13
  Expression: m == (inv(f))(f(m))
   Evaluated: (-1, -3) == (1, -3)
ERROR: There was an error during testing

The test is taken from AbstractAlgebra's test/generic/InvariantFactorDecomposition-test.jl.

The correct fix for this is not to define div/divrem for fmpz differently to Julia (we already agreed they should concur with Julia's definition) but to do the following:

• do not import div/divrem from Base in Nemo
• export the current definitions of div/divrem for fmpz in Nemo as Base.div and Base.divrem only (so that this is what the user of Nemo will see) (*)
• define div/divrem inside Nemo to agree with the internal AbstractAlgebra definitions for Julia integers (which agree with mod not rem); these must extend AbstractAlgebra.div/divrem
• import Nemo.div/divrem not Base.div/divrem in Oscar (**)
• fix all the bugs in Nemo/Hecke/Oscar that arise as a result of this change
• write up some developer documentation somewhere that explains this (***)

[wbhart]

It goes without saying that this issue could hit anyone at any time (and probably does already without them knowing it). It has nothing to do with modules per se.

[wbhart]

If all goes well time-wise, I intend to close this issue today, unless there is an absolute uproar about what I proposed.

I strongly believe this is an important issue, perhaps the one thing that bugs me the most about Nemo.

I am not intent on (*) if there are objections.

I will not personally do (**) as the owners of the Oscar package will want to do what works for them.

The point (***) would need a central location for developer documentation/decisions.

[We have had many "canonical" places for such decisions to be written up, none of which I can now recall the URL's for. We should decide on a single canonical place (I suggest in the actual documentation of Nemo and/or Oscar, whichever is more appropriate for the particular issue, as this is where other developers will actually look, instead of some random collaborative pad, a wiki, meeting minutes, random text files, on the oscar-dev list, GitHub tickets, a blog, the Oscar website, group emails or worse still in a purely verbal agreement, etc. I'm not personally opposed to a PEP system, so long as it gets written up somewhere that is accessible from the documentation and doesn't make use of a technology here today and gone tomorrow.)]

[thofma]

Does this mean Base.div(::fmpz, ::fmpz) and Nemo.div(::fmpz, ::fmpz) are different? So in Nemo and in the REPL the behaviors are different? This makes working with the only useful integer type that we have more cumbersome, just to cope with some questionable design desicions for BigInt, which no one should use anyway. We should make working with Nemo types more easy, not more complicated.

If we need a different function, just give it a different name. But having this subtle difference between Base.div and Nemo.div is just asking for trouble. I think this will introduce more confusion down the road. It seems that the Module code assumes more than the mathematical definition of division with remainder promises (or what the people introducing Base.div imagined). If one needs canonical representatives, it should be a different function.

To emphasize: Apparently, the people who invented these methods of Base.div for julia integer types thought that this is the interface for integers and this is the expected behaviour. I personally find it a bit confusing if we introduce a second layer of functions with the same names(!) with slightly different behavior and it will forever be hit and miss if you get the correct one.

Edit: If not a different function name, how about a keyword argument or something like that. div(a, b, unique = true) or div(a, b, canonical = true)? Not sure if this is better than a different function name, just bouncing off some ideas here.

[wbhart]

Oh dear. Sorry, but I am not going to go over this for the hundredth time.

div behaves differently in AbstractAlgebra and Base already. div and divrem are made to agree with each other in AbstractAlgebra, and this is required for the generic code to work.

This has nothing to do with Modules. This bug will hit you in many different places.

If we go with your definition, we will have to change EVERY div and divrem in AbstractAlgebra, Nemo, Oscar to div(a, b, true).

The definition must agree with what Flint does for division with remainder, and we cannot rewrite Flint from scratch just because Julia has an inconsistency.

[wbhart]

Why didn't you complain about this before I implemented it? Or on one of the many other tickets where I have painstakingly gone through all our options?

[thofma]

Can you point me to the issue where we agreed on this "fix"?

[thofma]

The main problem is that we try to use julia BigInt to do something useful for symbolic mathematics, for which they are clearly not designed. Now we have to bend here and there and make everything everywhere (!) else worse, just so that Modules or Polynomials or whatever over BigInt can work, which no one should be using in the first place.

Are we really going to make Nemo experience worse at the REPL by requiring people to do Nemo.div(a, b) instead of div(a, b)?

[wbhart]

In AbstractAlgebra. It was already fixed there. We've been using different internal versions of a number of functions that Julia messes up for ages. All I'm doing here is rolling that change out to Nemo. I even went to the trouble of finding an explicit bug that is caused by it (one of many I am sure).

Here is one ticket where we discussed many permutations of things that could fix the div/rem/mod/divrem issues (there are others):

Nemocas/AbstractAlgebra.jl#101

There simply aren't any other workable solutions. I've exhausted myself trying to find one!

[wbhart]

You can export the Nemo version if you want, but it won't match Base. That means AbstractAlgebra and Nemo will behave differently for the user. If that's what you want, I'll make the change.

Or, as I offered above a few days ago, we can not export it at all.

The only thing I care about is that the definition of div and divrem inside Nemo match the definitions inside AbstractAlgebra, so that the generic code in AbstractAlgebra works with the Nemo integer type, AND so that it matches what Flint is doing everywhere, so I don't have to rewrite most of Flint.

The other option is to rename div and divrem EVERYWHERE to something else, e.g. quo and quorest (you can't use quo and quorem because rem does not agree with the definition that we require). If you have a better idea, let's hear it, but it would have been much better to bring it up months ago, before we changed AbstractAlgebra.

[wbhart]

This has nothing to do with BigInt! The definition for BigInt does NOT match the definition we use everywhere in AbstractAlgebra, which is the definition Flint uses.

The problem is, it doesn't match what Nemo uses (internally) hence the bugs. That's precisely what I'm trying to fix.

I can't change how Julia handles BigInt's outside of AbstractAlgebra, but we have changed it so it is usable inside AbstractAlgebra. Are we really going to go back and change everything all over again because now after months you suddenly don't like it again!?

[wbhart]

Here are your three options:

1. Change AbstractAlgebra back so it uses the Julia definitions of div and divrem. This will sacrifice uniqueness of HNF. Change every function in Flint which does a division with remainder over Z to return the same values as AbstractAlgebra. This includes rewriting the HNF in Flint, amongst many other things.

2. Rename div and divrem to something else that doesn't use the words div, rem or mod. You may want to think carefully before using quo as it's used for other things. And what are you going to call the remainder? rem is out. % is out. Then you will still have to fix the definition in Nemo, and it will be exactly what I'm trying to change it to!

3. Do what I propose here which is that we have our own definitions of div, rem and divrem in AbstractAlgebra/Nemo/Oscar. You can choose whether to export the Nemo definition, a definition which agrees with Base or not export anything. I don't care which.

All three options suck. But the one that is consistent with how we've changed AbstractAlgebra, and which requires the least rewriting of code is the one I picked, namely 3.

I really don't like 2 because no one wants to use a CAS where the Euclidean quotient and remainder functions are named something bizarre.

[thofma]

Just to make sure, inside Nemo div still has the same behavior as before, right?

So the only annoyance is that inside Nemo, for a, b of type fmpz the expression div(a, b) has a different meaning than in the REPL, right?

(P.S.: Exporting AbstractAlgebra.div would not work, as then the user is greeted by `Both Base and AbstractAlgebra export div. Use Base.div or AbstractAlgebra.div)

[wbhart]

For your first question: sort of. If you are doing div of anything other than fmpz integers, yes. If you are doing div of integers, it now uses the definition that is the same as the one used inside AbstractAlgebra and internally in Flint, which it wasn't before. So it does change the behaviour of div for fmpz inside Nemo.

The answer to your second question is yes. But as I said, I don't care what behaviour you pick for the REPL. It can be the same or different or nonexistent. I just chose what would give exactly the same behaviour as AbstractAlgebra currently has (which we have no choice about, as it depends on what Julia does in Base).

As for your P.S. yes, correct, we cannot export the "better" definitions or the user will be greeted with unhappiness.

[thofma]

Great. This means that div for fmpz (finally) agrees with div in Magma.

I think what I meant to say is that we have no choice for my second question.

[wbhart]

No, we have a choice. You can choose to make the Base definition for fmpz match the Nemo definition. It just won't match what Base does for its own integers at the REPL (and hence won't match the behaviour of AbstractAlgebra).

I am not 100% sure my choice is the better one on that issue. We can certainly change this in the future if we decide it is the wrong choice. However, we've lived with this in AbstractAlgebra for ages now and so far it seems to work, so I am cautiously optimistic.

[wbhart]

For completeness, here is another discussion we had on the same issue, basically going back and forth over the same things and not really reaching a workable conclusion. I don't suggest re-reading it now. Most of what is there is just wrong. The later ticket is "more correct", but still incorrect. I'm only giving it here for future reference.

#699

I think there is another somewhere, where I spelled out all our options and what is wrong with each of them. I'll try to find it, but so far am having some trouble doing so.

[wbhart]

This is related, including for completeness:

#688

[wbhart]

I found the other "discussion" we had (basically me explaining over and over the same points in different ways). It was an email discussion. I reproduce part of it here (with some modifications so it flows better), since it is not otherwise online (I have made sure there are no personal details). I have added emphasis ******* at one point below:

I decided to have a look at #699

To remind you of the problem, so you don't have to read all that:

• In Julia rem(-2, 3) = -2 != 1 = rem(1, 3), i.e. there is no unique
  set of representatives

• Julia makes div, rem and divrem a triple that agree with one another over Z

• Julia also has mod, which is defined differently to rem over Z, i.e.
  mod(-2, 3) = 1 = mod(1, 3)

• Nemo agrees with all Julia's definitions over Z

• Internally in AbstractAlgebra we define divrem to return the
  quotient and remainder that you would get if you had used mod, not rem

• This issue has nothing to do with the modulus being positive or not.
  Note rem(-2, 3) != rem(1, 3) and the modulus is positive!

My reasoning in #699 is that if we eventually want to have a
meaningful Euclidean ring interface (which includes Z), we should
switch over to using the functions div, rem and divrem for Euclidean
division since we can't redefine div or divrem over Z and Julia
defines div and divrem to return the result you would get if you used
rem not mod for remainder.

At first this seemed like a simple task. All I have to do is remove
our special definition of div and divrem inside AbstractAlgebra and
change every case where we currently use mod for Euclidean division
over to use rem instead. Every time we define a function called mod,
just rename it to rem.

In fact the only place in AbstractAlgebra where we define a mod
function other than Z is univariate polynomials. So this is all
actually easy to do.

The problem is the tests fail. There are two places where we really
rely on using mod over Z instead of rem:

• HNF
• residue rings mod Z

In HNF we currently have the nice property that HNF is unique over Z.
And when doing comparisons in Z/nZ you really don't want -2 and 1 to
be different mod 3. Flint also uses mod instead of rem in these two
cases, so there's no point just breaking the unique residues thing in
AbstractAlgebra and weakening all the HNF tests. The results won't
agree with Flint if we do that anyway.

There's clearly no way to fix all this.

Obviously in every ring other than Z we can just set rem and mod to do
the same thing.

This all strongly argues for continuing to do all Euclidean divisions
using our special hacked div and divrem functions in AbstractAlgebra
over Z and keep reducing everything using mod, not rem.

But I have no idea what we do for the general Euclidean interface from
the user perspective (e.g. in Oscar). We certainly can't use div, mod
and divrem as the three Euclidean functions because div and divrem are
defined incompatibly with mod in Julia over Z.

We can't use div, rem and divrem for the reasons I outlined above.

Of course Euclidean division doesn't require a unique set of
representatives. That's not the problem here. The problem is that we
need to eventually design a Euclidean interface, so that all generic
functions that work over Euclidean domains work over any ring that (is
Euclidean and) defines functions for Euclidean division. That means we
need to pick three that are compatible for that interface. The problem
is HNF and residue rings will have to work over rings that satisfy
that interface.

We have two choices:

• lose uniqueness of HNF over Z and lose uniqueness of residues in
  Z/nZ (and polynomials and matrices over Z/nZ, etc) and use div, rem
  and divrem as the Euclidean functions in our Euclidean interface and
  rewrite a good proportion of Flint.

• Tell the user that div, rem, divrem and mod are all functions for
  Euclidean division and that div and divrem are defined in terms of
  rem, but that internally AbstractAlgebra/Nemo/Hecke/Oscar use a triple
  of functions for Euclidean division that are defined in terms of mod
  instead of rem, which can be accessed via Oscar.div, Oscar.mod and
  Oscar.divrem.

(Actually, accessing them via Oscar.blah might not be a good idea
since people won't be able to resist importing Base.div, Base.rem and
Base.divrem into Oscar which will break everything all the time. It
might be better to access them via Oscar.AbstractAlgebra.blah.)

Obviously the first option above is bad and impractical. The second is
the only option remaining as I see it. It's messy, but it's going to
have to do.

A third option would be to define a hacked version of rem over Z in
AbstractAlgebra which does what our current mod does, then use these
internal div, rem and divrem functions as the triple for Euclidean
division everywhere in Oscar.

The reason this is pointless is Julia also defines invmod and we have
mulmod and powmod in addition to this. So just using rem instead of
mod so that the names divrem and rem match is pointless. Changing from
mod to rem would break the nice symmetry between mod, invmod, mulmod
and powmod, all of which are defined in terms of our mod.

But just for fun:

julia> using Nemo

julia> invmod(-2, -3) # julia's invmod, defined in terms of mod it seems
-2

julia> invmod(ZZ(-2), ZZ(-3)) # Nemo's
ERROR: DomainError with -3:
Modulus must be non-negative

I've no idea if that is good or bad.

The upshot of all of this is we are going to have to use
AbstractAlgebra.div and AbstractAlgebra.divrem everywhere in Oscar
without fail, and enforce that on all our developers, otherwise we are
going to get broken code over Z on a regular basis.

I guess in Oscar we can just import the AbstractAlgebra definitions of
div, mod and divrem instead of the Base ones and put a note above the
imports to not import the Julia definitions of these functions. This
means our developers can use div, mod and divrem everywhere in Oscar
without having to qualify it and it'll be the correct one.

Actually, what I said above isn't correct either. Whilst we have a
special definition of div and divrem in AbstractAlgebra so that div,
divrem and mod are a triple of functions for Euclidean division that
agree with Flint, these only apply to Julia integer types inside
AbstractAlgebra.

Inside Nemo we use the Nemo definitions of div and divrem, which agree
with the Julia definitions (and rem). This means that if any Nemo (or
Oscar) code ever uses generic Abstract algebra code which uses div or
divrem and applies it to ZZ, it will get the wrong (Nemo) definitions.

It's just luck that this hasn't happened to us so far in any critical
way. This is likely because we use Flint's hnf over Z (which already
does reduction using the equivalent of Julia's mod instead of rem) and
for Z/nZ we have Flint modules that do everything for us, rather than
using generic AbstractAlgebra code.

Even over Z/nZ for multiprecision n we'd be using Flint's fmpz_mod for
(generic AbstractAlgebra) matrices over Z/nZ, so we probably always
get the right answer.

At the present moment I can't see any way to trip up Nemo by making it
call the wrong divrem or div, but that's more or less an accident, not
by design.

We will of course have the same problem in Oscar.

I think the only strictly correct way of doing things would be to
define div and divrem for fmpz internally in Nemo to agree with the
internal AbstractAlgebra definitions for Julia integers and only
export the current definition for use by the user (for consistency
with Julia).

Then Oscar could internally import both the internal
AbstractAlgebra and Nemo definitions of div and divrem so that the
results will always internally be consistent with answers you'd get
from Flint, but again from the user's point of view the Base
definitions of div and divrem for fmpz would be the ones exported by
Nemo (which agree with Julia).

Hecke would have to do the same thing, i.e. not import div and divrem
from Base but from Nemo.

The biggest risk would be code in Nemo/Hecke/Oscar that is written for
Z/nZ but lifts to Z, uses divrem and then uses the result expecting
the result to be a non-negative integer. I can't think of a
computation that does that. But such a computation could certainly
exist.

I can't think of a way to modify the Oscar documentation to make a
note of all this. I think I will wait until I write the documentation
for Z/nZ to explain that reduction is always done using the equivalent
of Julia's mod. We'll also just have to take care to explain when we
get to setting out the Euclidean domain infrastructure that the
functions that need to be defined are div, divrem and mod and that
these need to be compatible, and that unless we fix Nemo in the way I
suggested above, ZZ is not a Euclidean domain for the purposes of
Oscar.

Actually, this is not a good example:

"The biggest risk would be code in Nemo/Hecke/Oscar that is written for
Z/nZ but lifts to Z, uses divrem and then uses the result expecting
the result to be a non-negative integer. I can't think of a
computation that does that. But such a computation could certainly
exist."

Of course if you lifted from Z/nZ to Z the result would actually be
non-negative. So divrem would do exactly what you expected.

It's actually hard to come up with a situation where you would end up
running into trouble. But I know we will eventually with the current
precarious inconsistency.

I found an example. The invariant factors function snf is generic for
modules over ZZ. It uses divrem internally.

So the following test function for modules over ZZ fails if you use
Nemo instead of AbstractAlgebra (so that ZZ is fmpz):

julia> for R in [ZZ, QQ]
for iter = 1:100
M = rand_module(R, -10:10)

            I, f = snf(M)

            K, g = kernel(f)

                     @test length(invariant_factors(K)) == 0

            m = rand(I, -10:10)

            @test m == inv(f)(f(m))
                  end
      end

Test Failed at REPL[17]:13
Expression: m == (inv(f))(f(m))
Evaluated: (-1, -3) == (1, -3)

I will open a ticket explaining the issue. Then we can decide if we
want to take the issue seriously.

I think I have a very simple way to explain the problem so that it is
absolutely clear:

• Nemo itself uses two different residue systems for Z. In some places
  it uses the equivalent of Julia's mod (least negative residue system)
  and in other places it uses the equivalent of Julia's rem (residue has
  same sign as dividend).

• It uses the former for compatibility with Flint and the latter for
  compatibility with Julia (both for good reasons which we have
  discussed ad infinitum before and agreed upon).

Quite obviously, people write code based on Nemo with the assumption
that even if one or the other residue system is used, that both are
not used at the same time! As my example shows, this is a natural
assumption. It is the assumption I made when writing the module code!

The only way to fix this is precisely what I proposed. Internal to
Nemo divrem and div must be compatible with mod, i.e. least positive
residue and external to Nemo they must be compatible with rem (least
absolute value with same sign as dividend). The alternative, to
completely rewrite about 100,000 lines of Flint code, is absolutely
out of the question.

There's no way you can make a system consistent if it changes the
Euclidean system half way through a computation (which it necessarily
does at present).

Some additional comments:

• all these problems have nothing to do with smod. This is not used anywhere that is
  relevant to this problem. rem is not smod (symmetric remainder
  system).

• they have nothing to do with hnf. Both Flint and AbstractAlgebra
  agree on the output of HNF over Z because of careful coding by others.

• this has nothing to do with the sign of the modulus. Even if all
  moduli everywhere were positive, we'd still have this exact same
  problem.

• these issues have nothing to do with canonical_units.

• as far as I can tell, Z/nZ is not relevant here. Both Flint and
  AbstractAlgebra agree on computations involving Z/nZ, though only by
  making the precise changes to AbstractAlgebra that I am now proposing
  we make to Nemo!

• Lifts aren't relevant at all. If
  I lift from Z/nZ to Z the result will always be least non-negative
  remainder and therefore is not contributing to the problem. The lift
  always takes d.data essentially, and we already agreed on d.data of a
  Z/nZ because that decision was made by Flint.

---

After this the discussion wandered into other (private) matters and various other irrelevancies.

[thofma]

Sorry for being a bit slow yesterday and rehashing the discussion. I agree with all you say and the fix you propose is the right approach.

I am also not 100% sure if is better to have Base.div(a::fmpz, b::fmpz) == Nemo.div(a, b) or Base.div(a::fmpz, b::fmpz) == Base.div(BigInt(a), BigInt(b)), but I am happy with the version that you chose. If we don't like it, we can change it later.

[wbhart]

The disadvantage of the first option is obviously that it may confuse people when debugging their code. They run something in the REPL and it gives different results to what is happening in the library.

The disadvantage of the second option is inconsistency with Julia and AbstractAlgebra which is not so nice for people using the system from the REPL.

I'm not sold on either solution, but at least the one I've chosen is consistent with something and easy to remember. But yeah, we can change it later if we want. It's not that critical from my perspective.

One thing the discussion yesterday highlights though is that we need a much better system for making proposals, updating proposals based on new considerations, making decisions and altering decisions if we come back to them later and feel that we've changed our minds.

@fingolfin has mentioned maybe using a PEP-like approach. But I think we should only use it for critical/breaking decisions on which there is not an obvious way forward that we agree on. And I think we should never make arbitrary decisions, but should instead let a preliminary decision sit for a while until everyone has forgotten about it, then revisit it to make a final decision. And of course a decision should only be made once the pros and cons are fully understood. In some cases a PEP should come with some prototype code. And as we've already discovered, PEP's should only be written when there is some chance of them being accepted.

The problem with the divrem/div stuff is I've been trying very hard to fix it for ages, but kept learning new reasons why my proposed approaches would not work. This got very confusing over a period of months, and eventually everyone got sick of it and perhaps even tuned out. So when it came time to make a decision, no one paid any attention. I therefore assumed no one really cared enough to complain.

There's much more div/divrem/mod/rem craziness still to come, but I think this is the last change where there are no perfect options. Most of the rest of the stuff has to do with consistency on when we raise exceptions, which at least @fieker and I seem to have been agreeing on. This div/divrem change is also the last really urgent change. If we left it much longer, we'd never be able to fix it.

[wbhart]

Some time next week I will try to see what breakage there is in Hecke, Singular and Oscar in the event this is merged in the future. I know some of the tests in Oscar will fail. But these can probably be removed and put into Nemo anyway. I plan to significantly improve the Nemo tests and docs in the coming months.

People can choose to cut'n'paste them for Oscar if they want. That's fine by me. But most of the low level docs and tests actually belong in the low level code, I think. Oscar really just needs to test that things are working mathematically the way the Oscar developers actually want them to work. They shouldn't have to test the functionality to see if it is broken or consistent, etc. That should be the job of the people writing the lower level code in the first place.

[thofma]

Don't worry about Hecke. I will take care of it.

[wbhart]

Ok.