add arb_poly_product_roots_complex; add test code and slight optimization for arb/acb_poly_product_roots

acb_poly/product_roots.c

@@ -31,6 +31,17 @@ _acb_poly_product_roots(acb_ptr poly, acb_srcptr xs, slong n, slong prec)
         acb_neg(poly + 1, poly + 1);
         acb_one(poly + 2);
     }
+    else if (n == 3)
+    {
+        acb_mul(poly + 1, xs, xs + 1, prec);
+        acb_mul(poly, poly + 1, xs + 2, prec);
+        acb_neg(poly, poly);
+        acb_add(poly + 2, xs, xs + 1, prec);
+        acb_addmul(poly + 1, poly + 2, xs + 2, prec);
+        acb_add(poly + 2, poly + 2, xs + 2, prec);
+        acb_neg(poly + 2, poly + 2);
+        acb_one(poly + 3);
+    }
     else
     {
         const slong m = (n + 1) / 2;

acb_poly/test/t-product_roots.c

@@ -0,0 +1,81 @@
+/*
+    Copyright (C) 2017 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "acb_poly.h"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("product_roots....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 10000 * arb_test_multiplier(); iter++)
+    {
+        slong prec, len, m, i;
+        acb_poly_t a, b, c, d;
+        acb_ptr r;
+
+        prec = 2 + n_randint(state, 200);
+        len = n_randint(state, 12);
+        if (len > 0)
+            m = n_randint(state, len);
+        else
+            m = 0;
+
+        acb_poly_init(a);
+        acb_poly_init(b);
+        acb_poly_init(c);
+        acb_poly_init(d);
+        r = _acb_vec_init(len);
+
+        for (i = 0; i < len; i++)
+            acb_randtest(r + i, state, 1 + n_randint(state, 200), 10);
+
+        acb_poly_product_roots(a, r, len, prec);
+        acb_poly_product_roots(b, r, m, prec);
+        acb_poly_product_roots(c, r + m, len - m, prec);
+        acb_poly_mul(d, b, c, prec);
+
+        if (!acb_poly_overlaps(a, d))
+        {
+            flint_printf("FAIL\n\n");
+
+            flint_printf("len = %wd, m = %wd\n\n", len, m);
+            for (i = 0; i < len; i++)
+            {
+                acb_printd(r + i, 15);
+                flint_printf("\n");
+            }
+
+            flint_printf("a = "); acb_poly_printd(a, 15); flint_printf("\n\n");
+            flint_printf("b = "); acb_poly_printd(b, 15); flint_printf("\n\n");
+            flint_printf("c = "); acb_poly_printd(c, 15); flint_printf("\n\n");
+            flint_printf("d = "); acb_poly_printd(d, 15); flint_printf("\n\n");
+
+            flint_abort();
+        }
+
+        acb_poly_clear(a);
+        acb_poly_clear(b);
+        acb_poly_clear(c);
+        acb_poly_clear(d);
+        _acb_vec_clear(r, len);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}

arb_poly.h

@@ -337,6 +337,12 @@ void _arb_poly_product_roots(arb_ptr poly, arb_srcptr xs, slong n, slong prec);
 
 void arb_poly_product_roots(arb_poly_t poly, arb_srcptr xs, slong n, slong prec);
 
+void _arb_poly_product_roots_complex(arb_ptr poly, arb_srcptr r, slong rn,
+    acb_srcptr c, slong cn, slong prec);
+
+void arb_poly_product_roots_complex(arb_poly_t poly,
+        arb_srcptr r, slong rn, acb_srcptr c, slong cn, slong prec);
+
 arb_ptr * _arb_poly_tree_alloc(slong len);
 
 void _arb_poly_tree_free(arb_ptr * tree, slong len);

arb_poly/product_roots.c

@@ -31,6 +31,17 @@ _arb_poly_product_roots(arb_ptr poly, arb_srcptr xs, slong n, slong prec)
         arb_neg(poly + 1, poly + 1);
         arb_one(poly + 2);
     }
+    else if (n == 3)
+    {
+        arb_mul(poly + 1, xs, xs + 1, prec);
+        arb_mul(poly, poly + 1, xs + 2, prec);
+        arb_neg(poly, poly);
+        arb_add(poly + 2, xs, xs + 1, prec);
+        arb_addmul(poly + 1, poly + 2, xs + 2, prec);
+        arb_add(poly + 2, poly + 2, xs + 2, prec);
+        arb_neg(poly + 2, poly + 2);
+        arb_one(poly + 3);
+    }
     else
     {
         const slong m = (n + 1) / 2;

arb_poly/product_roots_complex.c

@@ -0,0 +1,94 @@
+/*
+    Copyright (C) 2017 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "arb_poly.h"
+
+void
+_arb_poly_product_roots_complex(arb_ptr poly, arb_srcptr r, slong rn,
+    acb_srcptr c, slong cn, slong prec)
+{
+    if (rn == 0 && cn == 0)
+    {
+        arb_one(poly);
+    }
+    else if (rn == 1 && cn == 0)
+    {
+        arb_neg(poly, r);
+        arb_one(poly + 1);
+    }
+    else if (rn == 2 && cn == 0)
+    {
+        arb_mul(poly, r + 0, r + 1, prec);
+        arb_add(poly + 1, r + 0, r + 1, prec);
+        arb_neg(poly + 1, poly + 1);
+        arb_one(poly + 2);
+    }
+    else if (rn == 3 && cn == 0)
+    {
+        arb_mul(poly + 1, r, r + 1, prec);
+        arb_mul(poly, poly + 1, r + 2, prec);
+        arb_neg(poly, poly);
+        arb_add(poly + 2, r, r + 1, prec);
+        arb_addmul(poly + 1, poly + 2, r + 2, prec);
+        arb_add(poly + 2, poly + 2, r + 2, prec);
+        arb_neg(poly + 2, poly + 2);
+        arb_one(poly + 3);
+    }
+    else if (rn == 0 && cn == 1)
+    {
+        arb_mul(poly, acb_realref(c), acb_realref(c), prec);
+        arb_addmul(poly, acb_imagref(c), acb_imagref(c), prec);
+        arb_mul_2exp_si(poly + 1, acb_realref(c), 1);
+        arb_neg(poly + 1, poly + 1);
+        arb_one(poly + 2);
+    }
+    else if (rn == 1 && cn == 1)
+    {
+        arb_mul(poly + 1, acb_realref(c), acb_realref(c), prec);
+        arb_addmul(poly + 1, acb_imagref(c), acb_imagref(c), prec);
+        arb_mul(poly, poly + 1, r, prec);
+        arb_neg(poly, poly);
+        arb_mul_2exp_si(poly + 2, acb_realref(c), 1);
+        arb_addmul(poly + 1, poly + 2, r, prec);
+        arb_add(poly + 2, poly + 2, r, prec);
+        arb_neg(poly + 2, poly + 2);
+        arb_one(poly + 3);
+    }
+    else
+    {
+        slong rm, cm, rm2, cm2;
+        arb_ptr tmp, tmp2;
+
+        rm = (rn + 1) / 2;
+        cm = cn / 2;
+        rm2 = rn - rm;
+        cm2 = cn - cm;
+
+        tmp = _arb_vec_init(rn + 2 * cn + 2);
+        tmp2 = tmp + rm + (2 * cm) + 1;
+
+        _arb_poly_product_roots_complex(tmp, r, rm, c, cm, prec);
+        _arb_poly_product_roots_complex(tmp2, r + rm, rm2, c + cm, cm2, prec);
+        _arb_poly_mul_monic(poly, tmp, rm + 2 * cm + 1, tmp2, rm2 + 2 * cm2 + 1, prec);
+
+        _arb_vec_clear(tmp, rn + 2 * cn + 2);
+    }
+}
+
+void
+arb_poly_product_roots_complex(arb_poly_t poly,
+        arb_srcptr r, slong rn, acb_srcptr c, slong cn, slong prec)
+{
+    arb_poly_fit_length(poly, rn + 2 * cn + 1);
+    _arb_poly_product_roots_complex(poly->coeffs, r, rn, c, cn, prec);
+    _arb_poly_set_length(poly, rn + 2 * cn + 1);
+}
+

arb_poly/test/t-product_roots.c

@@ -0,0 +1,81 @@
+/*
+    Copyright (C) 2017 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "arb_poly.h"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("product_roots....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 3000 * arb_test_multiplier(); iter++)
+    {
+        slong prec, len, m, i;
+        arb_poly_t a, b, c, d;
+        arb_ptr r;
+
+        prec = 2 + n_randint(state, 200);
+        len = n_randint(state, 12);
+        if (len > 0)
+            m = n_randint(state, len);
+        else
+            m = 0;
+
+        arb_poly_init(a);
+        arb_poly_init(b);
+        arb_poly_init(c);
+        arb_poly_init(d);
+        r = _arb_vec_init(len);
+
+        for (i = 0; i < len; i++)
+            arb_randtest(r + i, state, 1 + n_randint(state, 200), 10);
+
+        arb_poly_product_roots(a, r, len, prec);
+        arb_poly_product_roots(b, r, m, prec);
+        arb_poly_product_roots(c, r + m, len - m, prec);
+        arb_poly_mul(d, b, c, prec);
+
+        if (!arb_poly_overlaps(a, d))
+        {
+            flint_printf("FAIL\n\n");
+
+            flint_printf("len = %wd, m = %wd\n\n", len, m);
+            for (i = 0; i < len; i++)
+            {
+                arb_printd(r + i, 15);
+                flint_printf("\n");
+            }
+
+            flint_printf("a = "); arb_poly_printd(a, 15); flint_printf("\n\n");
+            flint_printf("b = "); arb_poly_printd(b, 15); flint_printf("\n\n");
+            flint_printf("c = "); arb_poly_printd(c, 15); flint_printf("\n\n");
+            flint_printf("d = "); arb_poly_printd(d, 15); flint_printf("\n\n");
+
+            flint_abort();
+        }
+
+        arb_poly_clear(a);
+        arb_poly_clear(b);
+        arb_poly_clear(c);
+        arb_poly_clear(d);
+        _arb_vec_clear(r, len);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}

arb_poly/test/t-product_roots_complex.c

@@ -0,0 +1,97 @@
+/*
+    Copyright (C) 2017 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "arb_poly.h"
+#include "acb_poly.h"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("product_roots_complex....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 2000 * arb_test_multiplier(); iter++)
+    {
+        slong prec, len1, len2, i;
+        arb_poly_t a;
+        acb_poly_t b, c;
+        arb_ptr r;
+        acb_ptr r2, r3;
+
+        prec = 2 + n_randint(state, 200);
+        len1 = n_randint(state, 12);
+        len2 = n_randint(state, 8);
+
+        arb_poly_init(a);
+        acb_poly_init(b);
+        acb_poly_init(c);
+
+        r = _arb_vec_init(len1);
+        r2 = _acb_vec_init(len2);
+        r3 = _acb_vec_init(len1 + 2 * len2);
+
+        for (i = 0; i < len1; i++)
+            arb_randtest(r + i, state, 1 + n_randint(state, 200), 3);
+        for (i = 0; i < len2; i++)
+            acb_randtest(r2 + i, state, 1 + n_randint(state, 200), 3);
+
+        for (i = 0; i < len1; i++)
+            acb_set_arb(r3 + i, r + i);
+        for (i = 0; i < len2; i++)
+            acb_set(r3 + len1 + i, r2 + i);
+        for (i = 0; i < len2; i++)
+            acb_conj(r3 + len1 + len2 + i, r2 + i);
+
+        arb_poly_product_roots_complex(a, r, len1, r2, len2, prec);
+        acb_poly_product_roots(b, r3, len1 + 2 * len2, prec);
+
+        acb_poly_set_arb_poly(c, a);
+
+        if (!acb_poly_overlaps(b, c))
+        {
+            flint_printf("FAIL\n\n");
+
+            flint_printf("len1 = %wd, len2 = %wd\n\n", len1, len2);
+            for (i = 0; i < len1; i++)
+            {
+                arb_printd(r + i, 15);
+                flint_printf("\n");
+            }
+            for (i = 0; i < len2; i++)
+            {
+                acb_printd(r2 + i, 15);
+                flint_printf("\n");
+            }
+
+            flint_printf("a = "); arb_poly_printd(a, 15); flint_printf("\n\n");
+            flint_printf("b = "); acb_poly_printd(b, 15); flint_printf("\n\n");
+            flint_printf("c = "); acb_poly_printd(c, 15); flint_printf("\n\n");
+
+            flint_abort();
+        }
+
+        arb_poly_clear(a);
+        acb_poly_clear(b);
+        acb_poly_clear(c);
+        _arb_vec_clear(r, len1);
+        _acb_vec_clear(r2, len2);
+        _acb_vec_clear(r3, len1 + 2 * len2);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}

doc/source/arb_poly.rst

@@ -560,6 +560,26 @@ Product trees
 
     Generates the polynomial `(x-x_0)(x-x_1)\cdots(x-x_{n-1})`.
 
+.. function:: void _arb_poly_product_roots_complex(arb_ptr poly, arb_srcptr r, slong rn, acb_srcptr c, slong cn, slong prec)
+
+.. function:: void arb_poly_product_roots_complex(arb_poly_t poly, arb_srcptr r, slong rn, acb_srcptr c, slong cn, slong prec)
+
+    Generates the polynomial
+
+    .. math ::
+
+        \left(\prod_{i=0}^{rn-1} (x-r_i)\right) \left(\prod_{i=0}^{cn-1} (x-c_i)(x-\bar{c_i})\right)
+
+    having *rn* real roots given by the array *r* and having `2cn` complex roots
+    in conjugate pairs given by the length-*cn* array *c*.
+    Either *rn* or *cn* or both may be zero.
+
+    Note that only one representative from each complex conjugate pair
+    is supplied (unless a pair is supposed to
+    be repeated with higher multiplicity).
+    To construct a polynomial from complex roots where the conjugate pairs
+    have not been distinguished, use :func:`acb_poly_product_roots` instead.
+
 .. function:: arb_ptr * _arb_poly_tree_alloc(slong len)
 
     Returns an initialized data structured capable of representing a