Smith normal form seems broken for matrices of rational numbers

[nsajko]

The same works with Nemo:

julia> using AbstractAlgebra

julia> const nums = -10:10
-10:10

julia> const dens = ((-5:-1)..., (1:5)...)
(-5, -4, -3, -2, -1, 1, 2, 3, 4, 5)

julia> r() = rand(nums)//rand(dens)
r (generic function with 1 method)

julia> m = [r() for _ ∈ 1:4, _ ∈ 1:6]
4×6 Matrix{Rational{Int64}}:
  -1     9     5//2  -5    -8//5   1//4
 10//3   3      0    -1      0    -1//3
  3//2  4//3    3    -2     -1      0
  -1     2    -9//2  2//3  -3//2    2

julia> a = matrix(QQ, m)
[-1//1   9//1    5//2   -5//1   -8//5    1//4]
[10//3   3//1    0//1   -1//1    0//1   -1//3]
[ 3//2   4//3    3//1   -2//1   -1//1    0//1]
[-1//1   2//1   -9//2    2//3   -3//2    2//1]

julia> snf(a)
ERROR: MethodError: no method matching mul_red!(::Rational{BigInt}, ::BigInt, ::Rational{BigInt}, ::Bool)

Closest candidates are:
  mul_red!(::T, ::T, ::T, ::Bool) where T<:NCRingElement
   @ AbstractAlgebra ~/.julia/packages/AbstractAlgebra/dsta0/src/NCRings.jl:215

Stacktrace:
 [1] kb_clear_row!(S::AbstractAlgebra.Generic.MatSpaceElem{Rational{BigInt}}, K::AbstractAlgebra.Generic.MatSpaceElem{Rational{BigInt}}, i::Int64, with_trafo::Bool)
   @ AbstractAlgebra ~/.julia/packages/AbstractAlgebra/dsta0/src/Matrix.jl:5348
 [2] snf_kb!(S::AbstractAlgebra.Generic.MatSpaceElem{Rational{BigInt}}, U::AbstractAlgebra.Generic.MatSpaceElem{Rational{BigInt}}, K::AbstractAlgebra.Generic.MatSpaceElem{Rational{BigInt}}, with_trafo::Bool)
   @ AbstractAlgebra ~/.julia/packages/AbstractAlgebra/dsta0/src/Matrix.jl:5376
 [3] _snf_kb
   @ AbstractAlgebra ~/.julia/packages/AbstractAlgebra/dsta0/src/Matrix.jl:5325 [inlined]
 [4] snf_kb
   @ AbstractAlgebra ~/.julia/packages/AbstractAlgebra/dsta0/src/Matrix.jl:5306 [inlined]
 [5] snf(a::AbstractAlgebra.Generic.MatSpaceElem{Rational{BigInt}})
   @ AbstractAlgebra ~/.julia/packages/AbstractAlgebra/dsta0/src/Matrix.jl:5432
 [6] top-level scope
   @ REPL[7]:1

Using snf_with_transform instead of snf results in the same MethodError.

[fieker]

On Fri, Oct 27, 2023 at 06:42:03AM -0700, Neven Sajko wrote: The same works with Nemo: ```julia-repl julia> using AbstractAlgebra julia> const nums = -10:10 -10:10 julia> const dens = ((-5:-1)..., (1:5)...) (-5, -4, -3, -2, -1, 1, 2, 3, 4, 5) julia> r() = rand(nums)//rand(dens) r (generic function with 1 method) julia> m = [r() for _ ∈ 1:4, _ ∈ 1:6] 4×6 Matrix{Rational{Int64}}: -1 9 5//2 -5 -8//5 1//4 10//3 3 0 -1 0 -1//3 3//2 4//3 3 -2 -1 0 -1 2 -9//2 2//3 -3//2 2 julia> a = matrix(QQ, m) [-1//1 9//1 5//2 -5//1 -8//5 1//4] [10//3 3//1 0//1 -1//1 0//1 -1//3] [ 3//2 4//3 3//1 -2//1 -1//1 0//1] [-1//1 2//1 -9//2 2//3 -3//2 2//1] julia> snf(a) ERROR: MethodError: no method matching mul_red!(::Rational{BigInt}, ::BigInt, ::Rational{BigInt}, ::Bool) Closest candidates are: mul_red!(::T, ::T, ::T, ::Bool) where T<:NCRingElement @ AbstractAlgebra ~/.julia/packages/AbstractAlgebra/dsta0/src/NCRings.jl:215 Stacktrace: [1] kb_clear_row!(S::AbstractAlgebra.Generic.MatSpaceElem{Rational{BigInt}}, K::AbstractAlgebra.Generic.MatSpaceElem{Rational{BigInt}}, i::Int64, with_trafo::Bool) @ AbstractAlgebra ~/.julia/packages/AbstractAlgebra/dsta0/src/Matrix.jl:5348 [2] snf_kb!(S::AbstractAlgebra.Generic.MatSpaceElem{Rational{BigInt}}, U::AbstractAlgebra.Generic.MatSpaceElem{Rational{BigInt}}, K::AbstractAlgebra.Generic.MatSpaceElem{Rational{BigInt}}, with_trafo::Bool) @ AbstractAlgebra ~/.julia/packages/AbstractAlgebra/dsta0/src/Matrix.jl:5376 [3] _snf_kb @ AbstractAlgebra ~/.julia/packages/AbstractAlgebra/dsta0/src/Matrix.jl:5325 [inlined] [4] snf_kb @ AbstractAlgebra ~/.julia/packages/AbstractAlgebra/dsta0/src/Matrix.jl:5306 [inlined] [5] snf(a::AbstractAlgebra.Generic.MatSpaceElem{Rational{BigInt}}) @ AbstractAlgebra ~/.julia/packages/AbstractAlgebra/dsta0/src/Matrix.jl:5432 [6] top-level scope @ REPL[7]:1 ``` Using `snf_with_transform` instead of `snf` results in the same `MethodError`. -- Reply to this email directly or view it on GitHub: #1477 You are receiving this because you are subscribed to this thread. Message ID: ***@***.***>

While I agree that the error does not help, what is the expected result? returning a diagonal matrix with 0 and 1 would be legal? What is expected? The snf in general works over euclidean rings, Q is a field, hence euc, but a legal defintion of gcd is 1 for all non-trivial inputs...

[nsajko]

While I agree that the error does not help, what is the expected result?

Presumably the same result as with Nemo? TBH I'm not very interested in this issue, but it came up while I was trying to help someone on the Discourse: https://discourse.julialang.org/t/smith-normal-form/105477

returning a diagonal matrix with 0 and 1 would be legal?

Yeah, there in the thread linked above, someone, who happens to be one of the other developers in this project, I think, pointed out as much.

[fieker]

The core problem is a julia decision:

julia> gcdx(Rational{Int}(1, 2), Rational{Int}(2,3))
(1//6, -1, 1)

julia> typeof(ans)
Tuple{Rational{Int64}, Int64, Int64}

(and similar for big rationals) breaks the assumption that the Bezout coefficients are always in the same ring/ field as the input. In julia this might work as "everything" is sooner or later Real or Complex.
In AA /Nemo this behaviour is bad as e.g. gcdx for polynomials has to return polynomialsm "never" integers.
It works in Nemo as QQ in Nemo is a different type. For Nemo.QQ the Bezout coefficients are rationals, hence the HNF/ SNF code succeeds.