SNF improvements

[edgarcosta]

Building on #2596, this PR adds fmpz_mat_snf_transform (the feature request in flint2#204), fmpz_mat_elementary_divisors, and replaces the SNF dispatcher with a faster algorithm.

fmpz_mat_snf_transform(S, U, V, A) computes S = UAV where S is the Smith normal form and U, V are unimodular, for arbitrary m x n matrices. Either U or V (or both) may be NULL to skip computing the corresponding transform. The algorithm alternates row and column HNF until the matrix is diagonal, then fixes the divisibility chain with xgcd-based row/column operations (O(n) per pair rather than the O(n^3) matrix multiplications in the original acb_theta implementation by Jean Kieffer). This also removes ~100 lines of static duplicates from acb_theta/sp2gz_decompose.c.

fmpz_mat_elementary_divisors(ed, rank, A) computes d_1 | d_2 | ... | d_r without full SNF using Luebeck's algorithm: compute HNF, factor the product of the pivots, then for each prime p iteratively compute nullspace mod p to determine p-adic valuations. Falls back to full SNF if any prime factor does not fit in a ulong. This uses the new nmod_mat_left_nullspace, which computes the left nullspace directly by row-reducing [A | I] rather than transposing.

fmpz_mat_is_diagonal is promoted from a static function in acb_theta/sp2gz_decompose.c (which had a /* todo: move out? */ comment).

The fmpz_mat_snf dispatcher now uses the iterative Hermite algorithm instead of Kannan-Bachem/Iliopoulos, avoiding the cost of determinant computation and Iliopoulos' modular reduction. There is also a fast path for already-diagonal matrices (44-360x speedup).

Timings (us) for square nonsingular matrices, comparing the new dispatcher (snf), snf_transform without and with U/V, Kannan-Bachem (KB), Iliopoulos (Ilio), and elementary_divisors (ed). KB/Ilio skipped for n > 20, ed skipped for 200-bit entries where factorization dominates:

  n  bits     snf  transf  transf+UV        KB   Ilio      ed
  3    10       0       0          0        20      0       0
  3    50       0       0         20         0      0       0
  5    10      20       0          0         0     20       0
  5    50      20      20         40        20      0      40
  8    10       0      20          0        20      0       0
  8    50      20      40         40        60     40      60
  8   200      40      60         80       300     80       -
 10    10       0      20         20        60     20      40
 10    50      20      40         80      1260    100    5540
 10   200      80     160        220      1400    500       -
 15    10     100     100        100     28200    100     100
 15    50     100     200        200     75000    300    4500
 15   200     400     800       1100    680800   3500       -
 20    10     100     200        300   2049600    300     300
 20    50     200     400        500  20440800   1000    1300
 20   200    1200    2700       3100  37894700  10200       -
 30    10     333     667       1000         -      -    1000
 30    50     667    2000       2667         -      -    7333
 30   200    5000   14333      17667         -      -       -
 50    10   14000    5000       7000         -      -   10000
 50    50   16000   18000      23000         -      -   33000
 50   200   88000  142000     235000         -      -       -

The new dispatcher is 5-37000x faster than Kannan-Bachem and 3-9x faster than Iliopoulos at sizes 10-20. Passing NULL for U/V saves 15-40% at larger sizes. elementary_divisors is competitive for small-bit entries but depends on factorization.

Diagonal fast path (us):

   n  bits   new    old  speedup
  20    50     2     88      44x
  50   200    20   4220     211x
 100   200   100  35950     360x

Tests include Magma-verified hardcoded cases and randomized cross-checks between all five SNF implementations (dispatcher, Kannan-Bachem, Iliopoulos, snf_transform, elementary_divisors). Documentation added for all new functions.

[fredrik-johansson]

nmod_mat_left_nullspace seems to compute a reduced matrix of column vectors such that X^T A = 0, whereas nmod_mat_nullspace computes a non-reduced basis of column vectors such that A X = 0.

Wouldn't be more natural to return a matrix of row vectors such that X A = 0, since row matrices in any case are more convenient with our row-major storage? And why not return it non-reduced for users who don't need it to be reduced?

Also, this function can be implemented in terms of nmod_mat_nullspace plus transposition (and nmod_mat_rref if you want it reduced). This should be faster than augmenting the input matrix.

Maybe we should have a function nmod_mat_kernel that allows a combination of flags:

• left/right kernel
• transpose/not transpose (or: basis as rows/column vectors)
• reduced/nonreduced basis
• output zero-padded n x n matrix, or resize matrix to n x r / r x n

Not that this needs to be done in this PR, but it would at least be nice if nullspace/left_nullspace were consistent.

[fredrik-johansson]

rank_out could be the return value instead of being passed as a pointer?

[fredrik-johansson]

I guess the 32-bit test failure is due to overflow in the hardcoded test case slong constant 3395595701880.

[edgarcosta]

Thank you! I appreciate the feedback!
I will hopefully tackle all this soon.

Indeed, I first had the left nulsspace via transpose, but without thinking I thought I was being lazy, and rewrote it via the augmentation... but I can see its issues now.

[fredrik-johansson]

In fmpz_mat_elementary_divisors, why do you multiply together all the leading entries and call fmpz_factor on the big product? Surely it's better to factor each leading entry and merge the factors.

Also, since the factorisation isn't used when a prime doesn't fit in ulong, you could 1) abort early when you find a leading entry with a larger factor, 2) use fmpz_factor_smooth to avoid effectively hanging if there happen to be 1000-bit entries (and/or maybe just abort and fall back on SNF when you hit a number to factor that exceeds some bit limit, say 2*FLINT_BITS).

[edgarcosta]

the commits got a bit out of hand, but to avoid rewriting history right now, I let them be.
At the end if you would like more than one commit, I can squash many of them.