Nemocas / Nemo.jl

Julia bindings for the FLINT number theory C library
http://nemocas.github.io/Nemo.jl/
Other
191 stars 58 forks source link

ERROR: MethodError: no method matching //(::fmpq_mat, ::fmpq) #1281

Closed simonbrandhorst closed 2 years ago

simonbrandhorst commented 2 years ago

I would like to have this. It is inconvenient.

wbhart commented 2 years ago

// is for constructing a fraction from elements of an integral domain. Surely you want divexact.

simonbrandhorst commented 2 years ago

Coming from mathematics I expect that matrix(over a field)/field_element just works . For instance Magma allows it, sage allows it and I am sure (but did not test) that Mathematics, matlab, maple etc. allow it:

> G:= Matrix(Rationals(),1,1,[1]);
> G/2;
[1/2]

It is similar to why I prefer 2*M over multiply(2, M) or something similar.

simonbrandhorst commented 2 years ago
julia> [1 2;]//2
1×2 Matrix{Rational{Int64}}:
 1//2  1//1
wbhart commented 2 years ago

This is not doing what you think it is. It is taking a matrix over the integers and it is creating the matrix over the fraction field with given denominator.

If the // operator is to be used for "division" (since Julia uses / for floating point division) then what operator should we use for constructing the object whose type represents the fraction object?

E.g. should (x^2 + x)//(x + 1) return x or should it return x//1? Those have different types, even if they are mathematical able to be identified.

simonbrandhorst commented 2 years ago

For consistency with

julia> 2//2
1//1

It should probably return (x^2 + x)//(x + 1). But let me be more precise. I propose to set // to scalar divison when the base ring is a field. Then there is no ambiguity. But O.K. I see now that this is possibly a change in the current behavior.

julia> R,x = PolynomialRing(QQ,"x")
(Univariate Polynomial Ring in x over Rational Field, x)
julia> 2*x
2*x
julia> (2*x)
2*x
julia> (2*x)//2
x
julia> parent((2*x)//2)
Fraction field of Univariate Polynomial Ring in x over Rational Field

julia> (2*x)//QQ(2)
ERROR: MethodError: no method matching //(::fmpq_poly, ::fmpq)
wbhart commented 2 years ago

It's not a good idea to change the existing behaviour of the system, especially since it would mean introducing a new operator for constructing fractions in that special case.

The existing behaviour has been the behaviour for more than 6 years. I don't believe it is practical to make changes to suit individual taste, especially when that change introduces inconsistency of meaning in the system.

Unfortunately, computers care not just about the mathematical meaning of operations, but also the type. There are very real inconsistencies in other systems because of these sorts of changes, which I've personally long resisted. The different kinds of division must be distinguished in a typed language.

simonbrandhorst commented 2 years ago

//(::fmpq_poly, ::fmpq_poly) should return an element in the fraction field whereas //(::fmpq_poly, ::fmpq) should return an fmpq_poly Can you give an example of how this could be dangerous/inconsistent? Then // is still the operator to produce a fraction. And similar for vector spaces over a field (e.g. vectors, matrices, polynomials over a field).

I dare say that my taste here is shared by most mathematicians - so the system diverges from the users expectations and we should have a good reason for it.

julia> M
2×2 Matrix{Rational{Int64}}:
 1//1  2//1
 2//1  1//1

julia> M//(1//2)
2×2 Matrix{Rational{Int64}}:
 2//1  4//1
 4//1  2//1

So here we also seem to diverge from julia?

(If this has been discussed (and decided) over and over again, I can live with that,)

thofma commented 2 years ago

At the moment if two ring elements x and y do not have a common parent, all binary operations f satisfy the property that f(x, y) == f(x, parent(x)(y)) respectively f(x, y) == f(parent(y)(x), y). I think we should keep it that way so that x//y is always the same as x//parent(x)(y). I personally don't want x//2 to be different from x//(2*x^0). I would assume that there are users who would welcome this as well.

And yes, we purposefully diverge from julia. I guess you don't want us to introduce

julia> det([1 2; 3 4])
-2.0

or

julia> (2)/(2)
1.0

julia> (2//1)/(2//1)
1//1

Julia linear algebra is not a consistent system and does not provide useful default behaviour unless you live in a numerical bubble and consider everything to be a subset of floating point numbers.

simonbrandhorst commented 2 years ago

O.K. so basically we apply promotion. Then I can see the point.

wbhart commented 2 years ago

Correct. The coercion principle Tommy outlined is precisely the reason it was designed the way it is. You won't have a consistent computer algebra system if you don't follow that rule.

It might help to know that this is not arbitrary. I extensively studied the coercion systems of Sage, Magma, Pari/GP and Julia itself about 6 years ago to see what the different possible approaches are before making that very deliberate decision.

The issue becomes particularly acute as the coercion system of the CAS becomes more sophisticated, which it will inevitably.

fingolfin commented 2 years ago

The real problem here seems to be that there is no method matching //(::fmpq_mat, ::fmpq_mat). If there was, the example would work, for the same reason x * y works when x is an fmpq_mat and y is an fmpq.

Yet I can write x * y^-1 and it works in all combinations of fmpq and fmpq_mat....

thofma commented 2 years ago

For rings (non-fields), // creates something in the fraction field. For matrices this of course makes no sense but I don't think it is worthwhile to introduce a "mismatch" between //(fmpq_poly, ::fmpq) and //(fmpq_mat, ::fmpq). Sometimes you would stay in the same parent and sometimes not. I do not understand the "yet" part and the "...."? Can you spell out what you mean?

wbhart commented 2 years ago

I presume @fingolfin is trying to do this.

julia> M = matrix(QQ, [1 2; 3 4])
[1   2]
[3   4]

julia> divexact(M, M)
[1   0]
[0   1]