fmpz_mat.snf randomly hangs

[antieau]

This bug was identified by Claude Code in the course of extensive correctness testing between several SNF algorithms. The report below was produced by Claude Code as well. I have verified that this computation completes in almost no time in SAGE (running M.smith_form(transformation=False)) and seemingly hangs with FLINT. The C file used to check this by Claude Code is attached.

FLINT fmpz_mat.snf() hangs on certain non-square matrices

Summary

flint.fmpz_mat.snf() (python-flint 0.8.0) hangs indefinitely on certain
non-square integer matrices with entries in a moderate range. The same
computation completes in milliseconds with PARI/GP (cypari2).

Environment

• FLINT: 3.4.0 (Homebrew)
• python-flint: 0.8.0
• Python: 3.12.3
• Platform: macOS 14 (Darwin 23.1.0)

Reproducer

import flint

M = [
    [4, -8, 9, -9, 3, 0, 9, -2, -10, -8, -3, 8, 8, -10, -2, 8, -9, -5, 5, 6, 10, 4, -2, -5, 8],
    [3, 10, 5, -8, 5, 1, 3, 0, 0, -7, -5, 0, 3, 5, -1, 2, 7, -9, 4, -8, 0, -2, 0, -7, 2],
    [6, -10, 7, 4, 3, -9, -4, 6, 1, 9, 5, 10, 4, -9, -4, -2, 7, -6, -1, 4, 5, -7, -10, 10, 9],
    [-3, -5, -1, 7, -10, 7, 3, -8, -3, -7, 4, -7, 10, -6, 5, -1, 6, -2, 3, 5, 5, -3, 4, 7, -6],
    [2, -4, 9, 6, -6, -8, -2, 3, 0, 6, -2, -10, -1, -1, 8, 8, 5, -6, 4, 7, 5, 1, 0, 7, 7],
    [2, 4, 0, -4, -3, 8, 2, -3, 3, -9, 0, 5, 2, 2, 10, -6, 5, -9, -6, 6, 8, 0, -7, 4, -7],
    [6, 4, -10, -6, 3, 10, -6, -8, 5, -2, 0, 9, 2, 10, -8, 0, 7, 2, 0, 10, 5, 7, -9, 9, -8],
    [-3, 10, -1, -3, -8, 3, -7, 10, -7, 4, -5, -1, -10, -9, 0, -9, -1, 1, 1, 3, -6, -3, 6, 3, 8],
    [-5, -5, -5, -8, 9, 2, 9, -3, 5, 8, -6, -3, 4, 10, -2, 4, -2, -10, 4, -1, 7, -5, -8, 4, 1],
    [8, -1, 10, 3, -2, 4, -1, -4, 2, 5, -7, -3, 2, 8, 1, 8, -1, -1, -10, 2, -2, -10, 8, -9, 9],
    [5, -1, -3, 9, 1, -3, 10, -4, 9, -2, -6, 10, -7, 10, 10, -9, -1, 4, -9, 8, 1, -6, -8, -1, 0],
    [3, -5, -4, -6, 7, 1, 6, 6, -2, -5, -2, 5, -1, 0, -7, 4, -8, -6, -3, 2, 7, 1, -8, 2, -10],
    [-2, 7, -7, 4, 1, -2, 8, 2, 10, 1, -7, -3, 5, -10, 9, 7, 0, 9, -3, 10, -8, 10, 4, -1, 10],
    [3, -7, -6, -9, -9, -1, 5, -7, -7, -3, 7, -6, 2, 4, 1, 7, 3, 8, -6, 3, 10, -7, 5, 9, 3],
    [-2, -9, 1, -4, 4, 4, -3, 1, -7, 1, 7, 10, 1, -9, 2, -2, -4, -7, 4, -8, -4, 10, 10, 9, -10],
    [-9, 0, -3, -6, 8, -4, -8, 7, -4, 8, -4, -3, 0, -6, 9, -10, -2, -6, -6, 7, -2, -5, -7, -10, -6],
    [-10, 1, -3, 8, 0, -10, -5, -2, -9, -6, 3, 6, -7, -8, 5, 4, 1, 6, 8, -7, 4, 6, -3, 9, -9],
    [6, -1, 4, 10, -10, -9, 5, 2, 3, -7, 5, 4, -8, -8, 0, 9, -6, -8, -6, -2, 9, 10, 8, 7, 0],
    [2, 9, 6, -1, 4, 6, 9, 3, -7, -7, 10, 10, 7, -4, 3, 4, -3, 3, 0, 4, 2, 3, -7, 0, 3],
    [0, -2, 1, -6, 5, -8, -8, -8, -8, 3, -7, 1, -6, 7, -9, 8, 7, 7, 0, -7, 3, 1, 3, -9, -1],
    [9, -1, 1, -7, 8, 6, -4, -6, 5, -3, -7, 1, 7, 1, -7, -2, 8, -3, 3, 7, 9, 9, 10, 7, -10],
    [9, -2, -10, -5, -2, -1, 0, 1, -10, -5, -6, 8, 2, -8, -6, 10, -10, -8, 6, -4, 2, 3, 4, 0, -5],
    [1, -1, 0, 8, 9, -8, -9, -6, -5, 9, -9, -8, -2, 4, 3, 5, 9, 4, 3, -2, -4, 6, -7, 1, 3],
    [-7, -1, 8, 5, 6, -1, -9, -3, 2, 9, -9, -10, -4, -1, -4, -6, -2, -1, 0, -7, -10, 5, 3, -5, -6],
    [2, 7, -3, 6, 7, 1, -8, 2, -9, 3, -10, 4, -8, 0, 8, 3, 8, 2, 10, 3, -1, -7, 2, -10, 0],
    [-5, 9, 4, 1, -8, 3, -7, -3, 3, 8, 2, 6, -8, 2, -1, 0, -3, 0, -5, -8, 6, 10, -7, 6, 6],
    [-4, 1, 1, 10, -6, -3, -7, -6, -2, -4, -5, 9, -6, 10, -8, -5, 10, 5, 4, 8, 8, 4, 8, 10, 10],
]

# This hangs indefinitely:
mat = flint.fmpz_mat(M)
snf = mat.snf()  # <-- never returns

Expected behaviour

snf should return promptly. The correct answer (verified with PARI/GP via
cypari2, and with MAGMA) is the 25×25 identity embedded in a 27×25 matrix — all
25 invariant factors equal 1:

invariant_factors = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]

cypari2 computes this in ~0.007 s.

Observations

• Square matrices are unaffected: fmpz_mat(n, n, ...).snf() completes
  correctly for all tested sizes (up to 30×30).
• Some non-square matrices are slow but finish: e.g. 20×25 can take ~0.2 s.
• The hang is input-dependent, not purely size-dependent: a different 27×25
  matrix may complete; this specific matrix does not.
• fmpz_mat.hnf() is unaffected: the same matrix returns instantly from
  .hnf().

How this was found

Discovered during random cross-backend consistency testing in
tests/test_random_crossbackend.py (snforacle project). The matrix is index 2
from _gen_matrices(1000, seed=42):

import random

def gen_matrices(seed=42):
    rng = random.Random(seed)
    for i in range(1000):
        nrows = rng.randint(10, 30)
        ncols = rng.randint(10, 30)
        entries = [[rng.randint(-10, 10) for _ in range(ncols)] for _ in range(nrows)]
        yield i, nrows, ncols, entries

i, nrows, ncols, entries = list(gen_matrices())[2]
# i=2, nrows=27, ncols=25

Root cause: bug is in FLINT itself, not the Python binding

A minimal C program calling fmpz_mat_snf() directly (no Python, no
python-flint) also hangs:

#include <flint/fmpz_mat.h>

fmpz_mat_t A, S;
fmpz_mat_init(A, 27, 25);
/* ... fill entries ... */
fmpz_mat_init(S, 27, 25);
fmpz_mat_snf(S, A);   /* hangs indefinitely */

Compiled and run against FLINT 3.4.0 (/opt/homebrew/Cellar/flint/3.4.0).
python-flint's fmpz_mat.snf() is just a thin Cython wrapper around this call,
so the binding is not at fault.

The bug should be reported at https://github.com/flintlib/flint/issues (the
FLINT C library), not the python-flint repository.

Not verified against FLINT versions beyond 3.4.0.

flint_snf_test.c