nmod_poly_sqrt: throw error for non-prime modulus

src/nmod_poly/invsqrt_series.c

@@ -10,15 +10,19 @@
     (at your option) any later version.  See <https://www.gnu.org/licenses/>.
 */
 
+#include "flint.h"
 #include "nmod_poly.h"
 #include "gr_poly.h"
 
 void
 _nmod_poly_invsqrt_series(nn_ptr g, nn_srcptr h, slong hlen, slong n, nmod_t mod)
 {
     gr_ctx_t ctx;
+
     _gr_ctx_init_nmod(ctx, &mod);
-    GR_MUST_SUCCEED(_gr_poly_rsqrt_series(g, h, hlen, n, ctx));
+    if (_gr_poly_rsqrt_series(g, h, hlen, n, ctx) != GR_SUCCESS)
+        flint_throw(FLINT_ERROR, "Exception (_nmod_poly_invsqrt_series). "
+            "Modulus is not prime, or square root does not exist.\n");
 }
 
 void

src/nmod_poly/sqrt.c

@@ -9,6 +9,8 @@
     (at your option) any later version.  See <https://www.gnu.org/licenses/>.
 */
 
+#include "flint.h"
+#include "nmod.h"
 #include "ulong_extras.h"
 #include "nmod_vec.h"
 #include "nmod_poly.h"
@@ -60,6 +62,9 @@ _nmod_poly_sqrt(nn_ptr s, nn_srcptr p, slong len, nmod_t mod)
         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");
     }
 
     if (len == 1)
@@ -70,23 +75,35 @@ _nmod_poly_sqrt(nn_ptr s, nn_srcptr p, slong len, nmod_t mod)
 
     slen = len / 2 + 1;
 
+    /* Even modulus: _nmod_poly_sqrt_series requires prime modulus */
+    if (mod.n % 2 == 0)
+        flint_throw(FLINT_ERROR, "Exception (nmod_poly_sqrt). "
+            "Modulus must be prime.\n");
+
     t = _nmod_vec_init(len);
 
     if (c == 1)
         _nmod_poly_sqrt_series(s, p, slen, slen, mod);
     else
     {
-       _nmod_vec_scalar_mul_nmod(t, p, slen, n_invmod(d, mod.n), mod);
+        ulong dinv, g;
+        g = n_gcdinv(&dinv, d, mod.n);
+        if (g != 1)
+            flint_throw(FLINT_ERROR, "Exception (nmod_poly_sqrt). "
+                "Modulus must be prime.\n");
+       _nmod_vec_scalar_mul_nmod(t, p, slen, dinv, mod);
         _nmod_poly_sqrt_series(s, t, slen, slen, mod);
     }
 
     if (c != 1)
         _nmod_vec_scalar_mul_nmod(s, s, slen, c, mod);
 
-    _nmod_poly_mulhigh(t, s, slen, s, slen, slen, mod);
-
+    /* Full verification: s^2 must equal p (all coefficients).
+       This catches incorrect results from Newton iteration
+       when the modulus is not prime. */
+    _nmod_poly_mul(t, s, slen, s, slen, mod);
 
-    result = _nmod_vec_equal(t + slen, p + slen, len - slen);
+    result = _nmod_vec_equal(t, p, len);
     _nmod_vec_clear(t);
     return result;
 }

src/nmod_poly/sqrt_series.c

@@ -10,15 +10,19 @@
     (at your option) any later version.  See <https://www.gnu.org/licenses/>.
 */
 
+#include "flint.h"
 #include "nmod_poly.h"
 #include "gr_poly.h"
 
 void
 _nmod_poly_sqrt_series(nn_ptr g, nn_srcptr h, slong hlen, slong n, nmod_t mod)
 {
     gr_ctx_t ctx;
+
     _gr_ctx_init_nmod(ctx, &mod);
-    GR_MUST_SUCCEED(_gr_poly_sqrt_series(g, h, hlen, n, ctx));
+    if (_gr_poly_sqrt_series(g, h, hlen, n, ctx) != GR_SUCCESS)
+        flint_throw(FLINT_ERROR, "Exception (_nmod_poly_sqrt_series). "
+            "Modulus is not prime, or square root does not exist.\n");
 }
 
 void