From d96a40f11eb0f1710819441518c9303ae393a99c Mon Sep 17 00:00:00 2001 From: William Hart Date: Tue, 25 Aug 2020 21:31:50 +0200 Subject: [PATCH 01/24] Add implementation of partial sqrt algorithm. --- nf_elem.h | 2 + nf_elem/sqrt.c | 303 ++++++++++++++++++++++++++++++++++++++++++ nf_elem/test/t-sqrt.c | 99 ++++++++++++++ 3 files changed, 404 insertions(+) create mode 100755 nf_elem/sqrt.c create mode 100644 nf_elem/test/t-sqrt.c diff --git a/nf_elem.h b/nf_elem.h index e9b7c0db1..f3a9ec353 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); + /****************************************************************************** Modular reduction diff --git a/nf_elem/sqrt.c b/nf_elem/sqrt.c new file mode 100755 index 000000000..a443d173e --- /dev/null +++ b/nf_elem/sqrt.c @@ -0,0 +1,303 @@ +/*============================================================================= + + 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 . + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2020 William Hart + +******************************************************************************/ + +#include "nf_elem.h" + +#define DEBUG 1 + +/* + TODO: + + * Handle aliasing + * try to reuse information from previous failed attempt + * improve bounds + * add LM bound termination for nonsquare case + * deal with algebraic integers with denominators + * add linear and quadratic cases +*/ + +slong _fmpz_poly_get_n_adic(fmpz * sqrt, slong len, fmpz_t z, fmpz_t n) +{ + fmpz_t n2, z2; + slong i, slen = 0; + + fmpz_init(n2); + fmpz_init(z2); + + fmpz_fdiv_q_2exp(n2, n, 1); + + if (fmpz_sgn(z) < 0) + fmpz_neg(z2, z); + else + fmpz_set(z2, z); + + for (i = 0; i < len; i++) + { + fmpz_mod(sqrt + i, z2, n); + + fmpz_sub(z2, z2, sqrt + i); + + if (fmpz_cmpabs(sqrt + i, n2) > 0) + { + fmpz_sub(sqrt + i, sqrt + i, n); + fmpz_add(z2, z2, n); + } + + fmpz_fdiv_q(z2, z2, n); + + if (!fmpz_is_zero(sqrt + i)) + slen = i + 1; + } + + fmpz_clear(n2); + fmpz_clear(z2); + + return slen; +} + +int nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) +{ + if (nf->flag & NF_LINEAR) + { + /* const fmpz * const bnum = LNF_ELEM_NUMREF(b); + const fmpz * const bden = LNF_ELEM_DENREF(b); + fmpz * const anum = LNF_ELEM_NUMREF(a); + fmpz * const aden = LNF_ELEM_DENREF(a); */ + + flint_printf("Sqrt for linear number fields not implemented yet\n"); + flint_abort(); + } else if (nf->flag & NF_QUADRATIC) + { + /* const fmpz * const bnum = QNF_ELEM_NUMREF(b); + const fmpz * const bden = QNF_ELEM_DENREF(b); + fmpz * const anum = QNF_ELEM_NUMREF(a); + fmpz * const aden = QNF_ELEM_DENREF(a); */ + + flint_printf("Sqrt for quadratic number fields not implemented yet\n"); + flint_abort(); + } else /* generic nf_elem */ + { + const slong lenb = NF_ELEM(b)->length, lenf = fmpq_poly_length(nf->pol); + slong nbits, bbits; + fmpq_t bnorm; + fmpz_factor_t fac; + flint_rand_t state; + fmpz_t disc, z, temp, n, m; + nf_elem_t sqr; + slong i, j, k; + fmpz * r, * mr; + int res = 0, factored; + + if (!fmpz_is_one(NF_ELEM_DENREF(b))) + { + flint_printf("Sqrt1 outside Z[alpha] not implemented yet\n"); + flint_abort(); + } + + if (lenb == 0) + { + nf_elem_zero(a, nf); + return 1; + } + + /* Step 1: compute norm and check it is square */ +#if DEBUG + flint_printf("Step 1\n"); +#endif + + fmpq_init(bnorm); + + nf_elem_norm(bnorm, b, nf); + + if (!fmpz_is_one(fmpq_denref(bnorm))) + { + flint_printf("Sqrt2 outside Z[alpha] not yet implemented yet\n"); + fmpq_clear(bnorm); + flint_abort(); + } + + if (!fmpz_is_square(fmpq_numref(bnorm))) + { + nf_elem_zero(a, nf); + fmpq_clear(bnorm); + return 0; + } + + /* + Step 2: compute number of bits for initial n + start by assuming sqrt k has coeffs of about b1 bits where + b1 = about half the bits of max_i{abs(coeff(b, i))} + we need nbits >= b1 + 1 (extra 1 for sign) and + bits(m) >= bits(k(n)) where m = f(n) for f = nf->pol + Note f has higher degree than k so this should be true + if bits(n) = b1 + 1. + */ +#if DEBUG + flint_printf("Step 2\n"); +#endif + + bbits = FLINT_ABS(_fmpz_vec_max_bits(NF_ELEM_NUMREF(b), lenb)); + nbits = (bbits + 1)/2 + 2; + + /* + Step 3: find a nbits bit prime such that z = f(n) is a product + at most five distinct primes p_i, none of which divide the + discriminant of f or the norm of b. This guarantees that + P_i = (p_i, x - n) are ideals of degree one in the number + 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(temp); + fmpz_init(n); + + do /* continue increasing nbits until square root found */ + { + fmpz_init(disc); + flint_randinit(state); + + _fmpz_poly_discriminant(disc, fmpq_poly_numref(nf->pol), lenf); + + factored = 0; + + while (!factored || fac->num > 5) /* no bound known for finding such a factorisation */ + { + fmpz_factor_clear(fac); + fmpz_factor_init(fac); + + fmpz_randprime(n, state, nbits, 0); + + _fmpz_poly_evaluate_fmpz(z, fmpq_poly_numref(nf->pol), lenf, n); + + factored = fmpz_factor_trial(fac, z, 3512); + + if (!factored) + factored = fmpz_is_probabprime(fac->p + fac->num - 1); + } + + flint_randclear(state); + fmpz_clear(disc); + + /* + Step 4: compute the square roots r_i of z = b(n) mod p_i for each + of the primes p_i dividing f(n). This allows us to compute + all of the square roots of b(n) mod m where m = prod_i p_i. + If m was sufficiently large, this allows us to retrieve the + square root of b from its n-adic expansion. (Note that z + here is not the same z as above!) + */ +#if DEBUG + flint_printf("Step 4\n"); +#endif + + r = _fmpz_vec_init(fac->num); + + _fmpz_poly_evaluate_fmpz(z, NF_ELEM_NUMREF(b), lenb, n); + + for (i = 0; i < fac->num; i++) + { + fmpz_mod(temp, z, fac->p + i); + + if (!fmpz_sqrtmod(r + i, temp, fac->p + i)) /* check it is square mod p_i */ + { + res = 0; + nf_elem_zero(a, nf); + goto cleanup; + } + } + + /* + Step 5: CRT recombination of the square roots r_i + We compute z congruent to r_i mod p_i, which is hopefully + k(n) where k is the square root of g. Of course we must + go through all possibilities of r_i and -r_i. (Note again + that z here has a different meaning to above.) + */ +#if DEBUG + flint_printf("Step 5\n"); +#endif + + nf_elem_init(sqr, nf); + + mr = _fmpz_vec_init(fac->num); /* compute -r_i mod p_i */ + fmpz_init(m); + + for (i = 0; i < fac->num; i++) + { + if (!fmpz_is_zero(r + i)) + fmpz_sub(mr + i, fac->p + i, r + i); + } + + for (j = 0; j < (1<num) && res != 1; j++) + { + k = j; + + fmpz_set_ui(z, 0); + fmpz_set_ui(m, 1); + nf_elem_zero(a, nf); + + for (i = 0; i < fac->num; i++) + { + if (k & 1) + fmpz_CRT(z, z, m, r + i, fac->p + i, 0); + else + fmpz_CRT(z, z, m, mr + i, fac->p + i, 0); + + fmpz_mul(m, m, fac->p + i); + + k >>= 1; + } + + /* Step 6: retrieve sqrt from n-adic expansion, check sqrt */ +#if DEBUG + flint_printf("Step 6\n"); +#endif + NF_ELEM(a)->length = _fmpz_poly_get_n_adic(NF_ELEM_NUMREF(a), + lenf - 1, z, n); + + nf_elem_mul(sqr, a, a, nf); + + res = nf_elem_equal(sqr, b, nf); + } + + fmpz_clear(m); + _fmpz_vec_clear(mr, fac->num); + + nf_elem_clear(sqr, nf); + + nbits = 1.1*nbits + 1; /* increase nbits and try again if sqrt not found */ + } while (res != 1); + +cleanup: + fmpz_clear(n); + fmpz_clear(temp); + fmpz_clear(z); + + _fmpz_vec_clear(r, fac->num); + + fmpz_factor_clear(fac); + + fmpq_clear(bnorm); + + return res; + } +} + diff --git a/nf_elem/test/t-sqrt.c b/nf_elem/test/t-sqrt.c new file mode 100644 index 000000000..fba2bc6ee --- /dev/null +++ b/nf_elem/test/t-sqrt.c @@ -0,0 +1,99 @@ +/*============================================================================= + + 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 . + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2020 William Hart + +******************************************************************************/ + +#include +#include "nf.h" +#include "nf_elem.h" + +int +main(void) +{ + int i, result; + flint_rand_t state; + + flint_printf("sqrt...."); + fflush(stdout); + + flint_randinit(state); + + /* test sqrt(a^2) */ + for (i = 0; i < 1000 * 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); + + 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 (nf->pol->length > 3 && (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(a, state, 20, 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); + + 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); + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return 0; +} From 54c8fcad328137a6620730a9d708694aa5386f5b Mon Sep 17 00:00:00 2001 From: William Hart Date: Tue, 25 Aug 2020 23:17:56 +0200 Subject: [PATCH 02/24] Speed up. --- nf_elem/sqrt.c | 17 +++++++++++++---- nf_elem/test/t-sqrt.c | 9 ++++----- 2 files changed, 17 insertions(+), 9 deletions(-) diff --git a/nf_elem/sqrt.c b/nf_elem/sqrt.c index a443d173e..c8b7fb166 100755 --- a/nf_elem/sqrt.c +++ b/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); diff --git a/nf_elem/test/t-sqrt.c b/nf_elem/test/t-sqrt.c index fba2bc6ee..8993f72d4 100644 --- a/nf_elem/test/t-sqrt.c +++ b/nf_elem/test/t-sqrt.c @@ -30,7 +30,7 @@ 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; @@ -38,7 +38,7 @@ main(void) 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); From d43075cd8329ec5883be2c99531cf5c197a99078 Mon Sep 17 00:00:00 2001 From: William Hart Date: Wed, 26 Aug 2020 12:21:35 +0200 Subject: [PATCH 03/24] Add to TODO list. --- nf_elem/sqrt.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nf_elem/sqrt.c b/nf_elem/sqrt.c index c8b7fb166..9c78d4d8f 100755 --- a/nf_elem/sqrt.c +++ b/nf_elem/sqrt.c @@ -179,7 +179,7 @@ int nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) factored = 0; iters = 0; - while (!factored || fac->num > 5) /* no bound known for finding such a factorisation */ + while (!factored || fac->num > 10) /* no bound known for finding such a factorisation */ { fmpz_factor_clear(fac); fmpz_factor_init(fac); From 2198263d642e9f9f8e4d37ffd6711fb5a31a144a Mon Sep 17 00:00:00 2001 From: William Hart Date: Thu, 27 Aug 2020 15:27:24 +0200 Subject: [PATCH 04/24] Deal with denominators. --- nf_elem/sqrt.c | 111 +++++++++++++++++++++++++++++++++++++----- nf_elem/test/t-sqrt.c | 4 +- 2 files changed, 101 insertions(+), 14 deletions(-) diff --git a/nf_elem/sqrt.c b/nf_elem/sqrt.c index 9c78d4d8f..214b384ab 100755 --- a/nf_elem/sqrt.c +++ b/nf_elem/sqrt.c @@ -27,8 +27,57 @@ * add LM bound termination for nonsquare case * deal with algebraic integers with denominators * 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 + * 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 */ +int _fmpq_poly_set_fmpz_poly_mod_fmpz(fmpq_poly_t X, + const fmpz * Xmod, slong Xlen, const fmpz_t mod) +{ + fmpz_t t, u, d; + fmpq_t Q; + slong i; + + int success = 1; + + fmpq_init(Q); + fmpz_init(d); + fmpz_init(t); + fmpz_init(u); + + fmpz_one(d); + + fmpq_poly_zero(X); + + for (i = 0; i < Xlen; i++) + { + fmpz_mul(t, d, Xmod + i); + fmpz_fdiv_qr(u, t, t, mod); + + success = _fmpq_reconstruct_fmpz(fmpq_numref(Q), fmpq_denref(Q), t, mod); + + if (!success) + goto cleanup; + + fmpz_mul(fmpq_denref(Q), fmpq_denref(Q), d); + fmpz_set(d, fmpq_denref(Q)); + + fmpq_poly_set_coeff_fmpq(X, i, Q); + } + +cleanup: + fmpq_clear(Q); + fmpz_clear(d); + fmpz_clear(t); + fmpz_clear(u); + + return success; +} + slong _fmpz_poly_get_n_adic(fmpz * sqrt, slong len, fmpz_t z, fmpz_t n) { fmpz_t n2, z2; @@ -98,14 +147,14 @@ int nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) fmpz_t disc, z, temp, n, m; nf_elem_t sqr; slong i, j, k; - fmpz * r, * mr; + fmpz * r, * mr, * bz; int res = 0, factored, iters; - if (!fmpz_is_one(NF_ELEM_DENREF(b))) + /* if (!fmpz_is_one(NF_ELEM_DENREF(b))) { flint_printf("Sqrt1 outside Z[alpha] not implemented yet\n"); flint_abort(); - } + }*/ if (lenb == 0) { @@ -113,7 +162,7 @@ int nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) return 1; } - /* Step 1: compute norm and check it is square */ + /* Step 1: compute norm and check it is square and rationalise denominator */ #if DEBUG flint_printf("Step 1\n"); #endif @@ -122,20 +171,41 @@ int nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) nf_elem_norm(bnorm, b, nf); - if (!fmpz_is_one(fmpq_denref(bnorm))) + /* if (!fmpz_is_one(fmpq_denref(bnorm))) { flint_printf("Sqrt2 outside Z[alpha] not yet implemented yet\n"); fmpq_clear(bnorm); flint_abort(); - } + } */ +/* if (!fmpz_is_square(fmpq_numref(bnorm))) { nf_elem_zero(a, nf); fmpq_clear(bnorm); +#if DEBUG + flint_printf("Norm does not have square numerator\n"); +#endif return 0; } + if (!fmpz_is_square(fmpq_denref(bnorm))) + { + nf_elem_zero(a, nf); + fmpq_clear(bnorm); +#if DEBUG + flint_printf("Norm does not have square denominator\n"); +#endif + return 0; + } +*/ + + fmpz_init(temp); + + /* get rid of denominator */ + bz = _fmpz_vec_init(NF_ELEM(b)->length); + _fmpz_vec_scalar_mul_fmpz(bz, NF_ELEM_NUMREF(b), NF_ELEM(b)->length, NF_ELEM_DENREF(b)); + /* Step 2: compute number of bits for initial n start by assuming sqrt k has coeffs of about b1 bits where @@ -149,7 +219,7 @@ int nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) flint_printf("Step 2\n"); #endif - bbits = FLINT_ABS(_fmpz_vec_max_bits(NF_ELEM_NUMREF(b), lenb)); + bbits = FLINT_ABS(_fmpz_vec_max_bits(bz, lenb)); nbits = (bbits + 1)/(2*lenf) + 2; /* @@ -166,7 +236,6 @@ int nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) fmpz_factor_init(fac); fmpz_init(z); - fmpz_init(temp); fmpz_init(n); do /* continue increasing nbits until square root found */ @@ -179,13 +248,15 @@ int nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) factored = 0; iters = 0; - while (!factored || fac->num > 10) /* no bound known for finding such a factorisation */ + while (!factored || fac->num > 6) /* no bound known for finding such a factorisation */ { fmpz_factor_clear(fac); fmpz_factor_init(fac); fmpz_randprime(n, state, nbits, 0); - + if (fmpz_sgn(n) < 0) + fmpz_neg(n, n); + _fmpz_poly_evaluate_fmpz(z, fmpq_poly_numref(nf->pol), lenf, n); factored = fmpz_factor_trial(fac, z, 3512); @@ -219,7 +290,7 @@ int nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) r = _fmpz_vec_init(fac->num); - _fmpz_poly_evaluate_fmpz(z, NF_ELEM_NUMREF(b), lenb, n); + _fmpz_poly_evaluate_fmpz(z, bz, lenb, n); for (i = 0; i < fac->num; i++) { @@ -282,9 +353,26 @@ int nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) NF_ELEM(a)->length = _fmpz_poly_get_n_adic(NF_ELEM_NUMREF(a), lenf - 1, z, n); + fmpz_set(NF_ELEM_DENREF(a), NF_ELEM_DENREF(b)); + nf_elem_mul(sqr, a, a, nf); res = nf_elem_equal(sqr, b, nf); + + if (!res) /* try rational coeffs */ + { + res = _fmpq_poly_set_fmpz_poly_mod_fmpz(NF_ELEM(a), + NF_ELEM_NUMREF(a), NF_ELEM(a)->length, n); + + if (res) + { + fmpz_mul(NF_ELEM_DENREF(a), NF_ELEM_DENREF(a), NF_ELEM_DENREF(b)); + + nf_elem_mul(sqr, a, a, nf); + + res = nf_elem_equal(sqr, b, nf); + } + } } fmpz_clear(m); @@ -301,6 +389,7 @@ int nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) fmpz_clear(z); _fmpz_vec_clear(r, fac->num); + _fmpz_vec_clear(bz, NF_ELEM(b)->length); fmpz_factor_clear(fac); diff --git a/nf_elem/test/t-sqrt.c b/nf_elem/test/t-sqrt.c index 8993f72d4..774fd9fd4 100644 --- a/nf_elem/test/t-sqrt.c +++ b/nf_elem/test/t-sqrt.c @@ -52,7 +52,7 @@ main(void) fmpz_poly_factor_clear(fac); - if (nf->pol->length > 3 && (nf->flag & NF_MONIC) && num_facs == 1) + if (nf->pol->length > 3 && nf->flag & NF_MONIC && num_facs == 1) { i++; @@ -63,8 +63,6 @@ main(void) nf_elem_randtest(a, state, 30, nf); - fmpz_set_ui(NF_ELEM_DENREF(a), 1); - nf_elem_mul(b, a, a, nf); is_square = nf_elem_sqrt(c, b, nf); From 04263bda596fb9ff1f6eff38021c0408b370003d Mon Sep 17 00:00:00 2001 From: William Hart Date: Thu, 27 Aug 2020 15:33:15 +0200 Subject: [PATCH 05/24] Add back norm test. --- nf_elem/sqrt.c | 2 -- 1 file changed, 2 deletions(-) diff --git a/nf_elem/sqrt.c b/nf_elem/sqrt.c index 214b384ab..6b4292d67 100755 --- a/nf_elem/sqrt.c +++ b/nf_elem/sqrt.c @@ -178,7 +178,6 @@ int nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) flint_abort(); } */ -/* if (!fmpz_is_square(fmpq_numref(bnorm))) { nf_elem_zero(a, nf); @@ -198,7 +197,6 @@ int nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) #endif return 0; } -*/ fmpz_init(temp); From 693e881488400087dfe0626181bcbc86b6fe9a63 Mon Sep 17 00:00:00 2001 From: William Hart Date: Thu, 27 Aug 2020 15:38:36 +0200 Subject: [PATCH 06/24] Cleanup. --- nf_elem/sqrt.c | 16 +++++----------- 1 file changed, 5 insertions(+), 11 deletions(-) diff --git a/nf_elem/sqrt.c b/nf_elem/sqrt.c index 6b4292d67..88dc4c5ba 100755 --- a/nf_elem/sqrt.c +++ b/nf_elem/sqrt.c @@ -25,7 +25,6 @@ * try to reuse information from previous failed attempt * improve bounds * add LM bound termination for nonsquare case - * deal with algebraic integers with denominators * add linear and quadratic cases * Tune the number of primes used in trial factoring * Use ECM and larger recombination for very large square roots @@ -33,6 +32,8 @@ * 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 + * Figure out why norm is square test fails if defining polynomial is not monic + * Remove small squares from denominator before rationalising */ int _fmpq_poly_set_fmpz_poly_mod_fmpz(fmpq_poly_t X, @@ -150,11 +151,11 @@ 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(NF_ELEM_DENREF(b))) + if (!fmpz_is_one(fmpq_poly_denref(nf->pol))) { - flint_printf("Sqrt1 outside Z[alpha] not implemented yet\n"); + flint_printf("Non-monic defining polynomial not supported in sqrt yet.\n"); flint_abort(); - }*/ + } if (lenb == 0) { @@ -170,13 +171,6 @@ int nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) fmpq_init(bnorm); nf_elem_norm(bnorm, b, nf); - - /* if (!fmpz_is_one(fmpq_denref(bnorm))) - { - flint_printf("Sqrt2 outside Z[alpha] not yet implemented yet\n"); - fmpq_clear(bnorm); - flint_abort(); - } */ if (!fmpz_is_square(fmpq_numref(bnorm))) { From 752cbebd613b2c1a8391b3f5dc6ef79c4353c3c2 Mon Sep 17 00:00:00 2001 From: William Hart Date: Thu, 27 Aug 2020 15:39:47 +0200 Subject: [PATCH 07/24] Add comment. --- nf_elem/sqrt.c | 1 + 1 file changed, 1 insertion(+) diff --git a/nf_elem/sqrt.c b/nf_elem/sqrt.c index 88dc4c5ba..51b5d3ce1 100755 --- a/nf_elem/sqrt.c +++ b/nf_elem/sqrt.c @@ -34,6 +34,7 @@ * Deal with primes dividing denominator of norm * Figure out why norm is square test fails if defining polynomial is not monic * Remove small squares from denominator before rationalising + * Move _fmpq_poly_set_fmpz_poly_mod_fmpz into Flint */ int _fmpq_poly_set_fmpz_poly_mod_fmpz(fmpq_poly_t X, From 41424a4b7d14da43d93f9c05ac777d3e8aa683d4 Mon Sep 17 00:00:00 2001 From: William Hart Date: Thu, 27 Aug 2020 15:53:43 +0200 Subject: [PATCH 08/24] Comment. --- nf_elem/sqrt.c | 1 + 1 file changed, 1 insertion(+) diff --git a/nf_elem/sqrt.c b/nf_elem/sqrt.c index 51b5d3ce1..eae2d8336 100755 --- a/nf_elem/sqrt.c +++ b/nf_elem/sqrt.c @@ -35,6 +35,7 @@ * Figure out why norm is square test fails if defining polynomial is not monic * Remove small squares from denominator before rationalising * Move _fmpq_poly_set_fmpz_poly_mod_fmpz into Flint + * Cache factorisation of f(n) on number field for future square roots */ int _fmpq_poly_set_fmpz_poly_mod_fmpz(fmpq_poly_t X, From e26a357fbd0c1b8544ed0575e3522344d7f77849 Mon Sep 17 00:00:00 2001 From: William Hart Date: Thu, 27 Aug 2020 17:21:17 +0200 Subject: [PATCH 09/24] More comments. --- nf_elem/sqrt.c | 1 + 1 file changed, 1 insertion(+) diff --git a/nf_elem/sqrt.c b/nf_elem/sqrt.c index eae2d8336..5e038bcc0 100755 --- a/nf_elem/sqrt.c +++ b/nf_elem/sqrt.c @@ -36,6 +36,7 @@ * Remove small squares from denominator before rationalising * Move _fmpq_poly_set_fmpz_poly_mod_fmpz into Flint * Cache factorisation of f(n) on number field for future square roots + * Deal with primality vs probable prime testing */ int _fmpq_poly_set_fmpz_poly_mod_fmpz(fmpq_poly_t X, From 07dcca885d5604deff469b62fc9be838d7a6bf27 Mon Sep 17 00:00:00 2001 From: William Hart Date: Wed, 16 Sep 2020 16:04:15 +0200 Subject: [PATCH 10/24] 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); From 3169ec57d4d64e3ff49a751a2f6a8b84e0e31c0e Mon Sep 17 00:00:00 2001 From: William Hart Date: Wed, 16 Sep 2020 16:21:27 +0200 Subject: [PATCH 11/24] Use fewer primes in trial factoring. --- nf_elem/sqrt.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nf_elem/sqrt.c b/nf_elem/sqrt.c index de447137e..f980b9e7d 100755 --- a/nf_elem/sqrt.c +++ b/nf_elem/sqrt.c @@ -259,7 +259,7 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) _fmpz_poly_evaluate_fmpz(z, fmpq_poly_numref(nf->pol), lenf, n); - factored = fmpz_factor_trial(fac, z, 3512); + factored = fmpz_factor_trial(fac, z, FLINT_MIN(log(nbits)*nbits, 3512)); if (!factored) { From cc337a3e1d1e0495fb5888d5a2c85bc0f4d01add Mon Sep 17 00:00:00 2001 From: William Hart Date: Wed, 16 Sep 2020 16:45:44 +0200 Subject: [PATCH 12/24] Add linear case. --- nf_elem/sqrt.c | 42 ++++++++++++++++++++++++++++++++---------- nf_elem/test/t-sqrt.c | 2 +- 2 files changed, 33 insertions(+), 11 deletions(-) diff --git a/nf_elem/sqrt.c b/nf_elem/sqrt.c index f980b9e7d..c9e7e6326 100755 --- a/nf_elem/sqrt.c +++ b/nf_elem/sqrt.c @@ -24,9 +24,7 @@ * 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 + * add quadratic 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) * Deal with number fields of degree 1 and 2 @@ -36,6 +34,8 @@ * Move _fmpq_poly_set_fmpz_poly_mod_fmpz into Flint * Cache factorisation of f(n) on number field for future square roots * Deal with primality vs probable prime testing + * add is_square function + * test non-squares */ int _fmpq_poly_set_fmpz_poly_mod_fmpz(fmpq_poly_t X, @@ -124,13 +124,31 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) { if (nf->flag & NF_LINEAR) { - /* const fmpz * const bnum = LNF_ELEM_NUMREF(b); + const fmpz * const bnum = LNF_ELEM_NUMREF(b); const fmpz * const bden = LNF_ELEM_DENREF(b); fmpz * const anum = LNF_ELEM_NUMREF(a); - fmpz * const aden = LNF_ELEM_DENREF(a); */ + fmpz * const aden = LNF_ELEM_DENREF(a); + fmpz_t r; + int res = 1; + + fmpz_init(r); - flint_printf("Sqrt for linear number fields not implemented yet\n"); - flint_abort(); + fmpz_sqrtrem(anum, r, bnum); + + if (!fmpz_is_zero(r)) + res = 0; + + fmpz_sqrtrem(aden, r, bden); + + if (!fmpz_is_zero(r)) + res = 0; + + if (!res) + nf_elem_zero(a, nf); + + fmpz_clear(r); + + return res; } else if (nf->flag & NF_QUADRATIC) { /* const fmpz * const bnum = QNF_ELEM_NUMREF(b); @@ -296,9 +314,13 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) if (!fmpz_sqrtmod(r + i, temp, fac->p + i)) /* check it is square mod p_i */ { - res = 0; - nf_elem_zero(a, nf); - goto cleanup; + /* we must prove primality of prime used to reject as nonsquare */ + if (fmpz_is_prime(fac->p + i)) + { + res = 0; + nf_elem_zero(a, nf); + goto cleanup; + } } } diff --git a/nf_elem/test/t-sqrt.c b/nf_elem/test/t-sqrt.c index 15e332b78..5e14297a5 100644 --- a/nf_elem/test/t-sqrt.c +++ b/nf_elem/test/t-sqrt.c @@ -57,7 +57,7 @@ main(void) fmpz_poly_factor_clear(fac); - if (nf->pol->length > 3 && nf->flag & NF_MONIC && num_facs == 1) + if (nf->pol->length != 3 && nf->flag & NF_MONIC && num_facs == 1) { i++; From 71f71c17309732fdee6f9ab9f824cf9f27d959a6 Mon Sep 17 00:00:00 2001 From: William Hart Date: Thu, 17 Sep 2020 03:02:36 +0200 Subject: [PATCH 13/24] Add code for square root in quadratic fields. --- nf_elem/sqrt.c | 190 ++++++++++++++++++++++++++++++++++++++++-- nf_elem/test/t-sqrt.c | 4 +- 2 files changed, 186 insertions(+), 8 deletions(-) diff --git a/nf_elem/sqrt.c b/nf_elem/sqrt.c index c9e7e6326..2f45f44b4 100755 --- a/nf_elem/sqrt.c +++ b/nf_elem/sqrt.c @@ -36,6 +36,9 @@ * Deal with primality vs probable prime testing * add is_square function * test non-squares + * deal with lifting of square roots mod prime factors of f(n) to same + exponent as in the factorisation of f(n) (e.g. the square root with inputs + f = x^8 - 8, g^2 = 4050*x^6 mod f fails currently) */ int _fmpq_poly_set_fmpz_poly_mod_fmpz(fmpq_poly_t X, @@ -151,13 +154,185 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) return res; } else if (nf->flag & NF_QUADRATIC) { - /* const fmpz * const bnum = QNF_ELEM_NUMREF(b); + const fmpz * const bnum = QNF_ELEM_NUMREF(b); const fmpz * const bden = QNF_ELEM_DENREF(b); fmpz * const anum = QNF_ELEM_NUMREF(a); - fmpz * const aden = QNF_ELEM_DENREF(a); */ + fmpz * const aden = QNF_ELEM_DENREF(a); + fmpz_t d, t, t2, t3; + fmpq_t r, s, r1, s1, tq, tq2; + nf_elem_t sqr; + int sq, ret = 1; + + /* + We use that (m + n*sqrt(d))^2 = m^2 + n^2*d + 2*m*n*sqrt(d). + Thus, finding (m^2)*(n^2*d) and (m^2) + (n^2*d) allows us to + retrieve m^2 and n^2*d as the roots of (Y - m^2)(Y - n^2*d), + from which we can retrieve m and n. + + Here d is the discriminant of the defining polynomial of the + number field. + */ + fmpq_init(r); + fmpq_init(s); + fmpq_init(r1); + fmpq_init(s1); + fmpq_init(tq); + fmpq_init(tq2); + fmpz_init(t); + fmpz_init(t2); + fmpz_init(t3); + fmpz_init(d); + nf_elem_init(sqr, nf); + + /* compute discriminant d of the defining polynomial */ + fmpz_mul(d, fmpq_poly_numref(nf->pol) + 1, fmpq_poly_numref(nf->pol) + 1); + fmpz_mul(t, fmpq_poly_numref(nf->pol) + 0, fmpq_poly_numref(nf->pol) + 2); + fmpz_mul_2exp(t, t, 2); + fmpz_sub(d, d, t); + + /* write x mod nf->pol in the form r1 + s1*sqrt(d) */ + fmpz_neg(fmpq_numref(r1), fmpq_poly_numref(nf->pol) + 1); + fmpz_mul_2exp(fmpq_denref(r1), fmpq_poly_numref(nf->pol) + 2, 1); + fmpz_set(fmpq_denref(s1), fmpq_denref(r1)); + fmpz_set_ui(fmpq_numref(s1), 1); + fmpq_canonicalise(r1); + fmpq_canonicalise(s1); + + /* write b in the form r + s*sqrt(d) */ + fmpq_set_fmpz_frac(tq, bnum + 1, bden); + fmpq_mul(r, r1, tq); + fmpq_mul(s, s1, tq); + fmpq_set_fmpz_frac(tq, bnum + 0, bden); + fmpq_add(r, r, tq); + + /* compute m^2*n^2*d and m^2 + n^2*d as above */ + fmpq_div_2exp(s, s, 1); + fmpq_mul(s, s, s); + fmpq_mul_fmpz(s, s, d); + + /* compute m^2 and n^2*d as above */ + fmpq_mul(tq, r, r); + fmpq_mul_2exp(tq2, s, 2); + fmpq_sub(tq, tq, tq2); + fmpz_sqrtrem(fmpq_numref(s), t, fmpq_numref(tq)); + if (!fmpz_is_zero(t)) + { + ret = 0; + nf_elem_zero(a, nf); + + goto quadratic_cleanup; + } + fmpz_sqrtrem(fmpq_denref(s), t, fmpq_denref(tq)); + if (!fmpz_is_zero(t)) + { + ret = 0; + nf_elem_zero(a, nf); + + goto quadratic_cleanup; + } + fmpq_div_2exp(r, r, 1); + fmpq_div_2exp(s, s, 1); - flint_printf("Sqrt for quadratic number fields not implemented yet\n"); - flint_abort(); + /* + we now have m^2 = r + s and n^2*d = r - s, or vice versa, so + compute m^2 and n^2 + */ + fmpq_add(tq, r, s); + fmpq_sub(s, r, s); + fmpq_swap(r, tq); + fmpq_div_fmpz(s, s, d); + + /* compute m and +/- n */ + fmpz_sqrtrem(t2, t, fmpq_numref(r)); + sq = fmpz_is_zero(t); + if (sq) + { + fmpz_sqrtrem(t3, t, fmpq_denref(r)); + sq = fmpz_is_zero(t); + } + if (!sq) + { + fmpq_mul_fmpz(s, s, d); + fmpq_div_fmpz(r, r, d); + fmpq_swap(r, s); + fmpz_sqrtrem(t2, t, fmpq_numref(r)); + sq = fmpz_is_zero(t); + if (sq) + { + fmpz_sqrtrem(t3, t, fmpq_denref(r)); + sq &= fmpz_is_zero(t); + } + if (!sq) + { + ret = 0; + nf_elem_zero(a, nf); + + goto quadratic_cleanup; + } + } + fmpz_swap(fmpq_numref(r), t2); + fmpz_swap(fmpq_denref(r), t3); + fmpz_sqrtrem(fmpq_numref(s), t, fmpq_numref(s)); + if (!fmpz_is_zero(t)) + { + ret = 0; + nf_elem_zero(a, nf); + + goto quadratic_cleanup; + } + fmpz_sqrtrem(fmpq_denref(s), t, fmpq_denref(s)); + if (!fmpz_is_zero(t)) + { + ret = 0; + nf_elem_zero(a, nf); + + goto quadratic_cleanup; + } + + /* write m + n*sqrt(d) as alpha*x + beta */ + fmpq_div(tq, s, s1); + fmpq_mul(tq2, tq, r1); + fmpq_sub(tq2, r, tq2); + + fmpz_gcd(t, fmpq_denref(tq), fmpq_denref(tq2)); + fmpz_mul(aden, fmpq_denref(tq), fmpq_denref(tq2)); + fmpz_tdiv_q(aden, aden, t); + fmpz_tdiv_q(anum + 1, fmpq_denref(tq2), t); + fmpz_tdiv_q(anum + 0, fmpq_denref(tq), t); + fmpz_mul(anum + 1, anum + 1, fmpq_numref(tq)); + fmpz_mul(anum + 0, anum + 0, fmpq_numref(tq2)); + + nf_elem_mul(sqr, a, a, nf); + + if (!nf_elem_equal(sqr, b, nf)) + { + fmpq_mul_2exp(r, r, 1); + fmpq_sub(tq2, tq2, r); + + fmpz_gcd(t, fmpq_denref(tq), fmpq_denref(tq2)); + fmpz_mul(aden, fmpq_denref(tq), fmpq_denref(tq2)); + fmpz_tdiv_q(aden, aden, t); + fmpz_tdiv_q(anum + 1, fmpq_denref(tq2), t); + fmpz_tdiv_q(anum + 0, fmpq_denref(tq), t); + fmpz_mul(anum + 1, anum + 1, fmpq_numref(tq)); + fmpz_mul(anum + 0, anum + 0, fmpq_numref(tq2)); + } + +quadratic_cleanup: + + fmpq_clear(r); + fmpq_clear(s); + fmpq_clear(r1); + fmpq_clear(s1); + fmpq_clear(tq); + fmpq_clear(tq2); + fmpz_clear(t); + fmpz_clear(t2); + fmpz_clear(t3); + fmpz_clear(d); + nf_elem_clear(sqr, nf); + + return ret; } else /* generic nf_elem */ { const slong lenb = NF_ELEM(b)->length, lenf = fmpq_poly_length(nf->pol); @@ -261,6 +436,8 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) while (!factored || fac->num > 14) /* no bound known for finding such a factorisation */ { + slong primes; + fmpz_factor_clear(fac); fmpz_factor_init(fac); @@ -277,14 +454,15 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) _fmpz_poly_evaluate_fmpz(z, fmpq_poly_numref(nf->pol), lenf, n); - factored = fmpz_factor_trial(fac, z, FLINT_MIN(log(nbits)*nbits, 3512)); + primes = FLINT_MIN((slong) log(nbits)*nbits + 2, 3512); + factored = fmpz_factor_trial(fac, z, primes); if (!factored) { fmpz_set(fac1, fac->p + fac->num - 1); fac->num--; - factored = fmpz_factor_smooth(fac, fac1, FLINT_MIN(20, nbits/5 + 1), 0); + factored = fmpz_factor_smooth(fac, fac1, FLINT_MIN(20, fmpz_bits(z)/100 + 1), 0); } iters++; diff --git a/nf_elem/test/t-sqrt.c b/nf_elem/test/t-sqrt.c index 5e14297a5..2a7f61745 100644 --- a/nf_elem/test/t-sqrt.c +++ b/nf_elem/test/t-sqrt.c @@ -57,7 +57,7 @@ main(void) fmpz_poly_factor_clear(fac); - if (nf->pol->length != 3 && nf->flag & NF_MONIC && num_facs == 1) + if (nf->flag & NF_MONIC && num_facs == 1) { i++; @@ -69,7 +69,7 @@ main(void) nf_elem_randtest(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); From 9793d163b6f1e067e016a53e40225d7cb7a3a7f0 Mon Sep 17 00:00:00 2001 From: William Hart Date: Thu, 17 Sep 2020 03:07:34 +0200 Subject: [PATCH 14/24] Revise todos. --- nf_elem/sqrt.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nf_elem/sqrt.c b/nf_elem/sqrt.c index 2f45f44b4..07f5e9b61 100755 --- a/nf_elem/sqrt.c +++ b/nf_elem/sqrt.c @@ -24,7 +24,6 @@ * try to reuse information from previous failed attempt * improve bounds * add LM bound termination for nonsquare case - * add quadratic 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) * Deal with number fields of degree 1 and 2 @@ -39,6 +38,7 @@ * deal with lifting of square roots mod prime factors of f(n) to same exponent as in the factorisation of f(n) (e.g. the square root with inputs f = x^8 - 8, g^2 = 4050*x^6 mod f fails currently) + * fix bug in fmpz_factor_trial #843 */ int _fmpq_poly_set_fmpz_poly_mod_fmpz(fmpq_poly_t X, From ee2820ac12cd1889385a916a6565a6929ca98523 Mon Sep 17 00:00:00 2001 From: William Hart Date: Tue, 29 Sep 2020 11:48:43 +0200 Subject: [PATCH 15/24] Fix four bugs in algorithm. --- nf_elem/sqrt.c | 287 +++++++++++++++++++++++++++++------------- nf_elem/test/t-sqrt.c | 214 ++++++++++++++++++++++++++++++- 2 files changed, 413 insertions(+), 88 deletions(-) diff --git a/nf_elem/sqrt.c b/nf_elem/sqrt.c index 07f5e9b61..43086891d 100755 --- a/nf_elem/sqrt.c +++ b/nf_elem/sqrt.c @@ -23,66 +23,20 @@ * try to reuse information from previous failed attempt * improve bounds - * add LM bound termination for nonsquare case + * 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) - * Deal with number fields of degree 1 and 2 - * Deal with primes dividing denominator of norm - * Figure out why norm is square test fails if defining polynomial is not monic + * Fix bug in antic norm function * Remove small squares from denominator before rationalising - * Move _fmpq_poly_set_fmpz_poly_mod_fmpz into Flint * Cache factorisation of f(n) on number field for future square roots - * Deal with primality vs probable prime testing * add is_square function * test non-squares - * deal with lifting of square roots mod prime factors of f(n) to same - exponent as in the factorisation of f(n) (e.g. the square root with inputs - f = x^8 - 8, g^2 = 4050*x^6 mod f fails currently) * fix bug in fmpz_factor_trial #843 + * valgrind */ -int _fmpq_poly_set_fmpz_poly_mod_fmpz(fmpq_poly_t X, - const fmpz * Xmod, slong Xlen, const fmpz_t mod) -{ - fmpz_t t, u, d; - fmpq_t Q; - slong i; - - int success = 1; - - fmpq_init(Q); - fmpz_init(d); - fmpz_init(t); - fmpz_init(u); - - fmpz_one(d); - - fmpq_poly_zero(X); - - for (i = 0; i < Xlen; i++) - { - fmpz_mul(t, d, Xmod + i); - fmpz_fdiv_qr(u, t, t, mod); - - success = _fmpq_reconstruct_fmpz(fmpq_numref(Q), fmpq_denref(Q), t, mod); - - if (!success) - goto cleanup; - - fmpz_mul(fmpq_denref(Q), fmpq_denref(Q), d); - fmpz_set(d, fmpq_denref(Q)); - - fmpq_poly_set_coeff_fmpq(X, i, Q); - } - -cleanup: - fmpq_clear(Q); - fmpz_clear(d); - fmpz_clear(t); - fmpz_clear(u); - - return success; -} +#define ROT(u,v,t) \ + do { fmpz _t = *u; *u = *v; *v = *t; *t = _t; } while (0); slong _fmpz_poly_get_n_adic(fmpz * sqrt, slong len, fmpz_t z, fmpz_t n) { @@ -340,10 +294,10 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) fmpq_t bnorm; fmpz_factor_t fac; flint_rand_t state; - fmpz_t disc, z, temp, n, m; + fmpz_t disc, z, temp, n, m, az, d; nf_elem_t sqr; slong i, j, k; - fmpz * r, * mr, * bz; + fmpz * r, * mr, * bz, * bz1, * modulus; int res = 0, factored, iters; if (lenb == 0) @@ -384,9 +338,10 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) fmpz_init(temp); /* get rid of denominator */ + bz1 = _fmpz_vec_init(NF_ELEM(b)->length); bz = _fmpz_vec_init(NF_ELEM(b)->length); - _fmpz_vec_scalar_mul_fmpz(bz, NF_ELEM_NUMREF(b), NF_ELEM(b)->length, NF_ELEM_DENREF(b)); - + _fmpz_vec_scalar_mul_fmpz(bz1, NF_ELEM_NUMREF(b), NF_ELEM(b)->length, NF_ELEM_DENREF(b)); + /* Step 2: compute number of bits for initial n start by assuming sqrt k has coeffs of about b1 bits where @@ -405,9 +360,8 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) /* Step 3: find a nbits bit prime such that z = f(n) is a product - at most five distinct primes p_i, none of which divide the - discriminant of f or the norm of b. This guarantees that - P_i = (p_i, x - n) are ideals of degree one in the number + of powers of at most 14 distinct primes p_i. This ensures the + P_i = (p_i, x - n) are ideals of degree one in the number 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. */ @@ -415,6 +369,8 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) fmpz_factor_init(fac); fmpz_init(z); fmpz_init(n); + fmpz_init(d); + fmpz_set_ui(d, 1); fmpz_init(disc); flint_randinit(state); @@ -427,6 +383,8 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) fmpz_init(fac1); + _fmpz_vec_set(bz, bz1, NF_ELEM(b)->length); + factored = 0; iters = 0; @@ -438,6 +396,7 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) { slong primes; + fmpz_set_ui(d, 1); fmpz_factor_clear(fac); fmpz_factor_init(fac); @@ -448,7 +407,7 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) nbits++; } - fmpz_randprime(n, state, nbits, 0); + fmpz_randbits(n, state, nbits); if (fmpz_sgn(n) < 0) fmpz_neg(n, n); @@ -459,24 +418,78 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) if (!factored) { + fmpz_factor_t fac2; + fmpz_set(fac1, fac->p + fac->num - 1); fac->num--; - factored = fmpz_factor_smooth(fac, fac1, FLINT_MIN(20, fmpz_bits(z)/100 + 1), 0); + fmpz_factor_init(fac2); + factored = fmpz_factor_smooth(fac2, fac1, FLINT_MIN(20, fmpz_bits(z)/100 + 1), 0); + _fmpz_factor_concat(fac, fac2, fac->exp[fac->num]); + fmpz_factor_clear(fac2); + } + + if (factored) + { + /* check for primes dividing the discriminant */ + j = 0; + for (i = 0; i < fac->num; i++) + { + fmpz_mod(temp, disc, fac->p + i); + + if (fmpz_is_zero(temp)) + { + if (fmpz_cmp_ui(fac->p + i, 2*lenf) > 0) + { + factored = 0; + + break; + } else + { + if (fac->exp[i] & 1) + fac->exp[i]++; + + fac->exp[i] >>= 1; + + fmpz_pow_ui(temp, fac->p + i, fac->exp[i]); + fmpz_mul(d, d, temp); + fmpz_mul(temp, temp, temp); + + fmpz_zero(fac->p + i); + } + } else + { + if (i != j) + { + fmpz_set(fac->p + j, fac->p + i); + fac->exp[j] = fac->exp[i]; + } + j++; + } + } + + if (i == fac->num) + fac->num = j; } iters++; } + if (!fmpz_is_one(d)) + { + _fmpz_vec_scalar_mul_fmpz(bz, bz, NF_ELEM(b)->length, d); + _fmpz_vec_scalar_mul_fmpz(bz, bz, NF_ELEM(b)->length, d); + } + fmpz_clear(fac1); /* - Step 4: compute the square roots r_i of z = b(n) mod p_i for each - of the primes p_i dividing f(n). This allows us to compute - all of the square roots of b(n) mod m where m = prod_i p_i. - If m was sufficiently large, this allows us to retrieve the - square root of b from its n-adic expansion. (Note that z - here is not the same z as above!) + Step 4: compute the square roots r_i of z = b(n) mod p_i^k_i for + each of the prime powers p_i^k_i dividing f(n). This allows + us to compute all of the square roots of b(n) mod m where + m = prod_i p_i^k_i. If m was sufficiently large, this + allows us to retrieve the square root of b from its n-adic + expansion. (Note that z here is not the same z as above!) */ #if DEBUG flint_printf("Step 4\n"); @@ -486,6 +499,10 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) _fmpz_poly_evaluate_fmpz(z, bz, lenb, n); + fmpz_init(m); + fmpz_init(az); + modulus = _fmpz_vec_init(fac->num); + for (i = 0; i < fac->num; i++) { fmpz_mod(temp, z, fac->p + i); @@ -497,17 +514,47 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) { res = 0; nf_elem_zero(a, nf); + _fmpz_vec_clear(modulus, fac->num); goto cleanup; } } + + fmpz_set(modulus + i, fac->p + i); /* initial modulus m = p_i */ + + /* + Step 4a: lift the square roots mod p_i to roots mod p_i^k_i + We require the lift to be unique so must exclude 2 + and the case where the root is 0 mod p_i from lifting + */ + if (fac->exp[i] != 1 && fmpz_cmp_ui(fac->p + i, 2) != 0 && !fmpz_is_zero(r + i)) + { + for (j = 1; j < fac->exp[i]; j++) + { + fmpz_add(az, r + i, r + i); + if (fmpz_cmp(az, modulus + i) >= 0) + fmpz_sub(az, az, modulus + i); + + fmpz_mul(m, r + i, r + i); + fmpz_sub(m, m, z); + fmpz_mul(modulus + i, modulus + i, fac->p + i); + fmpz_mod(m, m, modulus + i); + fmpz_mul(m, m, az); + fmpz_mod(m, m, modulus + i); + fmpz_sub(r + i, r + i, m); + if (fmpz_sgn(r + i) < 0) + fmpz_add(r + i, r + i, modulus + i); + } + } } + fmpz_clear(az); + /* Step 5: CRT recombination of the square roots r_i - We compute z congruent to r_i mod p_i, which is hopefully - k(n) where k is the square root of g. Of course we must - go through all possibilities of r_i and -r_i. (Note again - that z here has a different meaning to above.) + We compute z congruent to r_i mod p_i^k_i, which is + hopefully k(n) where k is the square root of g. Of course + we must go through all possibilities of r_i and -r_i. (Note + again that z here has a different meaning to above.) */ #if DEBUG flint_printf("Step 5\n"); @@ -515,13 +562,12 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) nf_elem_init(sqr, nf); - mr = _fmpz_vec_init(fac->num); /* compute -r_i mod p_i */ - fmpz_init(m); + mr = _fmpz_vec_init(fac->num); /* compute -r_i mod p_i^k_i */ for (i = 0; i < fac->num; i++) { if (!fmpz_is_zero(r + i)) - fmpz_sub(mr + i, fac->p + i, r + i); + fmpz_sub(mr + i, modulus + i, r + i); } for (j = 0; j < (1<num) && res != 1; j++) @@ -535,16 +581,19 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) for (i = 0; i < fac->num; i++) { if (k & 1) - fmpz_CRT(z, z, m, r + i, fac->p + i, 0); + fmpz_CRT(z, z, m, r + i, modulus + i, 0); else - fmpz_CRT(z, z, m, mr + i, fac->p + i, 0); + fmpz_CRT(z, z, m, mr + i, modulus + i, 0); - fmpz_mul(m, m, fac->p + i); + fmpz_mul(m, m, modulus + i); k >>= 1; } - /* Step 6: retrieve sqrt from n-adic expansion, check sqrt */ + /* + Step 6: retrieve sqrt from n-adic expansion, check sqrt. + First try assuming coeffs are integers. + */ #if DEBUG flint_printf("Step 6\n"); #endif @@ -552,30 +601,94 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) lenf - 1, z, n); fmpz_set(NF_ELEM_DENREF(a), NF_ELEM_DENREF(b)); + + if (!fmpz_is_one(d)) + { + fmpz_mul(NF_ELEM_DENREF(a), NF_ELEM_DENREF(a), d); + fmpq_poly_canonicalise(NF_ELEM(a)); + } nf_elem_mul(sqr, a, a, nf); res = nf_elem_equal(sqr, b, nf); - if (!res) /* try rational coeffs */ + /* + Step 6a. Now try rational coeffs. We loop through all + possible denominators of the *evaluation* k(n), + and try to compute the numerator as an integer + polynomial. + */ + if (!res) { - res = _fmpq_poly_set_fmpz_poly_mod_fmpz(NF_ELEM(a), - NF_ELEM_NUMREF(a), NF_ELEM(a)->length, n); - - if (res) + fmpz_t q, r, s, t, g, num, den; + int sgn; + + fmpz_init(q); + fmpz_init(r); + fmpz_init(s); + fmpz_init(t); + fmpz_init(g); + fmpz_init(num); + fmpz_init(den); + + fmpz_set(r, m); fmpz_zero(s); + fmpz_set(num, z); fmpz_one(den); + + while (!res && !fmpz_is_zero(num)) /* try rational coeffs */ { - fmpz_mul(NF_ELEM_DENREF(a), NF_ELEM_DENREF(a), NF_ELEM_DENREF(b)); + fmpz_fdiv_q(q, r, num); + fmpz_mul(t, q, num); fmpz_sub(t, r, t); ROT(r, num, t); + fmpz_mul(t, q, den); fmpz_sub(t, s, t); ROT(s, den, t); + + sgn = fmpz_sgn(den) < 0; + + if (sgn) + { + fmpz_neg(num, num); + fmpz_neg(den, den); + } + + fmpz_gcd(g, num, den); + res = fmpz_is_one(g); + + if (res) + { + NF_ELEM(a)->length = _fmpz_poly_get_n_adic(NF_ELEM_NUMREF(a), + lenf - 1, num, n); + + fmpz_mul(NF_ELEM_DENREF(a), NF_ELEM_DENREF(b), den); + + if (!fmpz_is_one(d)) + { + fmpz_mul(NF_ELEM_DENREF(a), NF_ELEM_DENREF(a), d); + fmpq_poly_canonicalise(NF_ELEM(a)); + } + + nf_elem_mul(sqr, a, a, nf); - nf_elem_mul(sqr, a, a, nf); + res = nf_elem_equal(sqr, b, nf); + } - res = nf_elem_equal(sqr, b, nf); + if (sgn) + { + fmpz_neg(num, num); + fmpz_neg(den, den); + } } + + fmpz_clear(q); + fmpz_clear(r); + fmpz_clear(s); + fmpz_clear(t); + fmpz_clear(g); + fmpz_clear(num); + fmpz_clear(den); } } fmpz_clear(m); _fmpz_vec_clear(mr, fac->num); - + _fmpz_vec_clear(modulus, fac->num); nf_elem_clear(sqr, nf); nbits = 1.1*nbits + 1; /* increase nbits and try again if sqrt not found */ @@ -586,11 +699,13 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) fmpz_clear(disc); fmpz_clear(n); + fmpz_clear(d); fmpz_clear(temp); fmpz_clear(z); _fmpz_vec_clear(r, fac->num); _fmpz_vec_clear(bz, NF_ELEM(b)->length); + _fmpz_vec_clear(bz1, NF_ELEM(b)->length); fmpz_factor_clear(fac); diff --git a/nf_elem/test/t-sqrt.c b/nf_elem/test/t-sqrt.c index 2a7f61745..2260b69a3 100644 --- a/nf_elem/test/t-sqrt.c +++ b/nf_elem/test/t-sqrt.c @@ -29,7 +29,218 @@ main(void) flint_randinit(state); - /* test sqrt(a^2) */ + /* Regression test */ + { + nf_t nf; + nf_elem_t a, b, c, d; + int is_square; + fmpq_poly_t f; + + /* f = x^3 - 1953*x^2 - x + 1 */ + fmpq_poly_init(f); + fmpq_poly_set_coeff_si(f, 0, 1); + fmpq_poly_set_coeff_si(f, 1, -1); + fmpq_poly_set_coeff_si(f, 2, -1953); + fmpq_poly_set_coeff_si(f, 3, 1); + + nf_init(nf, f); + + nf_elem_init(a, nf); + nf_elem_init(b, nf); + nf_elem_init(c, nf); + nf_elem_init(d, nf); + + /* a = -26977/6*x^2+40549/6 */ + fmpq_poly_set_coeff_si(NF_ELEM(a), 0, 40549); + fmpq_poly_set_coeff_si(NF_ELEM(a), 2, -26977); + fmpz_set_ui(NF_ELEM_DENREF(a), 6); + + 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); + + fmpq_poly_clear(f); + } + + /* Regression test */ + { + nf_t nf; + nf_elem_t a, b, c, d; + int is_square; + fmpq_poly_t f; + + /* f = x^3 + 44*x^2 - 3*x + 18 */ + fmpq_poly_init(f); + fmpq_poly_set_coeff_si(f, 0, 18); + fmpq_poly_set_coeff_si(f, 1, -3); + fmpq_poly_set_coeff_si(f, 2, 44); + fmpq_poly_set_coeff_si(f, 3, 1); + + nf_init(nf, f); + + nf_elem_init(a, nf); + nf_elem_init(b, nf); + nf_elem_init(c, nf); + nf_elem_init(d, nf); + + /* a = -1/80916*x^2 + 1/80916*x + 2/80916 */ + fmpq_poly_set_coeff_si(NF_ELEM(a), 0, 2); + fmpq_poly_set_coeff_si(NF_ELEM(a), 1, 1); + fmpq_poly_set_coeff_si(NF_ELEM(a), 2, -1); + fmpz_set_ui(NF_ELEM_DENREF(a), 80916); + + 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); + + fmpq_poly_clear(f); + } + + /* Regression test */ + { + nf_t nf; + nf_elem_t a, b, c, d; + int is_square; + fmpq_poly_t f; + + /* f = x^3 + 55*x^2 - 10*x - 1 */ + fmpq_poly_init(f); + fmpq_poly_set_coeff_si(f, 0, -1); + fmpq_poly_set_coeff_si(f, 1, -10); + fmpq_poly_set_coeff_si(f, 2, 55); + fmpq_poly_set_coeff_si(f, 3, 1); + + nf_init(nf, f); + + nf_elem_init(a, nf); + nf_elem_init(b, nf); + nf_elem_init(c, nf); + nf_elem_init(d, nf); + + /* a = 862/16*x^2 - 149/16*x - 13/16 */ + fmpq_poly_set_coeff_si(NF_ELEM(a), 0, -13); + fmpq_poly_set_coeff_si(NF_ELEM(a), 1, -149); + fmpq_poly_set_coeff_si(NF_ELEM(a), 2, 862); + fmpz_set_ui(NF_ELEM_DENREF(a), 16); + + 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); + + fmpq_poly_clear(f); + } + + /* Regression test */ + { + nf_t nf; + nf_elem_t a, b, c, d; + int is_square; + fmpq_poly_t f; + + /* f = x^8 - 8 */ + fmpq_poly_init(f); + fmpq_poly_set_coeff_si(f, 0, -8); + fmpq_poly_set_coeff_si(f, 8, 1); + + nf_init(nf, f); + + nf_elem_init(a, nf); + nf_elem_init(b, nf); + nf_elem_init(c, nf); + nf_elem_init(d, nf); + + /* a = 4050*x^6 */ + fmpq_poly_set_coeff_si(NF_ELEM(a), 6, 4050); + fmpz_set_ui(NF_ELEM_DENREF(a), 1); + + 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); + + fmpq_poly_clear(f); + } + + /* test sqrt(a^2) */ for (i = 0; i < 100 * antic_test_multiplier(); ) { nf_t nf; @@ -60,7 +271,6 @@ 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); From d864324b1f68c4cda621ccfe5d94c30f590ae41a Mon Sep 17 00:00:00 2001 From: William Hart Date: Tue, 29 Sep 2020 11:56:44 +0200 Subject: [PATCH 16/24] Remove memory leaks. --- nf_elem/sqrt.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/nf_elem/sqrt.c b/nf_elem/sqrt.c index 43086891d..637c0c878 100755 --- a/nf_elem/sqrt.c +++ b/nf_elem/sqrt.c @@ -32,7 +32,6 @@ * add is_square function * test non-squares * fix bug in fmpz_factor_trial #843 - * valgrind */ #define ROT(u,v,t) \ @@ -514,6 +513,7 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) { res = 0; nf_elem_zero(a, nf); + _fmpz_vec_clear(r, fac->num); _fmpz_vec_clear(modulus, fac->num); goto cleanup; } @@ -687,6 +687,7 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) } fmpz_clear(m); + _fmpz_vec_clear(r, fac->num); _fmpz_vec_clear(mr, fac->num); _fmpz_vec_clear(modulus, fac->num); nf_elem_clear(sqr, nf); @@ -703,7 +704,6 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) fmpz_clear(temp); fmpz_clear(z); - _fmpz_vec_clear(r, fac->num); _fmpz_vec_clear(bz, NF_ELEM(b)->length); _fmpz_vec_clear(bz1, NF_ELEM(b)->length); From 1d63916b72085aad563386d2a1c9e75d137166d6 Mon Sep 17 00:00:00 2001 From: William Hart Date: Tue, 29 Sep 2020 14:25:10 +0200 Subject: [PATCH 17/24] Test non-monic case. --- nf/init_randtest.c | 9 +++++ nf_elem.h | 3 ++ nf_elem/doc/nf_elem.txt | 9 ++++- nf_elem/randtest.c | 78 +++++++++++++++++++++++++++++++++++++++++ nf_elem/sqrt.c | 1 - nf_elem/test/t-sqrt.c | 70 ++++++++++++++++++++++++++++++++++-- 6 files changed, 166 insertions(+), 4 deletions(-) diff --git a/nf/init_randtest.c b/nf/init_randtest.c index 109218320..2159ea0ad 100644 --- a/nf/init_randtest.c +++ b/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); } diff --git a/nf_elem.h b/nf_elem.h index fc8c05335..55648bcbc 100644 --- a/nf_elem.h +++ b/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); diff --git a/nf_elem/doc/nf_elem.txt b/nf_elem/doc/nf_elem.txt index a032312b6..4b2270eb3 100644 --- a/nf_elem/doc/nf_elem.txt +++ b/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) diff --git a/nf_elem/randtest.c b/nf_elem/randtest.c index 028c78651..28d746351 100644 --- a/nf_elem/randtest.c +++ b/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) { diff --git a/nf_elem/sqrt.c b/nf_elem/sqrt.c index 637c0c878..be98418f5 100755 --- a/nf_elem/sqrt.c +++ b/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 diff --git a/nf_elem/test/t-sqrt.c b/nf_elem/test/t-sqrt.c index 2260b69a3..0c6eb42fd 100644 --- a/nf_elem/test/t-sqrt.c +++ b/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); From b17ea471903fd28c10bf4e98186c9569cd1805cb Mon Sep 17 00:00:00 2001 From: William Hart Date: Tue, 29 Sep 2020 15:11:09 +0200 Subject: [PATCH 18/24] Add is_square and test non-squares. --- nf_elem.h | 2 + nf_elem/doc/nf_elem.txt | 13 +++ nf_elem/is_square.c | 32 ++++++++ nf_elem/sqrt.c | 34 +++++--- nf_elem/test/t-is_square.c | 158 +++++++++++++++++++++++++++++++++++++ 5 files changed, 227 insertions(+), 12 deletions(-) create mode 100755 nf_elem/is_square.c create mode 100644 nf_elem/test/t-is_square.c diff --git a/nf_elem.h b/nf_elem.h index 55648bcbc..4e14a5fb0 100644 --- a/nf_elem.h +++ b/nf_elem.h @@ -941,6 +941,8 @@ 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 diff --git a/nf_elem/doc/nf_elem.txt b/nf_elem/doc/nf_elem.txt index 4b2270eb3..b258d7668 100644 --- a/nf_elem/doc/nf_elem.txt +++ b/nf_elem/doc/nf_elem.txt @@ -336,6 +336,19 @@ 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. + ******************************************************************************* Representation matrix diff --git a/nf_elem/is_square.c b/nf_elem/is_square.c new file mode 100755 index 000000000..677efa372 --- /dev/null +++ b/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 . + +=============================================================================*/ +/****************************************************************************** + + 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; +} diff --git a/nf_elem/sqrt.c b/nf_elem/sqrt.c index be98418f5..3dc4a81ec 100755 --- a/nf_elem/sqrt.c +++ b/nf_elem/sqrt.c @@ -22,14 +22,10 @@ TODO: * try to reuse information from previous failed attempt - * improve bounds - * 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) * Remove small squares from denominator before rationalising * Cache factorisation of f(n) on number field for future square roots - * add is_square function - * test non-squares * fix bug in fmpz_factor_trial #843 */ @@ -86,6 +82,9 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) fmpz_t r; int res = 1; + if (fmpz_sgn(bnum) < 0) + return 0; + fmpz_init(r); fmpz_sqrtrem(anum, r, bnum); @@ -166,14 +165,16 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) fmpq_mul(tq, r, r); fmpq_mul_2exp(tq2, s, 2); fmpq_sub(tq, tq, tq2); - fmpz_sqrtrem(fmpq_numref(s), t, fmpq_numref(tq)); - if (!fmpz_is_zero(t)) + if (fmpz_sgn(fmpq_numref(tq)) >= 0) + fmpz_sqrtrem(fmpq_numref(s), t, fmpq_numref(tq)); + if (fmpz_sgn(fmpq_numref(tq)) < 0 || !fmpz_is_zero(t)) { ret = 0; nf_elem_zero(a, nf); goto quadratic_cleanup; } + fmpz_sqrtrem(fmpq_denref(s), t, fmpq_denref(tq)); if (!fmpz_is_zero(t)) { @@ -195,8 +196,12 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) fmpq_div_fmpz(s, s, d); /* compute m and +/- n */ - fmpz_sqrtrem(t2, t, fmpq_numref(r)); - sq = fmpz_is_zero(t); + if (fmpz_sgn(fmpq_numref(r)) >= 0) + { + fmpz_sqrtrem(t2, t, fmpq_numref(r)); + sq = fmpz_is_zero(t); + } else + sq = 0; if (sq) { fmpz_sqrtrem(t3, t, fmpq_denref(r)); @@ -207,8 +212,12 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) fmpq_mul_fmpz(s, s, d); fmpq_div_fmpz(r, r, d); fmpq_swap(r, s); - fmpz_sqrtrem(t2, t, fmpq_numref(r)); - sq = fmpz_is_zero(t); + if (fmpz_sgn(fmpq_numref(r)) >= 0) + { + fmpz_sqrtrem(t2, t, fmpq_numref(r)); + sq = fmpz_is_zero(t); + } else + sq = 0; if (sq) { fmpz_sqrtrem(t3, t, fmpq_denref(r)); @@ -224,8 +233,9 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) } fmpz_swap(fmpq_numref(r), t2); fmpz_swap(fmpq_denref(r), t3); - fmpz_sqrtrem(fmpq_numref(s), t, fmpq_numref(s)); - if (!fmpz_is_zero(t)) + if (fmpz_sgn(fmpq_numref(s)) >= 0) + fmpz_sqrtrem(fmpq_numref(s), t, fmpq_numref(s)); + if (fmpz_sgn(fmpq_numref(s)) < 0 || !fmpz_is_zero(t)) { ret = 0; nf_elem_zero(a, nf); diff --git a/nf_elem/test/t-is_square.c b/nf_elem/test/t-is_square.c new file mode 100644 index 000000000..509801bff --- /dev/null +++ b/nf_elem/test/t-is_square.c @@ -0,0 +1,158 @@ +/*============================================================================= + + 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 . + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2020 William Hart + +******************************************************************************/ + +#include +#include "nf.h" +#include "nf_elem.h" + +int +main(void) +{ + int i, result; + flint_rand_t state; + + flint_printf("is_square...."); + fflush(stdout); + + flint_randinit(state); + + /* test a^2 is a square */ + for (i = 0; i < 100 * antic_test_multiplier(); ) + { + nf_t nf; + nf_elem_t a, b; + int 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_randtest_bounded(a, state, abits, nf); + + nf_elem_mul(b, a, a, nf); + + result = nf_elem_is_square(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"); + abort(); + } + + nf_elem_clear(a, nf); + nf_elem_clear(b, nf); + } + + nf_clear(nf); + } + + /* test non-squares */ + for (i = 0; i < 100 * antic_test_multiplier(); ) + { + nf_t nf; + nf_elem_t a, b, z; + int 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(z, nf); + + do { + nf_elem_randtest_bounded(z, state, abits, nf); + } while (nf_elem_is_square(z, nf)); + + do { + nf_elem_randtest_bounded(a, state, abits, nf); + } while (nf_elem_is_zero(a, nf)); + + nf_elem_mul(b, a, a, nf); + nf_elem_mul(b, b, z, nf); + + result = !nf_elem_is_square(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("z = "); nf_elem_print_pretty(z, nf, "x"); printf("\n"); + abort(); + } + + nf_elem_clear(a, nf); + nf_elem_clear(b, nf); + nf_elem_clear(z, nf); + } + + nf_clear(nf); + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return 0; +} From a65b06f1bdb5c08e7736d7442d5817b2fea50a5c Mon Sep 17 00:00:00 2001 From: William Hart Date: Tue, 29 Sep 2020 15:42:42 +0200 Subject: [PATCH 19/24] Rough bound improvement. --- nf_elem/sqrt.c | 9 ++++----- nf_elem/test/t-sqrt.c | 6 +++--- 2 files changed, 7 insertions(+), 8 deletions(-) diff --git a/nf_elem/sqrt.c b/nf_elem/sqrt.c index 3dc4a81ec..2d72b802a 100755 --- a/nf_elem/sqrt.c +++ b/nf_elem/sqrt.c @@ -21,9 +21,6 @@ /* TODO: - * try to reuse information from previous failed attempt - * 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) * Remove small squares from denominator before rationalising * Cache factorisation of f(n) on number field for future square roots * fix bug in fmpz_factor_trial #843 @@ -363,8 +360,10 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) flint_printf("Step 2\n"); #endif - bbits = FLINT_ABS(_fmpz_vec_max_bits(bz, lenb)); - nbits = (bbits + 1)/(20) + 2; + bbits = FLINT_ABS(_fmpz_vec_max_bits(bz1, lenb)); + bbits -= (lenf - 1)*(FLINT_ABS(fmpz_bits(fmpq_poly_numref(nf->pol) + lenf - 1)) - 1); + bbits = FLINT_MAX(bbits, 1); + nbits = (2*bbits + 1)/(lenf) + 2; /* Step 3: find a nbits bit prime such that z = f(n) is a product diff --git a/nf_elem/test/t-sqrt.c b/nf_elem/test/t-sqrt.c index 0c6eb42fd..a864a9b69 100644 --- a/nf_elem/test/t-sqrt.c +++ b/nf_elem/test/t-sqrt.c @@ -315,9 +315,9 @@ main(void) 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; + flen = n_randint(state, 15) + 2; + fbits = n_randint(state, 15) + 1; + abits = n_randint(state, 15) + 1; nf_init_randtest(nf, state, flen, fbits); From 017f6f140571c83f2952525203fbc5d98992c61d Mon Sep 17 00:00:00 2001 From: William Hart Date: Tue, 29 Sep 2020 16:25:10 +0200 Subject: [PATCH 20/24] Remove small squares from denominator. --- nf_elem/sqrt.c | 39 ++++++++++++++++++++++++++++++--------- nf_elem/test/t-sqrt.c | 2 +- 2 files changed, 31 insertions(+), 10 deletions(-) diff --git a/nf_elem/sqrt.c b/nf_elem/sqrt.c index 2d72b802a..3b6803de6 100755 --- a/nf_elem/sqrt.c +++ b/nf_elem/sqrt.c @@ -21,7 +21,6 @@ /* TODO: - * Remove small squares from denominator before rationalising * Cache factorisation of f(n) on number field for future square roots * fix bug in fmpz_factor_trial #843 */ @@ -295,11 +294,11 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) } else /* generic nf_elem */ { const slong lenb = NF_ELEM(b)->length, lenf = fmpq_poly_length(nf->pol); - slong nbits, bbits; + slong nbits, bbits, primes; fmpq_t bnorm; fmpz_factor_t fac; flint_rand_t state; - fmpz_t disc, z, temp, n, m, az, d; + fmpz_t disc, z, temp, n, m, az, d, bden; nf_elem_t sqr; slong i, j, k; fmpz * r, * mr, * bz, * bz1, * modulus; @@ -341,11 +340,34 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) } fmpz_init(temp); - + fmpz_init(bden); + + fmpz_set_ui(bden, 1); + /* get rid of denominator */ bz1 = _fmpz_vec_init(NF_ELEM(b)->length); bz = _fmpz_vec_init(NF_ELEM(b)->length); - _fmpz_vec_scalar_mul_fmpz(bz1, NF_ELEM_NUMREF(b), NF_ELEM(b)->length, NF_ELEM_DENREF(b)); + + fmpz_factor_init(fac); + nbits = FLINT_ABS(fmpz_bits(NF_ELEM_DENREF(b))); + primes = FLINT_MIN((slong) 4*log(nbits)*nbits + 2, 3512); + fmpz_factor_trial(fac, NF_ELEM_DENREF(b), primes); + + for (i = 0; i < fac->num; i++) + { + if (fac->exp[i] > 1) + { + fmpz_pow_ui(temp, fac->p + i, fac->exp[i]/2); + fmpz_mul(bden, bden, temp); + } + } + + fmpz_fdiv_q(temp, NF_ELEM_DENREF(b), bden); + fmpz_fdiv_q(temp, temp, bden); + + _fmpz_vec_scalar_mul_fmpz(bz1, NF_ELEM_NUMREF(b), NF_ELEM(b)->length, temp); + + fmpz_mul(bden, temp, bden); /* Step 2: compute number of bits for initial n @@ -400,9 +422,7 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) #endif while (!factored || fac->num > 14) /* no bound known for finding such a factorisation */ - { - slong primes; - + { fmpz_set_ui(d, 1); fmpz_factor_clear(fac); fmpz_factor_init(fac); @@ -608,7 +628,7 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) NF_ELEM(a)->length = _fmpz_poly_get_n_adic(NF_ELEM_NUMREF(a), lenf - 1, z, n); - fmpz_set(NF_ELEM_DENREF(a), NF_ELEM_DENREF(b)); + fmpz_set(NF_ELEM_DENREF(a), bden); if (!fmpz_is_one(d)) { @@ -711,6 +731,7 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) fmpz_clear(d); fmpz_clear(temp); fmpz_clear(z); + fmpz_clear(bden); _fmpz_vec_clear(bz, NF_ELEM(b)->length); _fmpz_vec_clear(bz1, NF_ELEM(b)->length); diff --git a/nf_elem/test/t-sqrt.c b/nf_elem/test/t-sqrt.c index a864a9b69..695027bd2 100644 --- a/nf_elem/test/t-sqrt.c +++ b/nf_elem/test/t-sqrt.c @@ -271,7 +271,7 @@ 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); From 9cb8267a2de9c30eb57bc8df5bba1144135e4002 Mon Sep 17 00:00:00 2001 From: William Hart Date: Fri, 2 Oct 2020 12:23:46 +0200 Subject: [PATCH 21/24] Add doc clarification suggested by Fredrik. --- nf_elem/doc/nf_elem.txt | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/nf_elem/doc/nf_elem.txt b/nf_elem/doc/nf_elem.txt index b258d7668..30ef70975 100644 --- a/nf_elem/doc/nf_elem.txt +++ b/nf_elem/doc/nf_elem.txt @@ -347,7 +347,8 @@ int nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) int nf_elem_is_square(const nf_elem_t b, const nf_t nf) - Return \code{1} if \code{b} is a square. + Return \code{1} if \code{b} is a square. This is currently not faster than + computing the square root. ******************************************************************************* From e8c959b989a560516182a716f3e78946ad8e480b Mon Sep 17 00:00:00 2001 From: William Hart Date: Fri, 2 Oct 2020 12:57:35 +0200 Subject: [PATCH 22/24] Make correction for non-monic case. --- nf_elem/sqrt.c | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/nf_elem/sqrt.c b/nf_elem/sqrt.c index 3b6803de6..af0011a04 100755 --- a/nf_elem/sqrt.c +++ b/nf_elem/sqrt.c @@ -406,6 +406,10 @@ 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); + /* denominator only has primes dividing the discriminant and gcd(lc(f), tc(f)) */ + fmpz_gcd(temp, fmpq_poly_numref(nf->pol) + 0, fmpq_poly_numref(nf->pol) + lenf - 1); + fmpz_mul(disc, disc, temp); + do /* continue increasing nbits until square root found */ { fmpz_t fac1; From 517290cf6c08668ba150f10d70e4fecc340a41b1 Mon Sep 17 00:00:00 2001 From: William Hart Date: Tue, 6 Oct 2020 00:07:48 +0200 Subject: [PATCH 23/24] Use proved bound for rational reconstruction and denominator and add another regression test. --- nf_elem/sqrt.c | 20 +++++++++---- nf_elem/test/t-sqrt.c | 67 +++++++++++++++++++++++++++++++++++++++---- 2 files changed, 76 insertions(+), 11 deletions(-) diff --git a/nf_elem/sqrt.c b/nf_elem/sqrt.c index af0011a04..8c851dc5f 100755 --- a/nf_elem/sqrt.c +++ b/nf_elem/sqrt.c @@ -298,7 +298,7 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) fmpq_t bnorm; fmpz_factor_t fac; flint_rand_t state; - fmpz_t disc, z, temp, n, m, az, d, bden; + fmpz_t disc, z, temp, n, m, az, d, bden, g, maxd; nf_elem_t sqr; slong i, j, k; fmpz * r, * mr, * bz, * bz1, * modulus; @@ -362,6 +362,8 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) } } + fmpz_factor_clear(fac); + fmpz_fdiv_q(temp, NF_ELEM_DENREF(b), bden); fmpz_fdiv_q(temp, temp, bden); @@ -402,14 +404,18 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) fmpz_set_ui(d, 1); fmpz_init(disc); + fmpz_init(g); + fmpz_init(maxd); flint_randinit(state); _fmpz_poly_discriminant(disc, fmpq_poly_numref(nf->pol), lenf); /* denominator only has primes dividing the discriminant and gcd(lc(f), tc(f)) */ - fmpz_gcd(temp, fmpq_poly_numref(nf->pol) + 0, fmpq_poly_numref(nf->pol) + lenf - 1); - fmpz_mul(disc, disc, temp); - + fmpz_gcd(g, fmpq_poly_numref(nf->pol) + 0, fmpq_poly_numref(nf->pol) + lenf - 1); + fmpz_pow_ui(temp, g, FLINT_MAX((lenf - 2)*(lenf - 3) - lenf/2, 0)); + fmpz_mul(maxd, disc, temp); + fmpz_mul(disc, disc, g); + do /* continue increasing nbits until square root found */ { fmpz_t fac1; @@ -666,7 +672,7 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) fmpz_set(r, m); fmpz_zero(s); fmpz_set(num, z); fmpz_one(den); - while (!res && !fmpz_is_zero(num)) /* try rational coeffs */ + while (!res && !fmpz_is_zero(num) && fmpz_cmpabs(den, maxd) <= 0) /* try rational coeffs */ { fmpz_fdiv_q(q, r, num); fmpz_mul(t, q, num); fmpz_sub(t, r, t); ROT(r, num, t); @@ -688,7 +694,7 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) NF_ELEM(a)->length = _fmpz_poly_get_n_adic(NF_ELEM_NUMREF(a), lenf - 1, num, n); - fmpz_mul(NF_ELEM_DENREF(a), NF_ELEM_DENREF(b), den); + fmpz_mul(NF_ELEM_DENREF(a), bden, den); if (!fmpz_is_one(d)) { @@ -730,6 +736,8 @@ int _nf_elem_sqrt(nf_elem_t a, const nf_elem_t b, const nf_t nf) cleanup: flint_randclear(state); fmpz_clear(disc); + fmpz_clear(g); + fmpz_clear(maxd); fmpz_clear(n); fmpz_clear(d); diff --git a/nf_elem/test/t-sqrt.c b/nf_elem/test/t-sqrt.c index 695027bd2..379b918a9 100644 --- a/nf_elem/test/t-sqrt.c +++ b/nf_elem/test/t-sqrt.c @@ -29,6 +29,62 @@ main(void) flint_randinit(state); + /* Regression test */ + { + nf_t nf; + nf_elem_t a, b, c, d; + int is_square; + fmpq_poly_t f; + + /* f = x^3 - 1953*x^2 - x + 1 */ + fmpq_poly_init(f); + fmpq_poly_set_coeff_si(f, 0, 7); + fmpq_poly_set_coeff_si(f, 1, -2); + fmpq_poly_set_coeff_si(f, 2, 1); + fmpq_poly_set_coeff_si(f, 3, -2); + fmpq_poly_set_coeff_si(f, 4, 7); + + nf_init(nf, f); + + nf_elem_init(a, nf); + nf_elem_init(b, nf); + nf_elem_init(c, nf); + nf_elem_init(d, nf); + + /* a = -26977/6*x^2+40549/6 */ + fmpq_poly_set_coeff_si(NF_ELEM(a), 0, 3); + fmpq_poly_set_coeff_si(NF_ELEM(a), 1, -97); + fmpq_poly_set_coeff_si(NF_ELEM(a), 2, 7); + fmpq_poly_set_coeff_si(NF_ELEM(a), 3, 30); + fmpz_set_ui(NF_ELEM_DENREF(a), 2); + + 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 1:\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); + + fmpq_poly_clear(f); + } + /* Regression test */ { nf_t nf; @@ -64,7 +120,7 @@ main(void) result = is_square && nf_elem_equal(d, b, nf); if (!result) { - printf("FAIL:\n"); + printf("FAIL 2:\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"); @@ -118,7 +174,7 @@ main(void) result = is_square && nf_elem_equal(d, b, nf); if (!result) { - printf("FAIL:\n"); + printf("FAIL 3:\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"); @@ -172,7 +228,7 @@ main(void) result = is_square && nf_elem_equal(d, b, nf); if (!result) { - printf("FAIL:\n"); + printf("FAIL 4:\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"); @@ -222,7 +278,7 @@ main(void) result = is_square && nf_elem_equal(d, b, nf); if (!result) { - printf("FAIL:\n"); + printf("FAIL 5:\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"); @@ -289,6 +345,7 @@ main(void) if (!result) { printf("FAIL:\n"); + printf("f = "); fmpq_poly_print_pretty(nf->pol, "x"); printf("\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"); @@ -320,7 +377,6 @@ main(void) abits = n_randint(state, 15) + 1; nf_init_randtest(nf, state, flen, fbits); - fmpz_poly_factor_init(fac); pol->coeffs = nf->pol->coeffs; @@ -354,6 +410,7 @@ main(void) if (!result) { printf("FAIL:\n"); + printf("f = "); fmpq_poly_print_pretty(nf->pol, "x"); printf("\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"); From 43c2c63f7446401ecf2d9c0c23627e19b943d26e Mon Sep 17 00:00:00 2001 From: wbhart Date: Mon, 7 Feb 2022 16:43:54 +0100 Subject: [PATCH 24/24] Update nf_elem/doc/nf_elem.txt MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: François Bobot --- nf_elem/doc/nf_elem.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nf_elem/doc/nf_elem.txt b/nf_elem/doc/nf_elem.txt index 30ef70975..422e8914e 100644 --- a/nf_elem/doc/nf_elem.txt +++ b/nf_elem/doc/nf_elem.txt @@ -338,7 +338,7 @@ void nf_elem_trace(fmpq_t res, const nf_elem_t a, const nf_t nf) 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 + If \code{b} 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)