Test non-monic case.

nf/init_randtest.c

@@ -22,6 +22,7 @@ void nf_init_randtest(nf_t nf, flint_rand_t state,
 {
     fmpq_poly_t pol;
     fmpz_poly_t q;
+    fmpz_t d;
 
     if (len < 2 || bits_in < 1)
     {
@@ -53,7 +54,15 @@ void nf_init_randtest(nf_t nf, flint_rand_t state,
         fmpz_randtest_not_zero(fmpq_poly_denref(pol), state, bits_in);
     fmpq_poly_canonicalise(pol);
 
+    fmpz_init(d);
+
+    _fmpz_vec_content(d, fmpq_poly_numref(pol), fmpq_poly_length(pol));
+
+    if (!fmpz_is_one(d))
+        _fmpz_vec_scalar_divexact_fmpz(fmpq_poly_numref(pol), fmpq_poly_numref(pol), fmpq_poly_length(pol), d);
+
     nf_init(nf, pol);
     fmpq_poly_clear(pol);
     fmpz_poly_clear(q);
+    fmpz_clear(d);
 }

nf_elem.h

@@ -83,6 +83,9 @@ FLINT_DLL void nf_elem_clear(nf_elem_t a, const nf_t nf);
 FLINT_DLL void nf_elem_randtest(nf_elem_t a, flint_rand_t state, 
                                               mp_bitcnt_t bits, const nf_t nf);
 
+FLINT_DLL void nf_elem_randtest_bounded(nf_elem_t a, flint_rand_t state,
+                                              mp_bitcnt_t bits, const nf_t nf);
+
 FLINT_DLL void nf_elem_randtest_not_zero(nf_elem_t a, flint_rand_t state, 
                                               mp_bitcnt_t bits, const nf_t nf);
 

nf_elem/doc/nf_elem.txt

@@ -34,7 +34,14 @@ void nf_elem_randtest(nf_elem_t a, flint_rand_t state,
                                                mp_bitcnt_t bits, const nf_t nf)
 
     Generate a random number field element $a$ in the number field \code{nf}
-    whose coefficients have up to the given number of bits.
+    whose rational coefficients have up to the given number of bits.
+
+void nf_elem_randtest_bounded(nf_elem_t a, flint_rand_t state,
+                                               mp_bitcnt_t bits, const nf_t nf)
+
+    Generate a random number field element $a$ in the number field \code{nf}
+    whose coefficients have up to the given number of bits, with combined
+    denominator also bounded by the given number of bits.
 
 void nf_elem_canonicalise(nf_elem_t a, nf_t nf)
 

nf_elem/randtest.c

@@ -66,6 +66,84 @@ void nf_elem_randtest(nf_elem_t a, flint_rand_t state,
     }
 }
 
+void nf_elem_randtest_bounded(nf_elem_t a, flint_rand_t state, 
+                                               mp_bitcnt_t bits, const nf_t nf)
+{
+    if (nf->flag & NF_LINEAR)
+    {
+        fmpz_randtest(LNF_ELEM_NUMREF(a), state, bits);
+
+        if (n_randint(state, 2))
+        {
+           fmpz_randtest_not_zero(LNF_ELEM_DENREF(a), state, bits);
+           fmpz_abs(LNF_ELEM_DENREF(a), LNF_ELEM_DENREF(a));
+
+           _fmpq_canonicalise(LNF_ELEM_NUMREF(a), LNF_ELEM_DENREF(a));
+        } else
+           fmpz_one(LNF_ELEM_DENREF(a));
+    } else if (nf->flag & NF_QUADRATIC)
+    {
+        fmpz_randtest(QNF_ELEM_NUMREF(a), state, bits);
+        fmpz_randtest(QNF_ELEM_NUMREF(a) + 1, state, bits);
+
+        if (n_randint(state, 2))
+        {
+           fmpz_t d;
+
+           fmpz_randtest_not_zero(QNF_ELEM_DENREF(a), state, bits);
+           fmpz_abs(QNF_ELEM_DENREF(a), QNF_ELEM_DENREF(a));
+
+           fmpz_init(d);
+           fmpz_gcd(d, QNF_ELEM_NUMREF(a), QNF_ELEM_NUMREF(a) + 1);
+           if (!fmpz_is_one(d))
+           {
+              fmpz_gcd(d, d, QNF_ELEM_DENREF(a));
+
+              if (!fmpz_is_one(d))
+              {
+                 _fmpz_vec_scalar_divexact_fmpz(QNF_ELEM_NUMREF(a), QNF_ELEM_NUMREF(a), 2, d);
+                 fmpz_divexact(QNF_ELEM_DENREF(a), QNF_ELEM_DENREF(a), d);
+              }
+           }
+        } else
+           fmpz_one(QNF_ELEM_DENREF(a));
+    }
+    else
+    {
+        slong i, lenf = nf->pol->length;
+        fmpz_t d;
+
+        for (i = 0; i < lenf - 1; i++)
+            fmpz_randtest(NF_ELEM_NUMREF(a) + i, state, bits);
+
+        if (n_randint(state, 2))                                                                                                {
+            fmpz_init(d);
+
+            fmpz_randtest_not_zero(NF_ELEM_DENREF(a), state, bits);
+            fmpz_abs(NF_ELEM_DENREF(a), NF_ELEM_DENREF(a));
+
+            _fmpz_vec_content(d, NF_ELEM_NUMREF(a), lenf - 1);
+
+            if (!fmpz_is_one(d))
+            {
+                fmpz_gcd(d, d, NF_ELEM_DENREF(a));
+
+                if (!fmpz_is_one(d))
+                {
+                    _fmpz_vec_scalar_divexact_fmpz(NF_ELEM_NUMREF(a), NF_ELEM_NUMREF(a), lenf - 1, d);
+                    fmpz_divexact(NF_ELEM_DENREF(a), NF_ELEM_DENREF(a), d);
+                }
+            }
+
+            fmpz_clear(d);
+        } else
+            fmpz_one(NF_ELEM_DENREF(a));
+
+        _fmpq_poly_set_length(NF_ELEM(a), lenf - 1);
+        _fmpq_poly_normalise(NF_ELEM(a));
+    }
+}
+
 void nf_elem_randtest_not_zero(nf_elem_t a, flint_rand_t state, 
                                                mp_bitcnt_t bits, const nf_t nf)
 {

nf_elem/sqrt.c

@@ -26,7 +26,6 @@
      * add termination bound for nonsquare case
      * Prove homomorphism to Z/pZ in all cases or exclude primes
      * Deal with lousy starting bounds (they are too optimistic if f is not monic)
-     * Fix bug in antic norm function
      * Remove small squares from denominator before rationalising
      * Cache factorisation of f(n) on number field for future square roots
      * add is_square function

nf_elem/test/t-sqrt.c

@@ -240,7 +240,7 @@ main(void)
         fmpq_poly_clear(f);
     }
 
-   /* test sqrt(a^2) */
+    /* test sqrt(a^2) monic defining poly */
     for (i = 0; i < 100 * antic_test_multiplier(); )
     {
         nf_t nf;
@@ -271,12 +271,78 @@ main(void)
         if (nf->flag & NF_MONIC && num_facs == 1)
         {
            i++;
+
+           nf_elem_init(a, nf);
+           nf_elem_init(b, nf);
+           nf_elem_init(c, nf);
+           nf_elem_init(d, nf);
+
+           nf_elem_randtest_bounded(a, state, abits, nf);
+
+           nf_elem_mul(b, a, a, nf);
+
+           is_square = nf_elem_sqrt(c, b, nf);
+
+           nf_elem_mul(d, c, c, nf);
+
+           result = is_square && nf_elem_equal(d, b, nf);
+           if (!result)
+           {
+              printf("FAIL:\n");
+              printf("a = "); nf_elem_print_pretty(a, nf, "x"); printf("\n");
+              printf("b = "); nf_elem_print_pretty(b, nf, "x"); printf("\n");
+              printf("c = "); nf_elem_print_pretty(c, nf, "x"); printf("\n");
+              printf("d = "); nf_elem_print_pretty(d, nf, "x"); printf("\n");
+              abort();
+           }
+
+           nf_elem_clear(a, nf);
+           nf_elem_clear(b, nf);
+           nf_elem_clear(c, nf);
+           nf_elem_clear(d, nf);
+        }
+
+        nf_clear(nf);
+    }
+
+     /* test sqrt(a^2) non-monic defining poly */
+    for (i = 0; i < 100 * antic_test_multiplier(); )
+    {
+        nf_t nf;
+        nf_elem_t a, b, c, d;
+        int is_square, num_facs;
+        slong flen, fbits, abits;
+        fmpz_poly_factor_t fac;
+        fmpz_poly_t pol; /* do not clear */
+
+        flen = n_randint(state, 10) + 2;
+        fbits = n_randint(state, 10) + 1;
+        abits = n_randint(state, 10) + 1;
+
+        nf_init_randtest(nf, state, flen, fbits);
+
+        fmpz_poly_factor_init(fac);
+
+        pol->coeffs = nf->pol->coeffs;
+        pol->length = nf->pol->length;
+        pol->alloc = nf->pol->alloc;
+
+        fmpz_poly_factor(fac, pol);
+
+        num_facs = fac->num*fac->exp[0];
+
+        fmpz_poly_factor_clear(fac);
+
+        if (num_facs == 1)
+        {
+           i++;
+
            nf_elem_init(a, nf);
            nf_elem_init(b, nf);
            nf_elem_init(c, nf);
            nf_elem_init(d, nf);
 
-           nf_elem_randtest(a, state, abits, nf);
+           nf_elem_randtest_bounded(a, state, abits, nf);
 
            nf_elem_mul(b, a, a, nf);