[WIP] Add implementation of partial sqrt algorithm

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);
 
@@ -934,6 +937,12 @@ FLINT_DLL void nf_elem_rep_mat(fmpq_mat_t res, const nf_elem_t a, const nf_t nf)
 
 FLINT_DLL void nf_elem_rep_mat_fmpz_mat_den(fmpz_mat_t res, fmpz_t den, const nf_elem_t a, const nf_t nf);
 
+FLINT_DLL int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf);
+
+FLINT_DLL int nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf);
+
+FLINT_DLL int nf_elem_is_square(const nf_elem_t b, const nf_t nf);
+
 /******************************************************************************
 
     Modular reduction

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)
 
@@ -329,6 +336,20 @@ void nf_elem_trace(fmpq_t res, const nf_elem_t a, const nf_t nf)
     Set \code{res} to the absolute trace of the given number field
     element $a$.
 
+int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf)
+
+    If \code{a} is a square, set \code{a} to the square root and return
+    \code{1}, else return \code{0}. No aliasing is allowed.
+
+int nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf)
+
+    If \code{a} is a square, set \code{a} to the square root and return                                                     \code{1}, else return \code{0}.
+
+int nf_elem_is_square(const nf_elem_t b, const nf_t nf)
+
+    Return \code{1} if \code{b} is a square. This is currently not faster than
+    computing the square root.
+
 *******************************************************************************
 
     Representation matrix

nf_elem/is_square.c

@@ -0,0 +1,32 @@
+/*=============================================================================
+
+    This file is part of Antic.
+
+    Antic is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version. See <http://www.gnu.org/licenses/>.
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2020 William Hart
+
+******************************************************************************/
+
+#include "nf_elem.h"
+
+int nf_elem_is_square(const nf_elem_t b, const nf_t nf)
+{
+   nf_elem_t t;
+
+   int ret;
+
+   nf_elem_init(t, nf);
+
+   ret = _nf_elem_sqrt(t, b, nf);
+
+   nf_elem_clear(t, nf);
+
+   return ret;
+}

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)
 {