`*` on matrices of custom number types used to work in Julia <=1.10 and does not anymore (now requires `zero` to be defined and maybe shouldn't)

[jmichel7]

The following (minimum definition of addition and multiplication of modular numbers) enabled to multiply matrices in Julia 1.0 to 1.10 (and probably before, I was not using Julia then)

               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.10.9 (2025-03-10)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

julia> struct Mod<:Number
         val::Int
         modulo::Int
         Mod(x::Int,y::Int)=new(x%y,y)
       end

julia> Base.:+(x::Mod, y::Mod)=Mod(x.val+y.val,x.modulo)

julia> Base.:*(x::Mod, y::Mod)=Mod(x.val*y.val,x.modulo)

julia> m=Mod.([4 5;5 4],20)
2×2 Matrix{Mod}:
 Mod(4, 20)  Mod(5, 20)
 Mod(5, 20)  Mod(4, 20)

julia> m*m
2×2 Matrix{Mod}:
 Mod(1, 20)  Mod(0, 20)
 Mod(0, 20)  Mod(1, 20)

but in julia 11 and 12 we get

julia> m*m
ERROR: MethodError: no method matching Mod(::Int64)
The type `Mod` exists, but no method is defined for this combination of argument types when trying to construct it.

Closest candidates are:
  Mod(::Int64, ::Int64)
   @ Main REPL[1]:4
  (::Type{T})(::T) where T<:Number
   @ Core boot.jl:965
  (::Type{T})(::Base.TwicePrecision) where T<:Number
   @ Base twiceprecision.jl:265
  ...

Stacktrace:
  [1] convert(::Type{Mod}, x::Int64)
    @ Base ./number.jl:7
  [2] zero(::Type{Mod})
    @ Base ./number.jl:314
  [3] _rmul_or_fill!
    @ ~/julia-1.12.1/share/julia/stdlib/v1.12/LinearAlgebra/src/generic.jl:169 [inlined]
  [4] _generic_matmatmul_nonadjtrans!
    @ ~/julia-1.12.1/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:1031 [inlined]
  [5] __generic_matmatmul!
    @ ~/julia-1.12.1/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:1027 [inlined]
  [6] _generic_matmatmul!(C::Matrix{Mod}, A::Matrix{Mod}, B::Matrix{Mod}, alpha::Bool, beta::Bool)
    @ LinearAlgebra ~/julia-1.12.1/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:1021
  [7] generic_matmatmul!
    @ ~/julia-1.12.1/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:1011 [inlined]
  [8] generic_matmatmul_wrapper!
    @ ~/julia-1.12.1/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:343 [inlined]
  [9] _mul!
    @ ~/julia-1.12.1/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:328 [inlined]
 [10] mul!
    @ ~/julia-1.12.1/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:297 [inlined]
 [11] mul!
    @ ~/julia-1.12.1/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:265 [inlined]
 [12] mul
    @ ~/julia-1.12.1/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:118 [inlined]
 [13] *(A::Matrix{Mod}, B::Matrix{Mod})
    @ LinearAlgebra ~/julia-1.12.1/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:114
 [14] top-level scope
    @ REPL[5]:1

The problem is the call zero(Mod), which should not be necessary to multiply two matrices. We can define

Base.zero(x::Mod)=Mod(0,x.modulo)

but not zero(Mod) since the type Mod does not contain the information modulo. I find this a serious regression in the generality of Julia. In mathematics, there are lots of similar numbers where zero(x) makes sense but not zero(typeof(x)), and Julia handles recently worse and worse this type of numbers (I found also that LinearAlgebra.det_bareiss(m), which used to work, does not anymore).

I think it is would be worth taking a minimum of care writing generic code in Linear algebra so more of it works with such objects. For example, one could use instead of zero(eltype(m)) a function

zeroentry(m)=isempty(m) ? zero(eltype(m)) : zero(first(m))

which would limit the problem to empty matrices.

[adienes]

I hope it's ok I've taken the liberty of moving & renaming since this is not an issue about Julia 1.11 or 1.12 specifically

[jmichel7]

Its OK, excepted that I am afraid that paradoxically less people care about the question in the LinearAlgebra issues that in the general issues, since here everybody is focused on float matrices and not on general algebra.

Also, the issue is somehow more basic (already apparent in Base.abstractarray.jl): for example

cat(m,m;dims=(1,2))

does not work while it should.

[dkarrasch]

for example

cat(m,m;dims=(1,2))

does not work while it should.

How should that work? With what would you fill the zero entries?

[jmichel7]

Note that there are zero entries only if m is non-empty. I would fill them with zero(first(m))

[dkarrasch]

The original issue is tricky. For instance, LinearAlgebra.haszero(Mod) returns true, because Mod is a concrete Number type. Moreover, we need to be able to multiply empty matrices, in which case the resulting matrix contains zeros, but there are no elements available to determine the concrete zero element, if not from the type.

[dkarrasch]

Here are a couple of operations that should work and perhaps be added as tests in the course of any attempt of fixing this:

m = Mod.(rand(0:19, 5, 0), 20)
m * copy(m')
m = Mod.([4 5; 5 4], 20)
m * m
m = Mod.(rand(0:19, 3, 3), 20)
m * m
m = Mod.(rand(0:19, 5, 5), 20)
m * m
m*m[:,1]

The 2x2 and 3x3 case are already working on the master branch.

[jishnub]

Perhaps haszero should be more conservative and return false here?

[jmichel7]

we need to be able to multiply empty matrices, in which case the resulting matrix contains zeros, but there are no elements available to determine the concrete zero element, if not from the type

Note that the product of an empty matrix by any matrix is an empty matrix. So returning an empty matrix of the same type is the perfectly right answer. It is a case where there is no need of zero elements.

[jmichel7]

Perhaps haszero should be more conservative and return false here?

A problem is the method in Base/number.jl

zero(::Type{T}) where {T<:Number} = convert(T,0)

which messes the issue. A possibility is to say that modular numbers are not <:Number even so this is somehow a lie.

[dkarrasch]

That's not correct:

julia> randn(5, 0)*randn(0,5)
5×5 Matrix{Float64}:
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0

[jmichel7]

That's not correct:

Yes, you are right. Well, this case should fail (it failed in julia1.10 and previous versions).

[dkarrasch]

Hm, I just checked and this has been working at least since v1.0.