Smith normal form over residue ring

[taboege]

The problem

The function snf_with_transform returns inconsistent results, at least over residue rings of ZZ. As documented, the returned tuple S,T,U = snf_with_transform(A) should satisfy T*A*U == S but doesn't in the following example:

julia> R, _ = residue_ring(ZZ, 8);

julia> A = R[4 2 6; 4 5 6; 8 8 2];

julia> S,T,U = snf_with_transform(A);

julia> S
[1   0   0]
[0   2   0]
[0   0   0]

julia> T*A*U
[1   0   0]
[0   2   0]
[0   0   4]

The function works fine over ZZ:

julia> S,T,U = snf_with_transform(ZZ.(A));

julia> T*A*U == S
true

julia> S
[1   0    0]
[0   2    0]
[0   0   12]

Note that a Smith normal form over ZZ reduces to a Smith normal form over ZZ/8 by reducing all the entries mod 8. This shows that the matrix S computed over ZZ/8 was wrong (its 3rd diagonal entry should be 4 instead of 0) and T*A*U was right.

Version information

julia> Oscar.versioninfo(full=true)
OSCAR version 1.2.0
  combining:
    AbstractAlgebra.jl   v0.43.9
    GAP.jl               v0.12.0
    Hecke.jl             v0.34.6
    Nemo.jl              v0.47.3
    Polymake.jl          v0.11.22
    Singular.jl          v0.23.10
  building on:
    FLINT_jll               v300.100.300+0
    GAP_jll                 v400.1300.102+2
    Singular_jll            v404.0.606+0
    libpolymake_julia_jll   v0.13.0+0
    libsingular_julia_jll   v0.46.0+0
    polymake_jll            v400.1300.2+0
See `]st -m` for a full list of dependencies.

Julia Version 1.11.1
Commit 8f5b7ca12ad (2024-10-16 10:53 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 96 × Intel(R) Xeon(R) Platinum 8168 CPU @ 2.70GHz
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, skylake-avx512)
Threads: 1 default, 0 interactive, 1 GC (on 96 virtual cores)
Official https://julialang.org/ release

[fingolfin]

Interesting. snf_with_transform is actually a function provided by AbstractAlgebra, but it is underdocumented and I don't think it was written with zero divisors in mind.

But your code snippet actually doesn't run in AA, nor in Nemo (entering it will throw various "not implemented" errors). It only runs in Hecke (and hence in OSCAR) because Hecke actually implements Base.divrem(n::T, m::T) where T <: Union{zzModRingElem,Nemo.ZZModRingElem}.

Anyway, Nemocas/AbstractAlgebra.jl#1899 fixes this particular instances. But I'd not be surprised if there were other issues with it.

[taboege]

Thanks, this is very interesting! I had tracked the SNF code to AbstractAlgebra but couldn't get my example to run, for the reason you mentioned. It seems like there is still something to be done in the implementation over in AbstractAlgebra to make it work in the presence of zero divisors. I'm not sure if this issue should be closed while the one in AA remains open?