Matrices known to break `hnfc!()`

[brainandforce]

Since I wasn't the one who implemented this algorithm, this might take a while to figure out.

The only example I have as of now is

[ 9  7 -7  7 -3  7
 -1  5  0  7 -4 -1
  2  2  9  7 -5  3
 -2  6  3 -9 -3  5]

which, when passed to hnfc(), fails because of RowHermite requiring that the U matrix be unimodular (the calculated determinant is 0, but it should be ±1).

[brainandforce]

Here's what happens when I run hnfc() and snf() on this matrix:

julia> test
4×6 Matrix{Int64}:
  9  7  -7   7  -3   7
 -1  5   0   7  -4  -1
  2  2   9   7  -5   3
 -2  6   3  -9  -3   5

julia> NormalForms.hnfc(test)
┌ Debug: Matrix state after zeroing along row or column k = 1:
│ A = 4×6 Matrix{Int64}:
│   1   0    0     0   0     0
│  23  52  161  -154  65  -162
│   2   4   23    -7   1   -11
│  30  68  213  -219  87  -205
│ U = 4×6 Matrix{Int64}:
│   1   0    0     0   0     0
│  23  52  161  -154  65  -162
│   2   4   23    -7   1   -11
│  30  68  213  -219  87  -205
└ @ NormalForms ~/git/NormalForms.jl/src/hermite.jl:145
┌ Debug: Matrix state after reducing off-diagonal at k = 1:
│ A = 4×6 Matrix{Int64}:
│   1   0    0     0   0     0
│  23  52  161  -154  65  -162
│   2   4   23    -7   1   -11
│  30  68  213  -219  87  -205
│ U = 4×6 Matrix{Int64}:
│   1   0    0     0   0     0
│  23  52  161  -154  65  -162
│   2   4   23    -7   1   -11
│  30  68  213  -219  87  -205
└ @ NormalForms ~/git/NormalForms.jl/src/hermite.jl:153
┌ Debug: Matrix state after zeroing along row or column k = 2:
│ A = 4×6 Matrix{Int64}:
│   1    0    0      0       0      0
│  23    1    0      0       0      0
│   2  223  552  34335  -14494  36115
│  30   53  128   7943   -3358   8381
│ U = 4×6 Matrix{Int64}:
│   1    0    0      0       0      0
│  23    1    0      0       0      0
│   2  223  552  34335  -14494  36115
│  30   53  128   7943   -3358   8381
└ @ NormalForms ~/git/NormalForms.jl/src/hermite.jl:145
┌ Debug: Matrix state after reducing off-diagonal at k = 2:
│ A = 4×6 Matrix{Int64}:
│      1    0    0      0       0      0
│      0    1    0      0       0      0
│  -5127  223  552  34335  -14494  36115
│  -1189   53  128   7943   -3358   8381
│ U = 4×6 Matrix{Int64}:
│      1    0    0      0       0      0
│      0    1    0      0       0      0
│  -5127  223  552  34335  -14494  36115
│  -1189   53  128   7943   -3358   8381
└ @ NormalForms ~/git/NormalForms.jl/src/hermite.jl:153
┌ Debug: Matrix state after zeroing along row or column k = 3:
│ A = 4×6 Matrix{Int64}:
│      1    0       0      0         0             0
│      0    1       0      0         0             0
│  -5127  223       1      0         0             0
│  -1189   53  452641  -3448  -1358016  -16347121334
│ U = 4×6 Matrix{Int64}:
│      1    0       0      0         0             0
│      0    1       0      0         0             0
│  -5127  223       1      0         0             0
│  -1189   53  452641  -3448  -1358016  -16347121334
└ @ NormalForms ~/git/NormalForms.jl/src/hermite.jl:145
┌ Debug: Matrix state after reducing off-diagonal at k = 3:
│ A = 4×6 Matrix{Int64}:
│           1           0       0      0         0             0
│           0           1       0      0         0             0
│           0           0       1      0         0             0
│  2320689218  -100938890  452641  -3448  -1358016  -16347121334
│ U = 4×6 Matrix{Int64}:
│           1           0       0      0         0             0
│           0           1       0      0         0             0
│           0           0       1      0         0             0
│  2320689218  -100938890  452641  -3448  -1358016  -16347121334
└ @ NormalForms ~/git/NormalForms.jl/src/hermite.jl:153
┌ Debug: Matrix state after zeroing along row or column k = 4:
│ A = 4×6 Matrix{Int64}:
│           1           0       0  0  0  0
│           0           1       0  0  0  0
│           0           0       1  0  0  0
│  2320689218  -100938890  452641  2  0  0
│ U = 4×6 Matrix{Int64}:
│           1           0       0  0  0  0
│           0           1       0  0  0  0
│           0           0       1  0  0  0
│  2320689218  -100938890  452641  2  0  0
└ @ NormalForms ~/git/NormalForms.jl/src/hermite.jl:145
ERROR: AssertionError: isunimodular(U)
Stacktrace:
 [1] hnf_ma!(A::Matrix{Int64}, R::RoundingMode{:Up})
   @ NormalForms ~/git/NormalForms.jl/src/hermite.jl:150
 [2] hnfc! (repeats 2 times)
   @ ~/git/NormalForms.jl/src/hermite.jl:159 [inlined]
 [3] hnfc(M::Matrix{Int64})
   @ NormalForms ~/git/NormalForms.jl/src/hermite.jl:182
 [4] top-level scope
   @ REPL[38]:1

julia> NormalForms.snf(test)
┌ Debug: Matrix state after zeroing row and column k = 1:
│ A = 4×6 Matrix{Int64}:
│  1   0    0     0   0     0
│  0  52  161  -154  65  -162
│  0   4   23    -7   1   -11
│  0  68  213  -219  87  -205
│ U = 4×6 Matrix{Int64}:
│  1   0    0     0   0     0
│  0  52  161  -154  65  -162
│  0   4   23    -7   1   -11
│  0  68  213  -219  87  -205
└ @ NormalForms ~/git/NormalForms.jl/src/smith.jl:106
┌ Debug: Matrix state after enforcing divisibility at k = 1:
│ A = 4×6 Matrix{Int64}:
│  1   0    0     0   0     0
│  0  52  161  -154  65  -162
│  0   4   23    -7   1   -11
│  0  68  213  -219  87  -205
│ U = 4×6 Matrix{Int64}:
│  1   0    0     0   0     0
│  0  52  161  -154  65  -162
│  0   4   23    -7   1   -11
│  0  68  213  -219  87  -205
└ @ NormalForms ~/git/NormalForms.jl/src/smith.jl:115
┌ Debug: Matrix state after zeroing row and column k = 2:
│ A = 4×6 Matrix{Int64}:
│  1  0    0      0       0      0
│  0  1    0      0       0      0
│  0  0  552  34335  -14494  36115
│  0  0  128   7943   -3358   8381
│ U = 4×6 Matrix{Int64}:
│  1  0    0      0       0      0
│  0  1    0      0       0      0
│  0  0  552  34335  -14494  36115
│  0  0  128   7943   -3358   8381
└ @ NormalForms ~/git/NormalForms.jl/src/smith.jl:106
┌ Debug: Matrix state after enforcing divisibility at k = 2:
│ A = 4×6 Matrix{Int64}:
│  1  0    0      0       0      0
│  0  1    0      0       0      0
│  0  0  552  34335  -14494  36115
│  0  0  128   7943   -3358   8381
│ U = 4×6 Matrix{Int64}:
│  1  0    0      0       0      0
│  0  1    0      0       0      0
│  0  0  552  34335  -14494  36115
│  0  0  128   7943   -3358   8381
└ @ NormalForms ~/git/NormalForms.jl/src/smith.jl:115
┌ Debug: Matrix state after zeroing row and column k = 3:
│ A = 4×6 Matrix{Int64}:
│  1  0  0      0         0             0
│  0  1  0      0         0             0
│  0  0  1      0         0             0
│  0  0  0  -3448  -1358016  -16347121334
│ U = 4×6 Matrix{Int64}:
│  1  0  0      0         0             0
│  0  1  0      0         0             0
│  0  0  1      0         0             0
│  0  0  0  -3448  -1358016  -16347121334
└ @ NormalForms ~/git/NormalForms.jl/src/smith.jl:106
┌ Debug: Matrix state after enforcing divisibility at k = 3:
│ A = 4×6 Matrix{Int64}:
│  1  0  0      0         0             0
│  0  1  0      0         0             0
│  0  0  1      0         0             0
│  0  0  0  -3448  -1358016  -16347121334
│ U = 4×6 Matrix{Int64}:
│  1  0  0      0         0             0
│  0  1  0      0         0             0
│  0  0  1      0         0             0
│  0  0  0  -3448  -1358016  -16347121334
└ @ NormalForms ~/git/NormalForms.jl/src/smith.jl:115
┌ Debug: Matrix state after zeroing row and column k = 4:
│ A = 4×6 Matrix{Int64}:
│  1  0  0  0  0  0
│  0  1  0  0  0  0
│  0  0  1  0  0  0
│  0  0  0  2  0  0
│ U = 4×6 Matrix{Int64}:
│  1  0  0  0  0  0
│  0  1  0  0  0  0
│  0  0  1  0  0  0
│  0  0  0  2  0  0
└ @ NormalForms ~/git/NormalForms.jl/src/smith.jl:106
ERROR: AssertionError: isunimodular(V)
Stacktrace:
 [1] snf_ma!(A::Matrix{Int64})
   @ NormalForms ~/git/NormalForms.jl/src/smith.jl:112
 [2] snf!
   @ ~/git/NormalForms.jl/src/smith.jl:117 [inlined]
 [3] snf(M::Matrix{Int64})
   @ NormalForms ~/git/NormalForms.jl/src/smith.jl:138
 [4] top-level scope
   @ REPL[39]:1

julia>

My suspicion is that the matrix entries are overflowing/underflowing. However, I do find it interesting that the matrix determinants end up being zero.

[brainandforce]

The fix/hnfc-failure branch has the debug statements and assertions used to make this (commit dfe8709)

[arnedecadt]

@brainandforce this is caused by integer overflow. It works if you use

test=Int128[ 9  7 -7  7 -3  7
        -1  5  0  7 -4 -1
         2  2  9  7 -5  3
        -2  6  3 -9 -3  5]

I guess an easy fix would be to redo the calculations after casting to BigInt if U is not unimodular.

[brainandforce]

I should have thought about that already...I'll update the tests to check that this works when using a wider integer type.