Nemocas / AbstractAlgebra.jl

Generic abstract algebra functionality in pure Julia (no C dependencies)
https://nemocas.github.io/AbstractAlgebra.jl/dev/index.html
Other
162 stars 62 forks source link

zero(f::PolyElem) is producing elements with the wrong parent #897

Closed thofma closed 3 years ago

thofma commented 3 years ago

From https://github.com/oscar-system/Oscar.jl/pull/469:

julia> Qx, x = PolynomialRing(QQ, "x", cached = false)
(Univariate Polynomial Ring in x over Rationals, x)

julia> parent(zero(x)) === parent(x)
false

I don't think this is intented.

Edit: Thanks @tthsqe12.

tthsqe12 commented 3 years ago

I think you mean this

julia> using AbstractAlgebra

julia> Qx, x = QQ["x"]
(Univariate Polynomial Ring in x over Rationals, x)

julia> parent(zero(x)) === parent(x)
true

julia> Qx, x = PolynomialRing(QQ, "x", cached = false)
(Univariate Polynomial Ring in x over Rationals, x)

julia> parent(zero(x)) === parent(x)
false
wbhart commented 3 years ago

This is apparently intended behaviour as it is actually tested for (before I touched it).

zero calls similar which makes a polynomial over the same base ring, and not necessarily with the same parent. You'll only get the same parent if it was cached.

I was surprised by this myself, but tests failed when I tried to do it the other way.

On Thu, 10 Jun 2021 at 17:48, tthsqe12 @.***> wrote:

I think you mean this

julia> using AbstractAlgebra

julia> Qx, x = QQ["x"] (Univariate Polynomial Ring in x over Rationals, x)

julia> parent(zero(x)) === parent(x)true

julia> Qx, x = PolynomialRing(QQ, "x", cached = false) (Univariate Polynomial Ring in x over Rationals, x)

julia> parent(zero(x)) === parent(x)false

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/Nemocas/AbstractAlgebra.jl/issues/897#issuecomment-858736469, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAB3HJX57NQHIO3O55EDF7DTSDNGHANCNFSM46OZLWNA .

thofma commented 3 years ago

Hm, I don't think this behavior is useful and it seems that you agree. I will have a look tomorrow.

tthsqe12 commented 3 years ago

If this is the indended behaviour of zero, then alot of our stuff is broken.

wbhart commented 3 years ago

@tthsqe12 No it's not. We scarcely use this zero if at all. You might be thinking of zero(R) for a ring R, which is widely used.

thofma commented 3 years ago

I cannot imagine a situation where this behavior of zero would ever be useful!

wbhart commented 3 years ago

@thofma I think this is a necessary consequence of the decisions we made for similar way back when we sorted that out.

I am too tired right now to piece together exactly why this is needed. But it has to do with the fact that similar isn't supposed to return everything the same, e.g. if you have a poly over number fields and you want a similar one over ZZ you don't want a generic polynomial but a Nemo one if you are using Nemo.

We defined zero in terms of similar, for reasons I don't precisely recall. Perhaps we can revisit that decision, but I do recall we did it very deliberately.

wbhart commented 3 years ago

There's a partial explanation of the rabbit hole here:

https://github.com/Nemocas/Nemo.jl/issues/1045

which includes a link to:

https://github.com/Nemocas/Nemo.jl/issues/651

Basically I think it may be done this way so that parent(zero(x)) returns the right thing when mapping coefficients and the like.

I agree that I don't like the behaviour, but we may be stuck with it.

thofma commented 3 years ago

We added the function here https://github.com/Nemocas/AbstractAlgebra.jl/pull/843. I can't find the discussion of why we would this behavior though. Do you remember where it is?

thofma commented 3 years ago

I will try zero(x) = zero(parent(x)) and see what goes wrong.

wbhart commented 3 years ago

Yes, we added it there, but that is modeled on other functions which already existed, e.g. for matrices and so on. The links I give should contain hints as to why it is like this.

I'm happy to change the behaviour of all zero(x)'s in the system to what you want, but the current test code will definitely break. We will have to come up with another way to handle some tricky issues, I believe.

wbhart commented 3 years ago

You have to change it for all types it is defined for to see it break something, I believe, including matrices from memory.

wbhart commented 3 years ago

Here is more required reading:

https://github.com/Nemocas/AbstractAlgebra.jl/issues/515

https://github.com/Nemocas/AbstractAlgebra.jl/issues/490

thofma commented 3 years ago

For matrices the parent discussion is not existent though.

wbhart commented 3 years ago

I think I understand how this all came about now. The function zero started with matrices. We initially wanted to be able to make a matrix similar to one that already existed, but with a different set of dimension, say, or a different base_ring. Thus similar was born.

But there were some problems. Firstly, suppose you have a matrix over a number field and you want to create a similar matrix over ZZ instead. The problem is, a matrix over a number field is generic, and so you definitely don't want the matrix over ZZ that gets returned to also be generic. You want it to be an fmpz_mat (say) if you are using Nemo.

We decided that similar should create the matrix most appropriate to the package in use. Therefore, similar considers only the base_ring R and not the parent of the matrix you pass to it. In essence it creates a matrix of type dense_matrix_type(R).

Originally, similar created a zero matrix by setting all entries to zero.

But at some point we noted that it was inefficient to create a matrix this way since you usually go back and change all the entries anyway and throw away all the zero entries you just created. Similar was changed to create a matrix with uninitialised entries.

zero was then introduced to create a matrix with the entries all set to zero (i.e. what similar originally did), but would otherwise act in precisely the same way as similar.

Now here's the thing. You could make zero steal the parent of the original matrix and use that if the base_ring requested is the same as that of the original matrix. But then the behaviour of zero is inconsistent. Start with a generic matrix over number fields and ask it for one over ZZ and it would give you an fmpz_mat. But start with a generic matrix over ZZ and ask it for one over ZZ and it will return a Generic.Mat{fmpz}. That's not good design.

Now, looking through our code so far, it seems we only ever use zero(M) for a matrix M in a couple of different scenarios:

1) we want a zero matrix of exactly the same kind as the input matrix M, e.g. because we want to return a zero matrix for a corner case of some algorithm 2) we want a zero matrix of exactly the same kind as the input matrix M but with a different number of rows and columns (e.g. nullspace) 3) we use it in cat, though I don't actually know what this code is doing

If there are actually any really funky uses of zero, they probably lie in Hecke, not AbstractAlgebra or Nemo.

Here is what I think. If we make zero use the same parent as the input (though allow different dimensions in the case of matrices) then it should not accept a (base) Ring at all. Otherwise the function is not consistent with itself (which is likely to result in bugs in more complex uses down the track) and with similar.

We actually don't need zero to take a base ring, as we already have zero_matrix for this precise purpose. The only advantage of the current zero over this is that you don't have to specify the number of rows and columns if you have a matrix of the right size already that you can pass to zero to steal the dimensions from.

In short, I would support the following:

I definitely don't support zero stealing the parent only in the case where the base_ring requested is the same as that of the supplied matrix/poly etc. That is inconsistent both with similar and with itself and likely to lead to bugs.

thofma commented 3 years ago

This is about polynomials and not matrices.

For matrices you can do M + zero(M) and it just works. For polynomials p + zero(p) does not work, which makes this a pretty useless function which should never ever be used anywhere. I am not talking about similar. I am just talking about the zero function.

julia> using AbstractAlgebra

julia> Qx, x = PolynomialRing(QQ, "x", cached = false)
(Univariate Polynomial Ring in x over Rationals, x)

julia> x + zero(x)
ERROR: Incompatible polynomial rings in polynomial operation
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:33
 [2] check_parent
   @ ~/.julia/packages/AbstractAlgebra/2tNNk/src/generic/Poly.jl:76 [inlined]
 [3] check_parent
   @ ~/.julia/packages/AbstractAlgebra/2tNNk/src/generic/Poly.jl:75 [inlined]
 [4] +(a::AbstractAlgebra.Generic.Poly{Rational{BigInt}}, b::AbstractAlgebra.Generic.Poly{Rational{BigInt}})
   @ AbstractAlgebra.Generic ~/.julia/packages/AbstractAlgebra/2tNNk/src/generic/Poly.jl:498
 [5] top-level scope
   @ REPL[4]:1
wbhart commented 3 years ago

Yes, but I modeled the zero for polynomials on that for matrices. And zero for matrices is using similar and has identical semantics except that it intialises all the entries to zero.

I really don't like the idea of having zero for matrices do one thing and zero for polynomials behave in a completely different way.

I understand the problem you are trying to solve. The question is, how do we solve it consistently?

wbhart commented 3 years ago

And the only reason that works for matrices is that we don't check parents in matrices because matrices don't store their parents.

thofma commented 3 years ago

This problem does not exist for matrices since they don't have parent checks. Otherwise it would be equally useless.

Can you describe where the current zero function for polynomials would be useful or where it could be used safely? It seems it must not be used in library code.

wbhart commented 3 years ago

I agree the current zero is not safe. I'm not arguing to keep the current behaviour. I'm arguing for finding a solution that isn't inconsistent with how zero works for matrices.

wbhart commented 3 years ago

I mean, it is not as if this is particularly satisfying:

julia> using AbstractAlgebra

julia> S = MatrixSpace(ZZ, 2, 2; cached=false)
Matrix Space of 2 rows and 2 columns over Integers

julia> M = rand(S, -10:10)
[1   7]
[8   8]

julia> parent(M) == S
false
rfourquet commented 3 years ago

(Just lurking here) Regarding your above example:

julia> parent(M) == S
false

Couldn't == be defined for matrix spaces? (this would check for equal dimensions and equal base_ring)

wbhart commented 3 years ago

One place where the current behaviour for polynomials might actually be useful is in converting Singular polynomials over Nemo coefficient rings to Nemo polynomials. Not that this is a good argument for keeping the current behaviour.

thofma commented 3 years ago

Yes, this whole similar thing was added for matrices in AbstractAlgebra with no regard to the parent system whatsoever. Now we have to live it.

Here two options to "fix" it for polynomials.

  1. If cached == true use the parent of the input. If cached == false create a new one like we do now.
    zero(p::PolyElem; cached = true) = !cached ? similar(p, cached = false) : zero(parent(p))
  2. Remove the cached argument from zero(p::PolyElem,...) (the one without the base ring argument) and just do
    zero(p::PolyElem) = zero(parent(p))

    (It think there is little use of zero(p, cached = false))

wbhart commented 3 years ago

@rfourquet Perhaps that will work. We might even need it, because MatrixAlgebras could be used in generic code that checks parents.

wbhart commented 3 years ago

I think we should define zero_polynomial(R::Ring) just as we have zero_matrix(R::Ring, r::Int, c::Int) and remove the ring argument to zero in both cases.

rfourquet commented 3 years ago

Couldn't == be defined for matrix spaces?

Actually this == logic is already implemented in check_parent(x, y) for matrices. So == could be defined, and check_parent could be modified to call out to ==.

wbhart commented 3 years ago

Yes, this whole similar thing was added for matrices in AbstractAlgebra with no regard to the parent system whatsoever. Now we have to live it.

That is simply not true. We agonised over the design, specifically with regard to the parent system. It is absolutely by design that it disregards the parent and only looks at the base ring. And if you look at the links I gave, it was actually you that pointed out the deficiencies in the original implementation of similar that forced us to make it operate the way it does.

I don't think the similar system is broken.

Where things broke is that we used zero for the version of similar that fills the entries with zeroes. That was a poor design choice because people now want to use zero a different (much more logical) way.

As I say, I think the best solution here is to remove the ring argument to zero altogether and force it to always use the parent (both for matrices and polynomials), differing dimensions excepted of course.

wbhart commented 3 years ago

Does Hecke ever use the ring argument to zero for matrices?

thofma commented 3 years ago

The same problem applies to similar. The call p + similar(p) is also broken.

wbhart commented 3 years ago

Yeah I know. This really isn't what similar is for. Similar is for creating the best type possible. E.g. if you pass in a generic polynomial over ZZ and ask for similar it will give you an fmpz_poly. This is our deliberate design choice.

Besides, I really don't think we can expect objects without cached parents to behave as though they had them.

Parent checking really serves the purpose that a type checker would serve if the type system allowed types parameterised on values. Allowing objects to not have consistent parents is really doing an end-run around that system. We should expect breakage. The question is just how to make it maximally useful and yet still consistent.

wbhart commented 3 years ago

One possible solution to the whole mess would be to define comparison of parents. E.g. for polynomials:

==(a::PolyRing{T}, b::PolyRing{T}) where T = var(a) == var(b) && base_ring(a) == base_ring(b)

function check_parent(a::PolyElem{T}, b::PolyElem{T}) where T
   if parent(a) !== parent(b)
      if parent(a) != parent(b)
         error("Parents not the same")
      end
   end
end

This way most of the time the comparison would be trivial. But it would still pass if the parents were not the same but the parameters were the same. Of course this would slow down generic code which didn't cache parents. But maybe it is worth it for the sake of consistency. I don't know.

wbhart commented 3 years ago

Of course this would be uber slow in the case of comparing a modulus for a residue ring. So maybe this would just be asking for trouble.

wbhart commented 3 years ago

On the other hand, we could do the same thing there and compare the moduli first for object equality. In the case where the second polynomial was created using similar it would have used precisely the same base_ring and so the comparison would be quick.

tthsqe12 commented 3 years ago

Assuming this is

==(a::PolyRing{T}, b::PolyRing{T}) where T = var(a) == var(b) && base_ring(a) == base_ring(b)

I think this is asking for trouble when arithmetic operations expect identical (===) parents

wbhart commented 3 years ago

Yes, typo on my part.

I'm saying check_parent (which is called by arithmetic operations) would use this equality check if the parents aren't identical.

What situation would that cause a problem for?

wbhart commented 3 years ago

The one thing that would be tricky is all the code we have that does things like:

parent(a) != base_ring(b) && error("Cannot coerce")

This would also have to be changed, though I suppose == could first check === rather than splitting the check as I do above, i.e. check_parent would just call == and nothing else.

wbhart commented 3 years ago

Anyhow, @thofma can you tell me if there is something wrong with my proposal:

The alternative would be parent comparison as @rfourquet and I have outlined.

I'm not in favour of making zero for polys behave differently to zero for matrices. The definition should be conceptually the same in both cases.