inverse of matrices over arbitrary fields stopped working in julia1.5.0-beta1

[jmichel7]

I haved used julia since 0.6, and at all times matrix inverse worked for fields I have defined (finite fields, cyclotomic fields, rational fractions...). It stooped working in 1.5.9-beta1

Here is a MWE where I define a matrix over a finite prime field

   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.4.2 (2020-05-23)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

julia> struct Mod{p} <: Number
          val::Int
          function Mod{p}(a::Integer) where {p}
                  new(mod(a,p))
          end
       end

julia> Base.zero(::Type{Mod{p}}) where p = Mod{p}(0)

julia> Base.:+(x::Mod{p}, y::Mod{p}) where p = Mod{p}(x.val+y.val)

julia> Base.:*(x::Mod{p}, y::Mod{p}) where p = Mod{p}(x.val*y.val)

julia> Base.:-(x::Mod{p}) where {p} = Mod{p}(-x.val)

julia> Base.:-(x::Mod{p}, y::Mod{p}) where p =x+(-y)

julia> Base.inv(x::Mod{p}) where p = Mod{p}(invmod(x.val,p))

julia> Base.:/(x::Mod{p}, y::Mod{p}) where p = x * inv(y)

julia> Base.isless(x::Mod{p}, y::Mod{p}) where p=x.val<y.val

julia> Base.abs(x::Mod)=x

julia> Base.conj(x::Mod)=x

julia> m=Mod{19}.([1 2 3;3 2 1;1 0 0])
3×3 Array{Mod{19},2}:
 Mod{19}(1)  Mod{19}(2)  Mod{19}(3)
 Mod{19}(3)  Mod{19}(2)  Mod{19}(1)
 Mod{19}(1)  Mod{19}(0)  Mod{19}(0)

julia> inv(m)*m
3×3 Array{Mod{19},2}:
 Mod{19}(1)  Mod{19}(0)  Mod{19}(0)
 Mod{19}(0)  Mod{19}(1)  Mod{19}(0)
 Mod{19}(0)  Mod{19}(0)  Mod{19}(1)

This stopped working in 1.5.0-beta1

julia> inv(m)
ERROR: MethodError: no method matching AbstractFloat(::Mod{19})
Closest candidates are:
  (::Type{T})(::T) where T<:Number at boot.jl:716
  AbstractFloat(::Bool) at float.jl:258
  AbstractFloat(::Int8) at float.jl:259
  ...
Stacktrace:
 [1] float(::Mod{19}) at ./float.jl:277
 [2] norm at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/generic.jl:648 [inlined] (repeats 2 times)
 [3] #generic_lufact!#143 at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lu.jl:143 [inlined]
 [4] #lu!#142 at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lu.jl:130 [inlined]
 [5] #lu#144 at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lu.jl:273 [inlined]
 [6] lu at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lu.jl:272 [inlined] (repeats 2 times)
 [7] inv(::Array{Mod{19},2}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/dense.jl:809
 [8] top-level scope at REPL[14]:1

But the code for inverting matrices over an arbitrary field must still be somewhere in the library since inverting matrices of rationals still work. How do I use the code for inverting matrices over arbitrary fields for my own types?

[andreasnoack]

In JuliaLang/julia#34575 we changed the LU to call norm instead of abs. So norm has a fallback for Number but, as your stacktrace shows, it calls float which we probably hadn't anticipated when we thought the change would have limited effect. I can see two possible solutions here. Either we revert JuliaLang/julia#34575 or we stop promoting in norm(::Number). The latter would mean that typeof(norm(1)) would be different from typeof(norm(1:1)) but I'd be fine with that.

[JeffBezanson]

I vote for just removing the call to float. Makes sense to me for norm(x::Number) to be the same type as x.

[lostella]

I vote for just removing the call to float. Makes sense to me for norm(x::Number) to be the same type as x.

Would that mean that norm(1.0 + 1.0im) == 1.4142135623730951 + 0.0im?

Edit: probably not, that wouldn’t rely on the Number fallback

[andreasnoack]

I'd forgotten #593. We'll have to revert JuliaLang/julia#34575 instead of changing the return type of norm(::Number).

I think we made the wrong call in JuliaLang/julia#1875. The norm code is such a mess because of all the functions bundled together. It's frustrating that we can't handle a useful case like the one in JuliaLang/julia#34575 because we need to consider the situation in #593 which has zero chance of showing up in real code.

[JeffBezanson]

Another option is to have norm(x::Number)::typeof(x), but add methods that promote integers to floats. In general anything defined on ::Number should have really minimal dependencies. We should not be sprinkling calls to float(x), especially not to fix issues specific to 2's complement signed integers.

[jebej]

Wouldn't the simple fix for this issue be to define Base.norm(x::Mod)=x, instead of Base.abs(x::Mod)=x?

[jebej]

Ah it's a compatibility issue I suppose.