diff --git a/doc/source/qqbar.rst b/doc/source/qqbar.rst index b0a4b6b440..9b845f391e 100644 --- a/doc/source/qqbar.rst +++ b/doc/source/qqbar.rst @@ -609,6 +609,41 @@ Polynomial roots of *mat* and then call :func:`qqbar_roots_fmpz_poly` with the same flags. +.. function:: int _qqbar_roots_poly_squarefree(qqbar_ptr roots, qqbar_srcptr coeffs, slong len, slong deg_limit, slong bits_limit) + + Writes to the vector *roots* the `d` roots of the polynomial + with algebraic number coefficients *coeffs* of length *len* (`d = len - 1`). + + Given the polynomial `f = a_0 + \ldots + a_d x^d` with coefficients in + `\overline{\mathbb{Q}}`, we construct an annihilating polynomial with + coefficients in `\mathbb{Q}` as + `g = \prod (\tilde a_0 + \ldots + \tilde a_d x^d)` + where the product is taken over all + combinations of algebraic conjugates `\tilde a_k` of the input + coefficients. + The polynomial `g` is subsequently factored to find candidate + roots. + + The leading coefficient `a_d` must be nonzero and the polynomial `f` + polynomial must be squarefree. + To compute roots of a general polynomial which may have repeated roots, + it is necessary to perform a squarefree factorization before calling + this function. + An option is to call :func:`gr_poly_roots` with a ``qqbar`` context object, + which wraps this function + and takes care of the initial squarefree factorization. + + Since the product `g` can explode in size very quickly, the *deg_limit* + and *bits_limit* parameters allow bounding the degree and working precision. + The function returns 1 on success and 0 on failure indicating that + such a limit has been exceeded. + Setting nonpositive values for the limits removes the restrictions; + however, the function can still fail and return 0 in that case if `g` + exceeds machine size. + +Note: to compute algebraic number roots of polynomials of various other +types, use :func:`gr_poly_roots_other`. + Roots of unity and trigonometric functions ------------------------------------------------------------------------------- diff --git a/src/gr/qqbar.c b/src/gr/qqbar.c index 373bb3c92e..0630550d78 100644 --- a/src/gr/qqbar.c +++ b/src/gr/qqbar.c @@ -12,6 +12,7 @@ #include "fmpz_poly.h" #include "fmpz_poly_factor.h" #include "fmpq.h" +#include "fmpq_vec.h" #include "fmpzi.h" #include "fexpr.h" #include "qqbar.h" @@ -1162,6 +1163,81 @@ TRIG3(acsc_pi) /* todo: exploit when we know that the field is real */ +/* todo: quickly skip nonreal roots over the real algebraic numbers */ +int +_gr_qqbar_poly_roots(gr_vec_t roots, gr_vec_t mult, const gr_poly_t poly, int flags, gr_ctx_t ctx) +{ + int status; + gr_ctx_t ZZ, Rx; + gr_vec_t fac, exp; + gr_ptr c; + slong i; + + if (poly->length == 0) + return GR_DOMAIN; + + /* todo: fast numerical check to avoid an exact squarefree factorization */ + + gr_ctx_init_fmpz(ZZ); + gr_ctx_init_gr_poly(Rx, ctx); + + gr_vec_set_length(roots, 0, ctx); + gr_vec_set_length(mult, 0, ZZ); + + c = gr_heap_init(ctx); + gr_vec_init(fac, 0, Rx); + gr_vec_init(exp, 0, ZZ); + + status = gr_poly_factor_squarefree(c, fac, exp, poly, ctx); + + if (status == GR_SUCCESS) + { + slong deg2; + qqbar_ptr croots; + int success; + slong j; + + for (i = 0; i < fac->length; i++) + { + gr_poly_struct * fac_i = gr_vec_entry_ptr(fac, i, Rx); + fmpz * exp_i = gr_vec_entry_ptr(exp, i, ZZ); + + deg2 = fac_i->length - 1; + + croots = _qqbar_vec_init(deg2); + + success = _qqbar_roots_poly_squarefree(croots, fac_i->coeffs, deg2 + 1, QQBAR_CTX(ctx)->deg_limit, QQBAR_CTX(ctx)->bits_limit); + if (!success) + { + status = GR_UNABLE; + break; + } + + for (j = 0; j < deg2; j++) + { + if (QQBAR_CTX(ctx)->real_only && !qqbar_is_real(croots + j)) + continue; + + GR_MUST_SUCCEED(gr_vec_append(roots, croots + j, ctx)); + GR_MUST_SUCCEED(gr_vec_append(mult, exp_i, ZZ)); + } + + _qqbar_vec_clear(croots, deg2); + } + } + + /* todo: qqbar_cmp_root_order, but must sort exponents as well */ + + gr_vec_clear(fac, Rx); + gr_vec_clear(exp, ZZ); + gr_heap_clear(c, ctx); + + gr_ctx_clear(ZZ); + gr_ctx_clear(Rx); + + return status; +} + /* todo: quickly skip nonreal roots over the real algebraic numbers */ int _gr_qqbar_poly_roots_other(gr_vec_t roots, gr_vec_t mult, const gr_poly_t poly, gr_ctx_t other_ctx, int flags, gr_ctx_t ctx) @@ -1220,6 +1296,50 @@ _gr_qqbar_poly_roots_other(gr_vec_t roots, gr_vec_t mult, const gr_poly_t poly, return status; } + /* Convert to the integer case */ + if (other_ctx->which_ring == GR_CTX_FMPQ) + { + fmpz_poly_t f; + fmpz_t den; + gr_ctx_t ZZ; + int status = GR_SUCCESS; + + gr_ctx_init_fmpz(ZZ); + fmpz_init(den); + fmpz_poly_init2(f, poly->length); + _fmpz_poly_set_length(f, poly->length); + _fmpq_vec_get_fmpz_vec_fmpz(f->coeffs, den, poly->coeffs, poly->length); + + status = _gr_qqbar_poly_roots_other(roots, mult, (gr_poly_struct *) f, ZZ, flags, ctx); + + fmpz_poly_clear(f); + fmpz_clear(den); + gr_ctx_clear(ZZ); + return status; + } + + if (other_ctx->which_ring == GR_CTX_REAL_ALGEBRAIC_QQBAR || other_ctx->which_ring == GR_CTX_COMPLEX_ALGEBRAIC_QQBAR) + return _gr_qqbar_poly_roots(roots, mult, poly, flags, ctx); + + /* Allow anything else that we can plausibly convert to qqbar */ + /* Todo: prefer computing the squarefree factorization in the original ring; this should + generally be more efficient. */ + if (other_ctx->which_ring == GR_CTX_FMPZI || + other_ctx->which_ring == GR_CTX_REAL_ALGEBRAIC_CA || + other_ctx->which_ring == GR_CTX_COMPLEX_ALGEBRAIC_CA || + other_ctx->which_ring == GR_CTX_RR_CA || + other_ctx->which_ring == GR_CTX_CC_CA) + { + gr_poly_t f; + int status = GR_SUCCESS; + gr_poly_init(f, ctx); + status = gr_poly_set_gr_poly_other(f, poly, other_ctx, ctx); + if (status == GR_SUCCESS) + status = _gr_qqbar_poly_roots(roots, mult, f, flags, ctx); + gr_poly_clear(f, ctx); + return status; + } + return GR_UNABLE; } @@ -1369,6 +1489,7 @@ gr_method_tab_input _qqbar_methods_input[] = {GR_METHOD_ASEC_PI, (gr_funcptr) _gr_qqbar_asec_pi}, {GR_METHOD_ACSC_PI, (gr_funcptr) _gr_qqbar_acsc_pi}, + {GR_METHOD_POLY_ROOTS, (gr_funcptr) _gr_qqbar_poly_roots}, {GR_METHOD_POLY_ROOTS_OTHER, (gr_funcptr) _gr_qqbar_poly_roots_other}, {0, (gr_funcptr) NULL}, diff --git a/src/python/flint_ctypes.py b/src/python/flint_ctypes.py index 8b5c401d3b..d718dab122 100644 --- a/src/python/flint_ctypes.py +++ b/src/python/flint_ctypes.py @@ -5690,6 +5690,10 @@ def eigenvalues(self, domain=None): ([Root a = 5.37228 of a^2-5*a-2, Root a = -0.372281 of a^2-5*a-2], [1, 1]) >>> Mat(ZZ)([[1,2],[3,4]]).eigenvalues(domain=RR) ([[-0.3722813232690143 +/- 3.01e-17], [5.372281323269014 +/- 3.31e-16]], [1, 1]) + >>> Mat(QQbar)([[1, 0, QQbar.i()], [0, 0, 1], [1, 1, 1]]).eigenvalues() + ([Root a = 1.94721 + 0.604643*I of a^6-4*a^5+4*a^4+2*a^3-3*a^2+1, Root a = 0.654260 - 0.430857*I of a^6-4*a^5+4*a^4+2*a^3-3*a^2+1, Root a = -0.601467 - 0.173786*I of a^6-4*a^5+4*a^4+2*a^3-3*a^2+1], [1, 1, 1]) + >>> Mat(ZZi)([[1, 0, ZZi.i()], [0, 0, 1], [1, 1, 1]]).eigenvalues(domain=QQbar) + ([Root a = 1.94721 + 0.604643*I of a^6-4*a^5+4*a^4+2*a^3-3*a^2+1, Root a = 0.654260 - 0.430857*I of a^6-4*a^5+4*a^4+2*a^3-3*a^2+1, Root a = -0.601467 - 0.173786*I of a^6-4*a^5+4*a^4+2*a^3-3*a^2+1], [1, 1, 1]) The matrix must be square: @@ -7715,6 +7719,21 @@ def test_set_str(): assert RRx("1 +/- 0") == RR(1) +def test_qqbar_roots(): + for R in [ZZ, QQ, ZZi, QQbar, AA, QQbar_ca, AA_ca, RR_ca, CC_ca]: + Rx = PolynomialRing(R) + assert Rx([-2,0,1]).roots(domain=AA) == ([AA(2).sqrt(), -AA(2).sqrt()], [1, 1]) + assert Rx([2,0,1]).roots(domain=AA) == ([], []) + assert Rx([2,0,1]).roots(domain=QQbar) == ([QQbar(-2).sqrt(), -QQbar(-2).sqrt()], [1, 1]) + assert (Rx([-2,0,1]) ** 2).roots(domain=AA) == ([AA(2).sqrt(), -AA(2).sqrt()], [2, 2]) + Rx = PolynomialRing(QQbar, "x") + x = Rx.gen() + g = -QQbar(3).sqrt() + x + f = 2 + QQbar(2).sqrt()*x + x**2 + h = g**2 * f + ((r1, r2, r3), (e1, e2, e3)) = h.roots(domain=QQbar) + assert (x-r1)**e1 * (x-r2)**e2 * (x-r3)**e3 == h + def test_ca_notebook_examples(): # algebraic number identity NumberI = fexpr("NumberI") diff --git a/src/qqbar.h b/src/qqbar.h index b31027f90c..a0b78ba165 100644 --- a/src/qqbar.h +++ b/src/qqbar.h @@ -373,8 +373,8 @@ int qqbar_evaluate_fmpz_mpoly(qqbar_t res, const fmpz_mpoly_t f, qqbar_srcptr x, #define QQBAR_ROOTS_UNSORTED 2 void qqbar_roots_fmpz_poly(qqbar_ptr res, const fmpz_poly_t poly, int flags); - void qqbar_roots_fmpq_poly(qqbar_ptr res, const fmpq_poly_t poly, int flags); +int _qqbar_roots_poly_squarefree(qqbar_ptr roots, qqbar_srcptr coeffs, slong len, slong deg_limit, slong bits_limit); void qqbar_eigenvalues_fmpz_mat(qqbar_ptr res, const fmpz_mat_t mat, int flags); diff --git a/src/qqbar/roots_poly_squarefree.c b/src/qqbar/roots_poly_squarefree.c new file mode 100644 index 0000000000..02aed64f80 --- /dev/null +++ b/src/qqbar/roots_poly_squarefree.c @@ -0,0 +1,344 @@ +/* + Copyright (C) 2024 Fredrik Johansson + + This file is part of FLINT. + + FLINT 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 3 of the License, or + (at your option) any later version. See . +*/ + +#include +#include "qqbar.h" +#include "acb_poly.h" +#include "arb_fmpz_poly.h" +#include "fmpz_poly_factor.h" + +/* +We convert an annihilating polynomial over QQbar to an annihilating polynomial +over Z by multiplying out all combinations of conjugates of the coefficients. + +TODO: + +* Consider actually allowing squareful polynomials. In general, an initial + squarefree factorization should be more efficient, but we *could* support + squareful polynomials also here with a few adjustments. + +* Consider using the quadratic formula for quadratics. + +* Compute tighter denominator bounds and/or attempt to rescale things to get + smaller coefficients. + +* Exploit when some coefficients lie in a small common number field to attempt + to reduce the degree. + +* Is it ever faster to verify roots with an exact computation instead of + numerics? + +* Guess some low-degree roots with LLL and remove them before running the full + factoring algorithm. + +* Compute tight precision estimates. + +*/ + +#define VERBOSE 0 + +/* Binary splitting product over all combinations of conjugates. */ +static void +bsprod(acb_poly_t res, slong a, slong b, slong deg, qqbar_srcptr coeffs, acb_ptr * conj_roots, slong prec) +{ + if (b - a == 1) + { + slong k, j; + acb_poly_fit_length(res, deg + 1); + _acb_poly_set_length(res, deg + 1); + k = a; + for (j = 0; j <= deg; j++) + { + acb_set(res->coeffs + j, conj_roots[j] + (k % qqbar_degree(coeffs + j))); + k /= qqbar_degree(coeffs + j); + } + } + else + { + acb_poly_t f, g; + + acb_poly_init(f); + acb_poly_init(g); + + bsprod(f, a, a + (b - a) / 2, deg, coeffs, conj_roots, prec); + bsprod(g, a + (b - a) / 2, b, deg, coeffs, conj_roots, prec); + acb_poly_mul(res, f, g, prec); + + acb_poly_clear(f); + acb_poly_clear(g); + } +} + +int +_qqbar_roots_poly_squarefree(qqbar_ptr roots, qqbar_srcptr coeffs, slong len, slong deg_limit, slong bits_limit) +{ + slong deg, deg_bound, conjugations, i; + slong prec, initial_prec; + acb_ptr * conj_roots; + acb_poly_t cpoly, ct; + fmpz_poly_t rpoly; + fmpz_t den, den_bound; + int success = 1; + + deg = len - 1; + + if (deg_limit < 0) + deg_limit = WORD_MAX; + + if (bits_limit < 0) + bits_limit = WORD_MAX / 8; + + if (deg <= 0) + return 1; + + FLINT_ASSERT(!qqbar_is_zero(coeffs + deg)); + + if (deg == 1) + { + /* todo: check limits */ + qqbar_div(roots, coeffs, coeffs + 1); + qqbar_neg(roots, roots); + return 1; + } + + conjugations = 1; + for (i = 0; i <= deg; i++) + { + ulong hi, lo; + umul_ppmm(hi, lo, conjugations, qqbar_degree(coeffs + i)); + if (hi != 0 || lo > deg_limit) + return 0; + conjugations = lo; + } + + fmpz_init(den_bound); + fmpz_init(den); + + fmpz_one(den_bound); + for (i = 0; i <= deg; i++) + { + qqbar_denominator(den, coeffs + i); + fmpz_lcm(den_bound, den_bound, den); + } + fmpz_pow_ui(den_bound, den_bound, conjugations); + + /* We have rational coefficients; run the specialized algorithm. */ + if (conjugations == 1) + { + fmpz_poly_t t; + fmpz_poly_init(t); + + fmpz_poly_fit_length(t, deg + 1); + _fmpz_poly_set_length(t, deg + 1); + + for (i = 0; i <= deg; i++) + { + fmpz_divexact(t->coeffs + i, den_bound, QQBAR_POLY(coeffs + i)->coeffs + 1); + fmpz_mul(t->coeffs + i, t->coeffs + i, QQBAR_POLY(coeffs + i)->coeffs); + fmpz_neg(t->coeffs + i, t->coeffs + i); + } + + qqbar_roots_fmpz_poly(roots, t, 0); + + fmpz_poly_clear(t); + fmpz_clear(den); + fmpz_clear(den_bound); + return 1; + } + + deg_bound = conjugations * deg; + if (deg_bound > deg_limit) + { + fmpz_clear(den); + fmpz_clear(den_bound); + return 0; + } + + conj_roots = flint_malloc(sizeof(acb_ptr) * (deg + 1)); + + for (i = 0; i <= deg; i++) + conj_roots[i] = _acb_vec_init(qqbar_degree(coeffs + i)); + + acb_poly_init(cpoly); + acb_poly_init(ct); + fmpz_poly_init(rpoly); + + fmpz_one(den_bound); + for (i = 0; i <= deg; i++) + { + qqbar_denominator(den, coeffs + i); + fmpz_lcm(den_bound, den_bound, den); + } + fmpz_pow_ui(den_bound, den_bound, conjugations); + +#if VERBOSE + flint_printf("deg_bound: %wd\n", deg_bound); + flint_printf("den_bound: %wd bits\n", fmpz_bits(den_bound)); +#endif + + initial_prec = 64 + fmpz_bits(den_bound) + 2 * FLINT_BIT_COUNT(conjugations); + initial_prec = FLINT_MAX(initial_prec, QQBAR_DEFAULT_PREC); + + for (prec = initial_prec; ; prec *= 2) + { +#if VERBOSE + flint_printf("prec = %wd\n", prec); +#endif + + for (i = 0; i <= deg; i++) + arb_fmpz_poly_complex_roots(conj_roots[i], QQBAR_POLY(coeffs + i), 0, prec); + + bsprod(cpoly, 0, conjugations, deg, coeffs, conj_roots, prec); + _acb_vec_scalar_mul_fmpz(cpoly->coeffs, cpoly->coeffs, cpoly->length, den_bound, prec); + +#if VERBOSE + acb_poly_printd(cpoly, 20); flint_printf("\n\n"); +#endif + + if (acb_poly_get_unique_fmpz_poly(rpoly, cpoly)) + break; + + /* todo: reuse enclosures (refine numerical roots) */ + /* todo: compute next precision from coefficient magnitude */ + + if (prec > bits_limit) + + + if (prec > 1000000) + flint_abort(); + } + +#if VERBOSE + flint_printf("rpoly: "); fmpz_poly_print(rpoly); flint_printf("\n\n"); +#endif + + slong found; + slong num_candidates = fmpz_poly_degree(rpoly); + slong num_candidates2; + + fmpz_poly_factor_t fac; + fmpz_poly_factor_init(fac); + fmpz_poly_factor(fac, rpoly); + + acb_t y; + acb_init(y); + + qqbar_ptr candidates; + slong * c_indices; + slong ci; + + candidates = _qqbar_vec_init(num_candidates); + c_indices = flint_malloc(sizeof(slong) * num_candidates); + + ci = 0; + for (i = 0; i < fac->num; i++) + { +#if VERBOSE + flint_printf("factor: "); + fmpz_poly_print(fac->p + i); flint_printf("\n"); +#endif + + qqbar_roots_fmpz_poly(candidates + ci, fac->p + i, QQBAR_ROOTS_IRREDUCIBLE | QQBAR_ROOTS_UNSORTED); + ci += fmpz_poly_degree(fac->p + i); + } + + num_candidates2 = ci; + +#if VERBOSE + for (ci = 0; ci < num_candidates2; ci++) + { + flint_printf("candidate %wd: ", ci); + qqbar_print(candidates + ci); flint_printf("\n"); + } +#endif + + initial_prec = QQBAR_DEFAULT_PREC; + for (prec = initial_prec; ; prec *= 2) + { + /* Check numerically for solutions to f(x) = 0. We will necessarily have + found >= deg. If found == deg, we have certainly isolated the + correct roots. Otherwise, we retry with higher precision. + This way we do not need exact qqbar arithmetic to validate f(x) = 0. */ + found = 0; + + acb_poly_fit_length(cpoly, deg + 1); + _acb_poly_set_length(cpoly, deg + 1); + for (i = 0; i <= deg; i++) + { + /* Todo: copy the original enclosures and reuse for refinement. */ + qqbar_get_acb(cpoly->coeffs + i, coeffs + i, prec); + } + + for (ci = 0; ci < num_candidates2; ci++) + { + /* Refine the candidate to the new precision */ + if (prec != initial_prec) + _qqbar_enclosure_raw(QQBAR_ENCLOSURE(candidates + ci), QQBAR_POLY(candidates + ci), QQBAR_ENCLOSURE(candidates + ci), prec); + + acb_poly_evaluate(y, cpoly, QQBAR_ENCLOSURE(candidates + ci), prec); + + if (acb_contains_zero(y)) + { + /* -- restart with higher precision */ + if (found == deg) + { + found = deg + 1; + break; + } + + c_indices[found] = ci; + found++; + } + } + +#if VERBOSE + flint_printf("prec %wd : found %wd / deg %wd\n", prec, found, deg); +#endif + + FLINT_ASSERT(found >= deg); + + if (found == deg) + { + for (i = 0; i < deg; i++) + qqbar_set(roots + i, candidates + c_indices[i]); + + break; + } + + if (prec > bits_limit) + { + success = 0; + break; + } + } + + _qqbar_vec_clear(candidates, fmpz_poly_degree(rpoly)); + flint_free(c_indices); + + fmpz_poly_factor_clear(fac); + acb_clear(y); + + for (i = 0; i <= deg; i++) + _acb_vec_clear(conj_roots[i], qqbar_degree(coeffs + i)); + + acb_poly_clear(cpoly); + acb_poly_clear(ct); + fmpz_poly_clear(rpoly); + fmpz_clear(den); + fmpz_clear(den_bound); + flint_free(conj_roots); + + /* This could be optional, but it's probably cheap compared to the actual computations... */ + qsort(roots, deg, sizeof(qqbar_struct), (int (*)(const void *, const void *)) qqbar_cmp_root_order); + + return success; +} + diff --git a/src/qqbar/test/main.c b/src/qqbar/test/main.c index 2173a2a3b2..a7aea93e31 100644 --- a/src/qqbar/test/main.c +++ b/src/qqbar/test/main.c @@ -61,6 +61,7 @@ #include "t-re_im.c" #include "t-root_of_unity.c" #include "t-roots_fmpz_poly.c" +#include "t-roots_poly_squarefree.c" #include "t-root_ui.c" #include "t-sec_pi.c" #include "t-set_d.c" @@ -121,6 +122,7 @@ test_struct tests[] = TEST_FUNCTION(qqbar_re_im), TEST_FUNCTION(qqbar_root_of_unity), TEST_FUNCTION(qqbar_roots_fmpz_poly), + TEST_FUNCTION(qqbar_roots_poly_squarefree), TEST_FUNCTION(qqbar_root_ui), TEST_FUNCTION(qqbar_sec_pi), TEST_FUNCTION(qqbar_set_d), diff --git a/src/qqbar/test/t-roots_poly_squarefree.c b/src/qqbar/test/t-roots_poly_squarefree.c new file mode 100644 index 0000000000..49dcba3a17 --- /dev/null +++ b/src/qqbar/test/t-roots_poly_squarefree.c @@ -0,0 +1,147 @@ +/* + Copyright (C) 2024 Fredrik Johansson + + This file is part of FLINT. + + FLINT 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 3 of the License, or + (at your option) any later version. See . +*/ + +#include "test_helpers.h" +#include "qqbar.h" +#include "gr.h" +#include "gr_poly.h" + +TEST_FUNCTION_START(qqbar_roots_poly_squarefree, state) +{ + slong iter; + slong count_success = 0; + + slong deg_limit = 200; + slong bits_limit = 1000; + + for (iter = 0; iter < 100 * flint_test_multiplier(); iter++) + { + gr_poly_t poly, fac; + gr_ctx_t ctx; + qqbar_ptr roots, computed_roots; + int success, status; + qqbar_t c; + slong deg, i, j; + + deg = n_randint(state, 4); + + roots = _qqbar_vec_init(deg); + computed_roots = _qqbar_vec_init(deg); + qqbar_init(c); + + gr_ctx_init_complex_qqbar(ctx); + _gr_ctx_qqbar_set_limits(ctx, deg_limit, bits_limit); + + gr_poly_init(poly, ctx); + gr_poly_init(fac, ctx); + + for (i = 0; i < deg; i++) + { + if (n_randint(state, 10) == 0) + qqbar_randtest(roots + i, state, 4, 4); + else + qqbar_randtest(roots + i, state, 2, 4); + + for (j = 0; j < i; j++) + { + if (qqbar_equal(roots + i, roots + j)) + { + i--; + break; + } + } + } + + if (n_randint(state, 5)) + { + qqbar_one(c); + } + else + { + do { + qqbar_randtest(c, state, 2, 4); + } while (qqbar_is_zero(c)); + } + + status = GR_SUCCESS; + status |= gr_poly_set_scalar(poly, c, ctx); + + for (i = 0; i < deg; i++) + { + status |= gr_poly_set_scalar(fac, roots + i, ctx); + status |= gr_poly_neg(fac, fac, ctx); + status |= gr_poly_set_coeff_si(fac, 1, 1, ctx); + status |= gr_poly_mul(poly, poly, fac, ctx); + } + + if (status == GR_SUCCESS) + { + success = _qqbar_roots_poly_squarefree(computed_roots, poly->coeffs, poly->length, deg_limit, bits_limit); + count_success += success; + + if (success) + { + for (i = 0; i < deg; i++) + { + int found = 0; + + for (j = 0; j < deg; j++) + { + if (qqbar_equal(roots + i, computed_roots + j)) + { + found = 1; + break; + } + } + + if (!found) + { + flint_printf("FAIL\n\n"); + flint_printf("deg = %wd\n", deg); + flint_printf("roots:\n"); + + for (j = 0; j < deg; j++) + { + qqbar_print(roots + i); + flint_printf("\n"); + } + + flint_printf("computed roots:\n"); + for (j = 0; j < deg; j++) + { + qqbar_print(roots + i); + flint_printf("\n"); + } + + flint_abort(); + } + } + } + +/* flint_printf("total: %wd / %wd\n", count_success, iter); */ + } + + gr_poly_clear(poly, ctx); + gr_poly_clear(fac, ctx); + gr_ctx_clear(ctx); + _qqbar_vec_clear(roots, deg); + _qqbar_vec_clear(computed_roots, deg); + qqbar_clear(c); + } + + if (iter > 100 && count_success < 0.1 * iter) + { + flint_printf("FAIL: only %wd / %wd succeeded\n", count_success, iter); + flint_abort(); + } + + TEST_FUNCTION_END(state); +}