Fix hang in nmod_mpoly_gcd

[fredrik-johansson]

I asked Fable 5 to fix #2581.

Its report is below. I did not ask it to follow up on its suggestion

A production fix should probably also bound the scales_not_found /
eval_point_not_found / eval_gcd_deg_too_high retries in both
gcdp_zippel inner loops (defense in depth — any other systematic failure
mode currently loops ~2^31 times for a 31-bit modulus), and could grow l
more carefully (e.g. +flen per round). The fmpz_mpoly Zippel path reduces to
these nmod/fq_nmod kernels, so both benefit.

since we currently have no other hanging examples, but someone else may want to investigate that.

Investigation of #2581 — nmod_mpoly_gcd hangs

Reproduced on current master (474f55e), p = 2^31 − 1, 12 variables. Baseline behavior:

call	original (g, h)	deflated pair (sg2, sh2)
nmod_mpoly_gcd	hang	—
nmod_mpoly_gcd_zippel	0.9 ms ✓	hang
nmod_mpoly_gcd_zippel2	hang	0.9 ms ✓
nmod_mpoly_gcd_hensel	hang	—
nmod_mpoly_gcd_brown	5.8 s ✓	—

Where the hang is

All of the hanging entry points funnel into the same loop. Stack sample of nmod_mpoly_gcd:

nmod_mpoly_gcd → _try_zippel2 → nmod_mpolyl_content (num_vars=2)
  → _nmod_mpoly_vec_content_mpoly → nmod_mpoly_gcd (recursive, on coefficients)
    → _try_zippel → nmod_mpolyu_gcdm_zippel
      → nmod_mpolyu_gcdp_zippel (var 7 → 6 → 5 → 4 → 3)
        → nmod_mpolyu_gcds_zippel → nmod_mat_rref   ← spinning here

_try_zippel2 (and _try_hensel) begin by removing content in the two main
variables; that content computation is a recursive nmod_mpoly_gcd, which
selects classic Zippel — so gcd, gcd_zippel2 and gcd_hensel all hang in
the same place. The "deflated polynomials" in the earlier comment are exactly
that inner content subproblem.

What the loop is doing

Instrumenting nmod_mpolyu_gcds_zippel (the LINZIP sparse-interpolation step)
shows that at recursion level var = 3 it returns
nmod_gcds_scales_not_found on every single call:

gcds call 7:     var=3 flen=3 ret=scales_not_found
gcds call 8:     var=3 flen=3 ret=scales_not_found
...
gcds call 20000: var=3 flen=3 ret=scales_not_found   (totals: 19994 × scales_not_found)

The caller nmod_mpolyu_gcdp_zippel treats scales_not_found as an unlucky
evaluation point and retries with a new α (goto inner_continue). Its loop
only terminates when α wraps all the way around the field — with p ≈ 2^31 and
an rref per iteration, that is an effective hang (months of CPU).

Root cause: l (number of images) is too small; the failure is structural, not bad luck

Dumping the failing subproblem (4 remaining variables; main variable X of
degree 2; coefficient supports of sizes 4, 8, 12) shows:

• the true gcd there is equal to B (B divides A),
• the assumed form f matches the gcd's support exactly,
• the gcd is primitive (content 1 in the non-main variables).

So neither the form nor content is the problem. Re-implementing the LINZIP
multiple-scaling linear system independently (same Kronecker-power evaluation
scheme, same per-block rref + stacked constraint matrix) shows that for this
support structure the system is generically rank-deficient at the image
count FLINT chooses:

l = 13 (FLINT's choice)  →  nullity = 3     (every random α; tested many)
l = 14                   →  nullity = 2
l = 15                   →  nullity = 1     (unique solution, as required)

FLINT computes l = (Σ len_j + flen − 3)/(flen − 1), max'd with the largest
block, plus one test image — a naive equation count that assumes the stacked
constraint rows are generically independent. For this input they are not: two
images' worth of equations are dependent, so the system stays underdetermined
no matter which evaluation point is chosen. Retrying with new points can
therefore never succeed, but nothing in the code path ever increases l
or gives up:

• nmod_mpolyu_gcds_zippel retries up to 10 fresh points internally, then
  returns scales_not_found;
• nmod_mpolyu_gcdp_zippel responds by picking yet another α and calling
  gcds again with the same form → the same l → the same nullity — for up
  to p − 2 iterations.

Two independent defects

1. No escalation on underdetermination. LINZIP's prescription for an
   underdetermined scale system is to add more images/equations. FLINT
   restarts with a fresh point but an identical system size.

2. Unbounded retry of a non-transient failure. The inner loop in
   gcdp_zippel bounds nothing except α wrapping the whole field. Capping the
   retries alone is not sufficient here: with only a retry cap, the failure
   propagates to nmod_mpolyu_gcdm_zippel's extension-field fallback, and the
   identical LINZIP code in fq_nmod_mpolyu_gcds_zippel hangs the same way
   (verified by stack sampling). The same duplicated logic lives in
   src/nmod_mpoly/mpolyu_gcdp_zippel.c and
   src/fq_nmod_mpoly/mpolyu_gcdp_zippel.c.

Proof-of-concept fix

flint-issue2581-poc-fix.patch modifies both copies of gcds_zippel: when
the system is still underdetermined after the existing two-point retry, free
the l-sized objects, grow l by ~50% (up to 4 times), and re-run — instead of
returning scales_not_found.

With the patch (nothing else changed), every previously hanging case
completes, all with the correct answer l1*a2 − l1*a3 − l2*a1 + l2*a3 + l3*a1 − l3*a2:

call	original (g, h)	deflated pair
nmod_mpoly_gcd	4.7 ms	—
gcd_zippel	1.0 ms	3.8 ms
gcd_zippel2	4.7 ms	1.2 ms
gcd_hensel	4.3 ms	—

Regression tests all pass: nmod_mpoly_gcd, gcd_zippel, gcd_zippel2,
gcd_hensel, gcd_cofactors, gcd_brown, and fq_nmod_mpoly_gcd,
fq_nmod_mpoly_gcd_zippel, fq_nmod_mpoly_gcd_cofactors.

A production fix should probably also bound the scales_not_found /
eval_point_not_found / eval_gcd_deg_too_high retries in both
gcdp_zippel inner loops (defense in depth — any other systematic failure
mode currently loops ~2^31 times for a 31-bit modulus), and could grow l
more carefully (e.g. +flen per round). The fmpz_mpoly Zippel path reduces to
these nmod/fq_nmod kernels, so both benefit.

Regression tests

flint-issue2581-fix-and-tests.patch supersedes the earlier PoC patch: it
contains the fix for both LINZIP copies plus regression tests in all six
affected test files. The test inputs are built compactly from the factored
forms (the deflated pair's 29-term common factor turns out to be exactly the
large cofactor of h, up to variable renaming):

test file	input	exercises
nmod_mpoly/test/t-gcd.c	original (g, h)	default algorithm cascade (zippel2 → content → zippel)
nmod_mpoly/test/t-gcd_cofactors.c	original (g, h)	cofactor variant of the cascade
nmod_mpoly/test/t-gcd_zippel2.c	original (g, h)	zippel2's content preprocessing
nmod_mpoly/test/t-gcd_hensel.c	original (g, h)	hensel's content preprocessing
nmod_mpoly/test/t-gcd_zippel.c	deflated pair	classic Zippel / nmod_mpolyu_gcds_zippel directly
fq_nmod_mpoly/test/t-gcd_zippel.c	deflated pair over F_p (deg-1 ctx)	the fq_nmod twin fq_nmod_mpolyu_gcds_zippel directly

Validation performed both ways:

• Unpatched library: all six tests hang (killed by a 12 s timeout) — they
  reliably catch the bug, including the fq_nmod branch that a fix to only the
  nmod copy would miss.
• Patched library: all six pass in normal test time (0.2–1.1 s per test
  function), along with nmod_mpoly_gcd_brown and fq_nmod_mpoly_gcd as
  broader regression checks.