Merge pull request #1834 from flintlib/qqbar

Implement qqbar roots of qqbar polynomials

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

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},

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

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