From 6177c63344319d8349c03aaca849ee233489100f Mon Sep 17 00:00:00 2001 From: William Hart Date: Wed, 16 Sep 2020 16:04:15 +0200 Subject: [PATCH] Small cleanup and deal with aliasing. --- nf_elem.h | 2 + nf_elem/sqrt.c | 76 +++++++++++++++++++++++------------ nf_elem/test/t-mul_div_fmpq.c | 2 +- nf_elem/test/t-sqrt.c | 11 +++-- 4 files changed, 62 insertions(+), 29 deletions(-) diff --git a/nf_elem.h b/nf_elem.h index f3a9ec353..fc8c05335 100644 --- a/nf_elem.h +++ b/nf_elem.h @@ -934,6 +934,8 @@ 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); /****************************************************************************** diff --git a/nf_elem/sqrt.c b/nf_elem/sqrt.c index 5e038bcc0..de447137e 100755 --- a/nf_elem/sqrt.c +++ b/nf_elem/sqrt.c @@ -21,14 +21,13 @@ /* TODO: - * Handle aliasing * try to reuse information from previous failed attempt * improve bounds * add LM bound termination for nonsquare case * add linear and quadratic cases * Tune the number of primes used in trial factoring * Use ECM and larger recombination for very large square roots - * Prove isomorphism to Z/pZ in all cases or exclude primes + * 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) * Deal with number fields of degree 1 and 2 * Deal with primes dividing denominator of norm @@ -121,7 +120,7 @@ slong _fmpz_poly_get_n_adic(fmpz * sqrt, slong len, fmpz_t z, fmpz_t n) return slen; } -int nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) +int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) { if (nf->flag & NF_LINEAR) { @@ -154,12 +153,6 @@ int nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) fmpz * r, * mr, * bz; int res = 0, factored, iters; - if (!fmpz_is_one(fmpq_poly_denref(nf->pol))) - { - flint_printf("Non-monic defining polynomial not supported in sqrt yet.\n"); - flint_abort(); - } - if (lenb == 0) { nf_elem_zero(a, nf); @@ -215,7 +208,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(bz, lenb)); - nbits = (bbits + 1)/(2*lenf) + 2; + nbits = (bbits + 1)/(20) + 2; /* Step 3: find a nbits bit prime such that z = f(n) is a product @@ -225,29 +218,41 @@ int nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) field K defined by f. Then O_K/P_i is isomorphic to Z/pZ via the map s(x) mod P_i -> s(n) mod p. */ -#if DEBUG - flint_printf("Step 3\n"); -#endif fmpz_factor_init(fac); fmpz_init(z); fmpz_init(n); + fmpz_init(disc); + flint_randinit(state); + + _fmpz_poly_discriminant(disc, fmpq_poly_numref(nf->pol), lenf); + do /* continue increasing nbits until square root found */ { - fmpz_init(disc); - flint_randinit(state); - - _fmpz_poly_discriminant(disc, fmpq_poly_numref(nf->pol), lenf); + fmpz_t fac1; + + fmpz_init(fac1); factored = 0; iters = 0; - while (!factored || fac->num > 6) /* no bound known for finding such a factorisation */ +#if DEBUG + flint_printf("Step 3\n"); +#endif + + while (!factored || fac->num > 14) /* no bound known for finding such a factorisation */ { fmpz_factor_clear(fac); fmpz_factor_init(fac); + /* ensure we don't exhaust all primes of the given size */ + if (nbits < 20 && iters >= (1<<(nbits-1))) + { + iters = 0; + nbits++; + } + fmpz_randprime(n, state, nbits, 0); if (fmpz_sgn(n) < 0) fmpz_neg(n, n); @@ -257,19 +262,17 @@ int nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) factored = fmpz_factor_trial(fac, z, 3512); if (!factored) - factored = fmpz_is_probabprime(fac->p + fac->num - 1); - - if (nbits < 20 && iters >= (1<<(nbits-1))) { - iters = 0; - nbits++; + fmpz_set(fac1, fac->p + fac->num - 1); + fac->num--; + + factored = fmpz_factor_smooth(fac, fac1, FLINT_MIN(20, nbits/5 + 1), 0); } iters++; } - flint_randclear(state); - fmpz_clear(disc); + fmpz_clear(fac1); /* Step 4: compute the square roots r_i of z = b(n) mod p_i for each @@ -379,6 +382,9 @@ int nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) } while (res != 1); cleanup: + flint_randclear(state); + fmpz_clear(disc); + fmpz_clear(n); fmpz_clear(temp); fmpz_clear(z); @@ -394,3 +400,23 @@ int nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) } } +int nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) +{ + nf_elem_t t; + + if (a == b) + { + int ret; + + nf_elem_init(t, nf); + + ret = _nf_elem_sqrt(t, b, nf); + nf_elem_swap(t, a, nf); + + nf_elem_clear(t, nf); + + return ret; + } + else + return _nf_elem_sqrt(a, b, nf); +} \ No newline at end of file diff --git a/nf_elem/test/t-mul_div_fmpq.c b/nf_elem/test/t-mul_div_fmpq.c index a90366fe0..8123eec42 100644 --- a/nf_elem/test/t-mul_div_fmpq.c +++ b/nf_elem/test/t-mul_div_fmpq.c @@ -21,7 +21,7 @@ int main(void) { - int i, result; + int i; flint_rand_t state; flint_randinit(state); diff --git a/nf_elem/test/t-sqrt.c b/nf_elem/test/t-sqrt.c index 774fd9fd4..15e332b78 100644 --- a/nf_elem/test/t-sqrt.c +++ b/nf_elem/test/t-sqrt.c @@ -35,10 +35,15 @@ main(void) 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 */ - nf_init_randtest(nf, state, 30, 30); + flen = n_randint(state, 30) + 2; + fbits = n_randint(state, 30) + 1; + abits = n_randint(state, 30) + 1; + + nf_init_randtest(nf, state, flen, fbits); fmpz_poly_factor_init(fac); @@ -60,8 +65,8 @@ main(void) nf_elem_init(b, nf); nf_elem_init(c, nf); nf_elem_init(d, nf); - - nf_elem_randtest(a, state, 30, nf); + + nf_elem_randtest(a, state, abits, nf); nf_elem_mul(b, a, a, nf);