Bug in `nmod_mpoly_factor`

[fredrik-johansson]

As observed in sagemath/sage#41812, FLINT throws an exception when attempting to factor the following (irreducible) polynomial:

sage: R.<x,y,z> = PolynomialRing(GF(3))
sage: P = x^18 + x^12*y^6 + x^6*y^12 + y^18 + x^12*y^4*z^2 - x^10*y^6*z^2 - x^6*y^10*z^2 + x^4*y^12*z^2 + x^12*y^2*z^4 - x^10*y^4*z^4 - x^4*y^10*z^4 + x^2*y^12*z^4 + x^12*z^6 - x^10*y^2*z^6 + x^6*y^6*z^6 - x^2*y^10*z^6 + y^12*z^6 - x^6*y^2*z^10 - x^4*y^4*z^10 - x^2*y^6*z^10 + x^6*z^12 + x^4*y^2*z^12 + x^2*y^4*z^12 + y^6*z^12 + z^18
sage: P.factor()

Reproducer in C:

#include <flint.h>
#include <nmod_mpoly.h>
#include <nmod_mpoly_factor.h>

int main(void)
{
    nmod_mpoly_ctx_t ctx;
    nmod_mpoly_ctx_init(ctx, 3, ORD_DEGLEX, 3);

    const char *vars[] = {"x", "y", "z"};

    nmod_mpoly_t f;
    nmod_mpoly_init(f, ctx);

    const char * poly_str =
        "x^18 + x^12*y^6 + x^6*y^12 + y^18"
        " + x^12*y^4*z^2 - x^10*y^6*z^2 - x^6*y^10*z^2 + x^4*y^12*z^2"
        " + x^12*y^2*z^4 - x^10*y^4*z^4 - x^4*y^10*z^4 + x^2*y^12*z^4"
        " + x^12*z^6 - x^10*y^2*z^6 + x^6*y^6*z^6 - x^2*y^10*z^6 + y^12*z^6"
        " - x^6*y^2*z^10 - x^4*y^4*z^10 - x^2*y^6*z^10"
        " + x^6*z^12 + x^4*y^2*z^12 + x^2*y^4*z^12 + y^6*z^12"
        " + z^18";

    nmod_mpoly_set_str_pretty(f, poly_str, vars, ctx);

    nmod_mpoly_factor_t fac;
    nmod_mpoly_factor_init(fac, ctx);

    if (!nmod_mpoly_factor(fac, f, ctx))
        flint_abort();

    nmod_mpoly_factor_clear(fac, ctx);
    nmod_mpoly_clear(f, ctx);
    nmod_mpoly_ctx_clear(ctx);

    return 0;
}

Tracing what happens leads to the block of code at lines 890 through 909 in

flint/src/nmod_mpoly_factor/n_bpoly_mod_factor_lgprime.c

Line 890 in 633ec99

next_lift_pow = curr_lift_pow + r;

    next_lift_pow = curr_lift_pow + r;
    next_lift_pow = FLINT_MIN(next_lift_pow, 2*curr_lift_pow);

    n_poly_mod_pow(p1, alpha, curr_lift_pow - prev_lift_pow, ctx);
    _hensel_lift_tree(-1, link, v, w, monicB, 2*r-4, prev_alpha_pow, p1, ctx);

    n_poly_mod_pow(p1, alpha, next_lift_pow - curr_lift_pow, ctx);

    n_poly_mod_mul(next_alpha_pow, next_alpha_pow, p1, ctx);
    n_bpoly_set(monicB, B);
    n_bpoly_mod_make_monic_mod(monicB, next_alpha_pow, ctx);

    _hensel_lift_tree(0, link, v, w, monicB, 2*r-4, curr_alpha_pow, p1, ctx);

    prev_lift_pow = curr_lift_pow;
    curr_lift_pow = next_lift_pow;
    n_poly_swap(prev_alpha_pow, curr_alpha_pow);
    n_poly_swap(curr_alpha_pow, next_alpha_pow);

    goto more;

According to codecov this block is not reached by our test code, so it's unsurprising that it has a bug.

The code looks obviously wrong because next_alpha_pow is the zero polynomial when we enter this block, and n_bpoly_mod_make_monic_mod then fails trying to invert modulo the zero polynomial. It's not obvious (to me) how to fix the logic. However, a workaround seems to be to simply replace the code block with:

    goto next_alpha;

I.e. instead of trying to make this alpha work, just try one with higher degree.

Before implementing this solution, I'll post this issue here in case anyone else has a better idea.

[edgarcosta]

Isn't the root cause the typo in:

flint/src/nmod_mpoly_factor/n_bpoly_mod_factor_lgprime.c

Line 898 in 633ec99

n_poly_mod_mul(next_alpha_pow, next_alpha_pow, p1, ctx);

I think it should be:

  n_poly_mod_mul(next_alpha_pow, curr_alpha_pow, p1, ctx);

The intended invariant is next_alpha_pow = alpha^next_lift_pow. But the code multiplies next_alpha_pow by itself, which is the zero polynomial on the first entry to more.

goto next_alpha would avoid the crash but unnecessarily abandons a valid alpha and increases the extension degree.

What do you think?

[fredrik-johansson]

That's logical, but that change by itself leads to a crash.