nmod_poly_sqrt: throw error for non-prime modulus

[cxzhong]

The nmod_poly_sqrt, nmod_poly_sqrt_series, and nmod_poly_invsqrt_series functions only work correctly for prime moduli because:

1. n_sqrtmod uses Tonelli-Shanks algorithm (requires prime)
2. gr_sqrt returns GR_UNABLE for non-fields
3. Division by 2 (gr_mul_2exp_si) fails for even moduli
4. Newton iteration requires invertible elements

Previously, using these functions with composite moduli like 16 would cause a crash with an uninformative GR_MUST_SUCCEED failure message.

Now they throw a clear error: 'Modulus must be prime.'

Fixes #2508

[thofma]

This is quite an expensive test. I think the usual flint philosophy is to state in the documentation that the modulus must be prime. If it is not prime, then "garbage in, garbage out" applies.

[cxzhong]

This is quite an expensive test. I think the usual flint philosophy is to state in the documentation that the modulus must be prime. If it is not prime, then "garbage in, garbage out" applies.

Thank you. I will close that. Then I think we need to write this in document. But I also do not see this in document.

[thofma]

Let's wait for some of the maintainers to chime in.

[fredrik-johansson]

As the primality test indeed can be relatively expensive, it would be better to try to compute the square root, and catch when this fails.

If p is not prime, then I think n_sqrtmod(a, p) still works, but either returns 0 or a bogus value. For nmod_poly_sqrt, I think this means you could just go ahead and compute the square root of the constant term and check the result with an nmod_mul. So it might work to change this block

    c = d = p[0];
    if (c != 1)
    {
        c = n_sqrtmod(c, mod.n);
        if (c == 0)
            return 0;
    }

into this (warning: untested):

    c = d = p[0];
    if (c != 1)
    {
        c = n_sqrtmod(c, mod.n);
        if (c == 0)
            return 0;
        if (nmod_mul(c, c, mod) != d)
            flint_throw(FLINT_ERROR, "Exception (nmod_poly_sqrt). "Modulus must be prime.\n");
    }

Then I'd guess you'd also need to intercept failure of the n_invmod below which converts to monic, and you'd need to check for even moduli where _nmod_poly_sqrt_series fails; all of these checks would be much cheaper than a primality test.

Even more user-friendly would be to fix n_sqrtmod to work with any modulus and have an n_sqrtmod_prime for when p is known to be prime.

For the functions wrapping gr algorithms the error handling is already implemented, so here you could just change

GR_MUST_SUCCEED(_gr_poly_sqrt_series(g, h, hlen, n, ctx))

to throw a more informative error, e.g.:

if (_gr_poly_sqrt_series(g, h, hlen, n, ctx) != GR_SUCCESS)
    flint_throw(FLINT_ERROR, "_nmod_poly_sqrt_series: modulus is not prime, or square root does not exist");

(Optionally, do some more checks when the computation actually failed to make the error message more specific.)

Actually, _gr_poly_sqrt_series and _gr_poly_rsqrt_series already end up doing a primality test if the constant term isn't 1, so your PR ends up doing the primality test twice. As noted above, this existing primality test could be optimized away by validating the output of n_sqrtmod.

These functions are really a bit of a mess at the moment; there's lots of things to clean up and optimize...

And as you noted in #2508, it would be nice if nmod_poly_sqrt supported composite moduli (e.g. CRT and Hensel lifting). But I agree that a more informative error message is a good intermediate solution.