Speed up.

nf_elem/sqrt.c

@@ -16,7 +16,7 @@
 
 #include "nf_elem.h"
 
-#define DEBUG 1
+#define DEBUG 0
 
 /*
    TODO:
@@ -99,7 +99,7 @@ int nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf)
       nf_elem_t sqr;
       slong i, j, k;
       fmpz * r, * mr;
-      int res = 0, factored;
+      int res = 0, factored, iters;
 
       if (!fmpz_is_one(NF_ELEM_DENREF(b)))
       {
@@ -150,7 +150,7 @@ int nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf)
 #endif
 
       bbits = FLINT_ABS(_fmpz_vec_max_bits(NF_ELEM_NUMREF(b), lenb));
-      nbits = (bbits + 1)/2 + 2;
+      nbits = (bbits + 1)/(2*lenf) + 2;
 
       /*
          Step 3: find a nbits bit prime such that z = f(n) is a product
@@ -177,7 +177,8 @@ int nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf)
          _fmpz_poly_discriminant(disc, fmpq_poly_numref(nf->pol), lenf);
 
          factored = 0;
-
+         iters = 0;
+
          while (!factored || fac->num > 5) /* no bound known for finding such a factorisation */
          {
             fmpz_factor_clear(fac);
@@ -191,6 +192,14 @@ int nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf)
 
             if (!factored)
                factored = fmpz_is_probabprime(fac->p + fac->num - 1);
+
+            if (nbits < 20 && iters >= (1<<(nbits-1)))
+            {
+               iters = 0;
+               nbits++;
+            }
+
+            iters++;
          }
 
          flint_randclear(state);

nf_elem/test/t-sqrt.c

@@ -30,15 +30,15 @@ main(void)
     flint_randinit(state);
 
     /* test sqrt(a^2) */
-    for (i = 0; i < 1000 * antic_test_multiplier(); )
+    for (i = 0; i < 100 * antic_test_multiplier(); )
     {
         nf_t nf;
         nf_elem_t a, b, c, d;
         int is_square, num_facs;
         fmpz_poly_factor_t fac;
         fmpz_poly_t pol; /* do not clear */
 
-        nf_init_randtest(nf, state, 10, 20);
+        nf_init_randtest(nf, state, 30, 30);
 
         fmpz_poly_factor_init(fac);
 
@@ -61,13 +61,12 @@ main(void)
            nf_elem_init(c, nf);
            nf_elem_init(d, nf);
 
-           nf_elem_randtest(a, state, 20, nf);
+           nf_elem_randtest(a, state, 30, nf);
 
            fmpz_set_ui(NF_ELEM_DENREF(a), 1);
 
-printf("a = "); fmpq_poly_print_pretty(NF_ELEM(a), "x"); printf("\n");
-
            nf_elem_mul(b, a, a, nf);
+
            is_square = nf_elem_sqrt(c, b, nf);
 
            nf_elem_mul(d, c, c, nf);