Resupdate

[gilvillard]

mbasis "resupdate" and dispatch by shape

nmod_mat_poly.h has long sketched a resupdate variant of mbasis
(residual maintained incrementally instead of recomputed from scratch each
iteration) and a dispatcher choosing between it and the current rescomp
variant by shape.

• nmod_mat_poly_mbasis (the current function) is renamed to
  nmod_mat_poly_mbasis_rescomp, unchanged internally (pure rename, same
  behaviour, same known limitations).
• nmod_mat_poly_mbasis_resupdate is added: maintains the vector of future
  residual coefficients R[d], d = ord..order-1, updating it via the same
  elementary row operations applied to appbas at each iteration, instead of
  recomputing the residual from appbas*matp from scratch.
• nmod_mat_poly_mbasis becomes the dispatcher: 2*cdim > rdim picks
  resupdate, else rescomp.
• Both variants apply identical row operations and the same nullspace pivot
  choice, so they always agree bit-for-bit on the same input — the dispatch
  is a pure timing decision, never a correctness one.
• Also includes a small update-step optimization (_mbasis_low_rank_addmul,
  scoped to the new resupdate function only): the
  per-iteration update bottom += nsbas*top is a rank-cdim update
  generically; below a small threshold (4, tuned in-context) it's done as
  explicit rank-1 updates instead of paying nmod_mat_addmul's per-entry
  verhead for a tiny inner dimension.

Portions developed with the assistance of Claude (Anthropic, Sonnet 5).

2*cdim > rdim vs cdim > rdim/2 + 1

rescomp's residual-recomputation cost is O(m n^2 order^2 / 2); resupdate's
residual-maintenance cost is O(m n (m-n) order^2 / 2) — a factor n/(m-n)
cheaper. The two costs are equal exactly at n = m-n, i.e. cdim = rdim/2.
2*cdim > rdim (strict) is the threshold validated by this data; at exactly
cdim = rdim/2 resupdate has no measured edge, so rescomp is used there.

A bug found (fixed here)

Found that the existing function has a pre-existing latent bug at order = 0:
the main loop never runs in that case, so nsbas is never initialized — but
the unconditional nmod_mat_clear(nsbas) at the end of the function is called
on it anyway.
Fixed in both rescomp and resupdate by giving nsbas a trivial 0 x 0
init before the loop, and changing the loop's own if (ord >= 1) nmod_mat_clear(nsbas) guard to an unconditional clear-before-reassign
(now always valid, since nsbas is always already initialized).

Test added

src/nmod_mat_poly_extra/test/t-mbasis_variants.c: targets
rescomp/resupdate/the dispatcher directly at the mat_poly level —
checks the two variants agree bit-for-bit, the dispatcher picks the right
one, and the result is a genuine approximant basis
(nmod_poly_mat_is_approximant_basis) via conversion. Includes an
explicit (non-random) order=0 regression case for the bug above.
Complements the existing nmod_poly_mat_extra/test/t-mbasis.c, which
exercises the nmod_poly_mat_mbasis wrapper and, since that wrapper calls
into this dispatcher, already provides transitive coverage across both
shape regimes via its own random rdim/cdim sweep.

What needs to be studied further

• The 2*cdim > rdim dispatcher.
• The NMOD_MAT_POLY_MBASIS_LOWRANK_THRES 4 low rank threshold.

Testing performed

• 200 trials specifically for the order = 0 fix, repeated 3x.
• Benchmarks for various input shapes, modulo 2 and modulo a random
  prime with 50 bits, see the files.

Riders

Few things to avoid warnings at compilation when using maple interface. Does
not modify anything if maple is not involved.
In particular, uses now CPPFLAGS="$CPPFLAGS -isystem $withval/include" in
configure.ac.
bench_mod_2.pdf
bench_50_bits.pdf

[vneiger]

Thanks! First comments on the easier parts:

• Existing bug for order 0: in fact approximant bases functions were written assuming that the order is strictly positive (this is somehow indicated in nmod_poly_mat_approximant.h, but this could be more clear). Well, I'm never sure whether supporting order 0 is a good idea (performance aspects versus user-friendliness) so I will open an issue about this to keep it in mind.
• The modifications to configure.ac were triggered by a compilation error, which is due to a missing guard around an include of machine_vectors.h (which should not be included for the detected nehalem architecture). I would rather fix this missing guard rather than do the configure modification, because the problem will still persist for architecture that indeed do not support machine_vectors.h .

[vneiger]

All changes look fine, thanks a lot! And nice speed-ups.

I'll just take time to re-read carefully the new variant, and then this should be ready to merge. (I'd like maybe to do very minor changes here and there but this can be kept for later.)

Minor comment concerning "portions developed with the assistance of Claude (Anthropic, Sonnet 5)": there was some discussion in a FLINT thread about mentioning this in copyrights at the top of the files. For the moment the approach was rather to mention the AI assistance in the pull request here, but leave just your name in the copyright.

[gilvillard]

All changes look fine, thanks a lot! And nice speed-ups.

I'll just take time to re-read carefully the new variant, and then this should be ready to merge. (I'd like maybe to do very minor changes here and there but this can be kept for later.)

Minor comment concerning "portions developed with the assistance of Claude (Anthropic, Sonnet 5)": there was some discussion in a FLINT thread about mentioning this in copyrights at the top of the files. For the moment the approach was rather to mention the AI assistance in the pull request here, but leave just your name in the copyright.

Ok, fine. Modified.

[vneiger]

Thanks!