add complex multipoint evaluation/interpolation and helper functions

fmpcb_poly.h

@@ -158,6 +158,15 @@ void _fmpcb_poly_mul(fmpcb_struct * C,
 void fmpcb_poly_mul(fmpcb_poly_t res, const fmpcb_poly_t poly1,
               const fmpcb_poly_t poly2, long prec);
 
+static __inline__ void
+_fmpcb_poly_mul_monic(fmpcb_struct * res, const fmpcb_struct * poly1, long len1,
+    const fmpcb_struct * poly2, long len2, long prec)
+{
+    if (len1 + len2 - 2 > 0)
+        _fmpcb_poly_mullow(res, poly1, len1, poly2, len2, len1 + len2 - 2, prec);
+    fmpcb_one(res + len1 + len2 - 2);
+}
+
 void _fmpcb_poly_inv_series(fmpcb_struct * Qinv, const fmpcb_struct * Q, long len, long prec);
 
 void fmpcb_poly_inv_series(fmpcb_poly_t Qinv, const fmpcb_poly_t Q, long n, long prec);
@@ -179,6 +188,75 @@ void _fmpcb_poly_rem(fmpcb_struct * R,
 void fmpcb_poly_divrem(fmpcb_poly_t Q, fmpcb_poly_t R,
                              const fmpcb_poly_t A, const fmpcb_poly_t B, long prec);
 
+void _fmpcb_poly_div_root(fmpcb_struct * Q, fmpcb_t R, const fmpcb_struct * A,
+    long len, const fmpcb_t c, long prec);
+
+void
+_fmpcb_poly_evaluate_vec_fast_precomp(fmpcb_struct * vs, const fmpcb_struct * poly,
+    long plen, fmpcb_struct ** tree, long len, long prec);
+
+void _fmpcb_poly_evaluate_vec_fast(fmpcb_struct * ys, const fmpcb_struct * poly, long plen,
+    const fmpcb_struct * xs, long n, long prec);
+
+void
+fmpcb_poly_evaluate_vec_fast(fmpcb_struct * ys,
+        const fmpcb_poly_t poly, const fmpcb_struct * xs, long n, long prec);
+
+void
+_fmpcb_poly_evaluate_vec_iter(fmpcb_struct * ys, const fmpcb_struct * poly, long plen,
+    const fmpcb_struct * xs, long n, long prec);
+
+void
+fmpcb_poly_evaluate_vec_iter(fmpcb_struct * ys,
+        const fmpcb_poly_t poly, const fmpcb_struct * xs, long n, long prec);
+
+void
+_fmpcb_poly_interpolate_barycentric(fmpcb_struct * poly,
+    const fmpcb_struct * xs, const fmpcb_struct * ys, long n, long prec);
+
+void
+fmpcb_poly_interpolate_barycentric(fmpcb_poly_t poly,
+    const fmpcb_struct * xs, const fmpcb_struct * ys, long n, long prec);
+
+void
+_fmpcb_poly_interpolation_weights(fmpcb_struct * w,
+    fmpcb_struct ** tree, long len, long prec);
+
+void
+_fmpcb_poly_interpolate_fast_precomp(fmpcb_struct * poly,
+    const fmpcb_struct * ys, fmpcb_struct ** tree, const fmpcb_struct * weights,
+    long len, long prec);
+
+void
+_fmpcb_poly_interpolate_fast(fmpcb_struct * poly,
+    const fmpcb_struct * xs, const fmpcb_struct * ys, long len, long prec);
+
+void
+fmpcb_poly_interpolate_fast(fmpcb_poly_t poly,
+        const fmpcb_struct * xs, const fmpcb_struct * ys, long n, long prec);
+
+void
+_fmpcb_poly_interpolate_newton(fmpcb_struct * poly, const fmpcb_struct * xs,
+    const fmpcb_struct * ys, long n, long prec);
+
+void
+fmpcb_poly_interpolate_newton(fmpcb_poly_t poly,
+    const fmpcb_struct * xs, const fmpcb_struct * ys, long n, long prec);
+
+void
+_fmpcb_poly_product_roots(fmpcb_struct * poly, const fmpcb_struct * xs, long n, long prec);
+
+void
+fmpcb_poly_product_roots(fmpcb_poly_t poly, fmpcb_struct * xs, long n, long prec);
+
+fmpcb_struct ** _fmpcb_poly_tree_alloc(long len);
+
+void _fmpcb_poly_tree_free(fmpcb_struct ** tree, long len);
+
+void
+_fmpcb_poly_tree_build(fmpcb_struct ** tree, const fmpcb_struct * roots, long len, long prec);
+
+
 void _fmpcb_poly_root_inclusion(fmpcb_t r, const fmpcb_t m,
     const fmpcb_struct * poly,
     const fmpcb_struct * polyder, long len, long prec);

fmpcb_poly/div_root.c

@@ -0,0 +1,60 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2011 William Hart
+    Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmpcb_poly.h"
+
+void
+_fmpcb_poly_div_root(fmpcb_struct * Q, fmpcb_t R, const fmpcb_struct * A,
+    long len, const fmpcb_t c, long prec)
+{
+    fmpcb_t r, t;
+    long i;
+
+    if (len < 2)
+    {
+        fmpcb_zero(R);
+        return;
+    }
+
+    fmpcb_init(r);
+    fmpcb_init(t);
+
+    fmpcb_set(t, A + len - 2);
+    fmpcb_set(Q + len - 2, A + len - 1);
+    fmpcb_set(r, Q + len - 2);
+
+    /* TODO: avoid the extra assignments (but still support aliasing)  */
+    for (i = len - 2; i > 0; i--)
+    {
+        fmpcb_mul(r, r, c, prec);
+        fmpcb_add(r, r, t, prec);
+        fmpcb_set(t, A + i - 1);
+        fmpcb_set(Q + i - 1, r);
+    }
+
+    fmpcb_mul(r, r, c, prec);
+    fmpcb_add(R, r, t, prec);
+}

fmpcb_poly/evaluate_vec_fast.c

@@ -0,0 +1,149 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmpcb_poly.h"
+
+/* This gives some speedup for small lengths. */
+static __inline__ void
+_fmpcb_poly_rem_2(fmpcb_struct * r, const fmpcb_struct * a, long al,
+    const fmpcb_struct * b, long bl, long prec)
+{
+    if (al == 2)
+    {
+        fmpcb_mul(r + 0, a + 1, b + 0, prec);
+        fmpcb_sub(r + 0, a + 0, r + 0, prec);
+    }
+    else
+    {
+        _fmpcb_poly_rem(r, a, al, b, bl, prec);
+    }
+}
+
+void
+_fmpcb_poly_evaluate_vec_fast_precomp(fmpcb_struct * vs, const fmpcb_struct * poly,
+    long plen, fmpcb_struct ** tree, long len, long prec)
+{
+    long height, i, j, pow, left;
+    long tree_height;
+    long tlen;
+    fmpcb_struct *t, *u, *swap, *pa, *pb, *pc;
+
+    /* avoid worrying about some degenerate cases */
+    if (len < 2 || plen < 2)
+    {
+        if (len == 1)
+        {
+            fmpcb_t tmp;
+            fmpcb_init(tmp);
+            fmpcb_neg(tmp, tree[0] + 0);
+            _fmpcb_poly_evaluate(vs + 0, poly, plen, tmp, prec);
+            fmpcb_clear(tmp);
+        }
+        else if (len != 0 && plen == 0)
+        {
+            _fmpcb_vec_zero(vs, len);
+        }
+        else if (len != 0 && plen == 1)
+        {
+            for (i = 0; i < len; i++)
+                fmpcb_set(vs + i, poly + 0);
+        }
+        return;
+    }
+
+    t = _fmpcb_vec_init(len);
+    u = _fmpcb_vec_init(len);
+
+    left = len;
+
+    /* Initial reduction. We allow the polynomial to be larger
+        or smaller than the number of points. */
+    height = FLINT_BIT_COUNT(plen - 1) - 1;
+    tree_height = FLINT_CLOG2(len);
+    while (height >= tree_height)
+        height--;
+    pow = 1L << height;
+
+    for (i = j = 0; i < len; i += pow, j += (pow + 1))
+    {
+        tlen = ((i + pow) <= len) ? pow : len % pow;
+        _fmpcb_poly_rem(t + i, poly, plen, tree[height] + j, tlen + 1, prec);
+    }
+
+    for (i = height - 1; i >= 0; i--)
+    {
+        pow = 1L << i;
+        left = len;
+        pa = tree[i];
+        pb = t;
+        pc = u;
+
+        while (left >= 2 * pow)
+        {
+            _fmpcb_poly_rem_2(pc, pb, 2 * pow, pa, pow + 1, prec);
+            _fmpcb_poly_rem_2(pc + pow, pb, 2 * pow, pa + pow + 1, pow + 1, prec);
+
+            pa += 2 * pow + 2;
+            pb += 2 * pow;
+            pc += 2 * pow;
+            left -= 2 * pow;
+        }
+
+        if (left > pow)
+        {
+            _fmpcb_poly_rem(pc, pb, left, pa, pow + 1, prec);
+            _fmpcb_poly_rem(pc + pow, pb, left, pa + pow + 1, left - pow + 1, prec);
+        }
+        else if (left > 0)
+            _fmpcb_vec_set(pc, pb, left);
+
+        swap = t;
+        t = u;
+        u = swap;
+    }
+
+    _fmpcb_vec_set(vs, t, len);
+    _fmpcb_vec_clear(t, len);
+    _fmpcb_vec_clear(u, len);
+}
+
+void _fmpcb_poly_evaluate_vec_fast(fmpcb_struct * ys, const fmpcb_struct * poly, long plen,
+    const fmpcb_struct * xs, long n, long prec)
+{
+    fmpcb_struct ** tree;
+
+    tree = _fmpcb_poly_tree_alloc(n);
+    _fmpcb_poly_tree_build(tree, xs, n, prec);
+    _fmpcb_poly_evaluate_vec_fast_precomp(ys, poly, plen, tree, n, prec);
+    _fmpcb_poly_tree_free(tree, n);
+}
+
+void
+fmpcb_poly_evaluate_vec_fast(fmpcb_struct * ys,
+        const fmpcb_poly_t poly, const fmpcb_struct * xs, long n, long prec)
+{
+    _fmpcb_poly_evaluate_vec_fast(ys, poly->coeffs,
+                                        poly->length, xs, n, prec);
+}

fmpcb_poly/evaluate_vec_iter.c

@@ -0,0 +1,44 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmpcb_poly.h"
+
+void
+_fmpcb_poly_evaluate_vec_iter(fmpcb_struct * ys, const fmpcb_struct * poly, long plen,
+    const fmpcb_struct * xs, long n, long prec)
+{
+    long i;
+
+    for (i = 0; i < n; i++)
+        _fmpcb_poly_evaluate(ys + i, poly, plen, xs + i, prec);
+}
+
+void
+fmpcb_poly_evaluate_vec_iter(fmpcb_struct * ys,
+        const fmpcb_poly_t poly, const fmpcb_struct * xs, long n, long prec)
+{
+    _fmpcb_poly_evaluate_vec_iter(ys, poly->coeffs,
+                                        poly->length, xs, n, prec);
+}

fmpcb_poly/interpolate_barycentric.c

@@ -0,0 +1,96 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmpcb_poly.h"
+
+void
+_fmpcb_poly_interpolate_barycentric(fmpcb_struct * poly,
+    const fmpcb_struct * xs, const fmpcb_struct * ys, long n, long prec)
+{
+    fmpcb_struct *P, *Q, *w;
+    fmpcb_t t;
+    long i, j;
+
+    if (n == 1)
+    {
+        fmpcb_set(poly, ys);
+        return;
+    }
+
+    P = _fmpcb_vec_init(n + 1);
+    Q = _fmpcb_vec_init(n);
+    w = _fmpcb_vec_init(n);
+    fmpcb_init(t);
+
+    _fmpcb_poly_product_roots(P, xs, n, prec);
+
+    for (i = 0; i < n; i++)
+    {
+        fmpcb_one(w + i);
+
+        for (j = 0; j < n; j++)
+        {
+            if (i != j)
+            {
+                fmpcb_sub(t, xs + i, xs + j, prec);
+                fmpcb_mul(w + i, w + i, t, prec);
+            }
+        }
+
+        fmpcb_inv(w + i, w + i, prec);
+    }
+
+    _fmpcb_vec_zero(poly, n);
+
+    for (i = 0; i < n; i++)
+    {
+        _fmpcb_poly_div_root(Q, t, P, n + 1, xs + i, prec);
+        fmpcb_mul(t, w + i, ys + i, prec);
+        _fmpcb_vec_scalar_addmul(poly, Q, n, t, prec);
+    }
+
+    _fmpcb_vec_clear(P, n + 1);
+    _fmpcb_vec_clear(Q, n);
+    _fmpcb_vec_clear(w, n);
+    fmpcb_clear(t);
+}
+
+void
+fmpcb_poly_interpolate_barycentric(fmpcb_poly_t poly,
+    const fmpcb_struct * xs, const fmpcb_struct * ys, long n, long prec)
+{
+    if (n == 0)
+    {
+        fmpcb_poly_zero(poly);
+    }
+    else
+    {
+        fmpcb_poly_fit_length(poly, n);
+        _fmpcb_poly_set_length(poly, n);
+        _fmpcb_poly_interpolate_barycentric(poly->coeffs,
+            xs, ys, n, prec);
+        _fmpcb_poly_normalise(poly);
+    }
+}

fmpcb_poly/interpolate_fast.c

@@ -0,0 +1,141 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmpcb_poly.h"
+
+void
+_fmpcb_poly_interpolation_weights(fmpcb_struct * w,
+    fmpcb_struct ** tree, long len, long prec)
+{
+    fmpcb_struct * tmp;
+    long i, n, height;
+
+    if (len == 0)
+        return;
+
+    if (len == 1)
+    {
+        fmpcb_one(w);
+        return;
+    }
+
+    tmp = _fmpcb_vec_init(len + 1);
+    height = FLINT_CLOG2(len);
+    n = 1L << (height - 1);
+
+    _fmpcb_poly_mul_monic(tmp, tree[height-1], n + 1,
+                        tree[height-1] + (n + 1), (len - n + 1), prec);
+
+    _fmpcb_poly_derivative(tmp, tmp, len + 1, prec);
+    _fmpcb_poly_evaluate_vec_fast_precomp(w, tmp, len, tree, len, prec);
+
+    for (i = 0; i < len; i++)
+        fmpcb_inv(w + i, w + i, prec);
+
+    _fmpcb_vec_clear(tmp, len + 1);
+}
+
+void
+_fmpcb_poly_interpolate_fast_precomp(fmpcb_struct * poly,
+    const fmpcb_struct * ys, fmpcb_struct ** tree, const fmpcb_struct * weights,
+    long len, long prec)
+{
+    fmpcb_struct *t, *u, *pa, *pb;
+    long i, pow, left;
+
+    if (len == 0)
+        return;
+
+    t = _fmpcb_vec_init(len);
+    u = _fmpcb_vec_init(len);
+
+    for (i = 0; i < len; i++)
+        fmpcb_mul(poly + i, weights + i, ys + i, prec);
+
+    for (i = 0; i < FLINT_CLOG2(len); i++)
+    {
+        pow = (1L << i);
+        pa = tree[i];
+        pb = poly;
+        left = len;
+
+        while (left >= 2 * pow)
+        {
+            _fmpcb_poly_mul(t, pa, pow + 1, pb + pow, pow, prec);
+            _fmpcb_poly_mul(u, pa + pow + 1, pow + 1, pb, pow, prec);
+            _fmpcb_vec_add(pb, t, u, 2 * pow, prec);
+
+            left -= 2 * pow;
+            pa += 2 * pow + 2;
+            pb += 2 * pow;
+        }
+
+        if (left > pow)
+        {
+            _fmpcb_poly_mul(t, pa, pow + 1, pb + pow, left - pow, prec);
+            _fmpcb_poly_mul(u, pb, pow, pa + pow + 1, left - pow + 1, prec);
+            _fmpcb_vec_add(pb, t, u, left, prec);
+        }
+    }
+
+    _fmpcb_vec_clear(t, len);
+    _fmpcb_vec_clear(u, len);
+}
+
+void
+_fmpcb_poly_interpolate_fast(fmpcb_struct * poly,
+    const fmpcb_struct * xs, const fmpcb_struct * ys, long len, long prec)
+{
+    fmpcb_struct ** tree;
+    fmpcb_struct *w;
+
+    tree = _fmpcb_poly_tree_alloc(len);
+    _fmpcb_poly_tree_build(tree, xs, len, prec);
+
+    w = _fmpcb_vec_init(len);
+    _fmpcb_poly_interpolation_weights(w, tree, len, prec);
+
+    _fmpcb_poly_interpolate_fast_precomp(poly, ys, tree, w, len, prec);
+
+    _fmpcb_vec_clear(w, len);
+    _fmpcb_poly_tree_free(tree, len);
+}
+
+void
+fmpcb_poly_interpolate_fast(fmpcb_poly_t poly,
+        const fmpcb_struct * xs, const fmpcb_struct * ys, long n, long prec)
+{
+    if (n == 0)
+    {
+        fmpcb_poly_zero(poly);
+    }
+    else
+    {
+        fmpcb_poly_fit_length(poly, n);
+        _fmpcb_poly_set_length(poly, n);
+        _fmpcb_poly_interpolate_fast(poly->coeffs, xs, ys, n, prec);
+        _fmpcb_poly_normalise(poly);
+    }
+}

fmpcb_poly/interpolate_newton.c

@@ -0,0 +1,119 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmpcb_poly.h"
+
+static void
+_interpolate_newton(fmpcb_struct * ys, const fmpcb_struct * xs, long n, long prec)
+{
+    fmpcb_t p, q, t;
+    long i, j;
+
+    fmpcb_init(p);
+    fmpcb_init(q);
+    fmpcb_init(t);
+
+    for (i = 1; i < n; i++)
+    {
+        fmpcb_set(t, ys + i - 1);
+
+        for (j = i; j < n; j++)
+        {
+            fmpcb_sub(p, ys + j, t, prec);
+            fmpcb_sub(q, xs + j, xs + j - i, prec);
+            fmpcb_set(t, ys + j);
+            fmpcb_div(ys + j, p, q, prec);
+        }
+    }
+
+    fmpcb_clear(p);
+    fmpcb_clear(q);
+    fmpcb_clear(t);
+}
+
+static void
+_newton_to_monomial(fmpcb_struct * ys, const fmpcb_struct * xs, long n, long prec)
+{
+    fmpcb_t t, u;
+    long i, j;
+
+    fmpcb_init(t);
+    fmpcb_init(u);
+
+    for (i = n - 2; i >= 0; i--)
+    {
+        fmpcb_set(t, ys + i);
+        fmpcb_set(ys + i, ys + i + 1);
+
+        for (j = i + 1; j < n - 1; j++)
+        {
+            fmpcb_mul(u, ys + j, xs + i, prec);
+            fmpcb_sub(ys + j, ys + j + 1, u, prec);
+        }
+
+        fmpcb_mul(u, ys + n - 1, xs + i, prec);
+        fmpcb_sub(ys + n - 1, t, u, prec);
+    }
+
+    _fmpcb_poly_reverse(ys, ys, n, n);
+
+    fmpcb_clear(t);
+    fmpcb_clear(u);
+}
+
+void
+_fmpcb_poly_interpolate_newton(fmpcb_struct * poly, const fmpcb_struct * xs,
+    const fmpcb_struct * ys, long n, long prec)
+{
+    if (n == 1)
+    {
+        fmpcb_set(poly, ys);
+    }
+    else
+    {
+        _fmpcb_vec_set(poly, ys, n);
+        _interpolate_newton(poly, xs, n, prec);
+        while (n > 0 && fmpcb_is_zero(poly + n - 1)) n--;
+        _newton_to_monomial(poly, xs, n, prec);
+    }
+}
+
+void
+fmpcb_poly_interpolate_newton(fmpcb_poly_t poly,
+    const fmpcb_struct * xs, const fmpcb_struct * ys, long n, long prec)
+{
+    if (n == 0)
+    {
+        fmpcb_poly_zero(poly);
+    }
+    else
+    {
+        fmpcb_poly_fit_length(poly, n);
+        _fmpcb_poly_set_length(poly, n);
+        _fmpcb_poly_interpolate_newton(poly->coeffs,
+            xs, ys, n, prec);
+        _fmpcb_poly_normalise(poly);
+    }
+}

fmpcb_poly/product_roots.c

@@ -0,0 +1,69 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2011 William Hart
+    Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmpcb_poly.h"
+
+void
+_fmpcb_poly_product_roots(fmpcb_struct * poly, const fmpcb_struct * xs, long n, long prec)
+{
+    if (n == 0)
+    {
+        fmpcb_one(poly);
+    }
+    else if (n == 1)
+    {
+        fmpcb_neg(poly, xs);
+        fmpcb_one(poly + 1);
+    }
+    else if (n == 2)
+    {
+        fmpcb_mul(poly, xs + 0, xs + 1, prec);
+        fmpcb_add(poly + 1, xs + 0, xs + 1, prec);
+        fmpcb_neg(poly + 1, poly + 1);
+        fmpcb_one(poly + 2);
+    }
+    else
+    {
+        const long m = (n + 1) / 2;
+        fmpcb_struct * tmp;
+
+        tmp = _fmpcb_vec_init(n + 2);
+
+        _fmpcb_poly_product_roots(tmp, xs, m, prec);
+        _fmpcb_poly_product_roots(tmp + m + 1, xs + m, n - m, prec);
+        _fmpcb_poly_mul_monic(poly, tmp, m + 1, tmp + m + 1, n - m + 1, prec);
+
+        _fmpcb_vec_clear(tmp, n + 2);
+    }
+}
+
+void
+fmpcb_poly_product_roots(fmpcb_poly_t poly, fmpcb_struct * xs, long n, long prec)
+{
+    fmpcb_poly_fit_length(poly, n + 1);
+    _fmpcb_poly_product_roots(poly->coeffs, xs, n, prec);
+    _fmpcb_poly_set_length(poly, n + 1);
+}

fmpcb_poly/test/t-evaluate_vec_fast.c

@@ -0,0 +1,104 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmpcb_poly.h"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("evaluate_vec_fast....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 10000; iter++)
+    {
+        long i, n, qbits1, qbits2, rbits1, rbits2, rbits3;
+        fmpq_poly_t F;
+        fmpq * X, * Y;
+        fmpcb_poly_t f;
+        fmpcb_struct * x, * y;
+
+        qbits1 = 2 + n_randint(state, 100);
+        qbits2 = 2 + n_randint(state, 100);
+        rbits1 = 2 + n_randint(state, 200);
+        rbits2 = 2 + n_randint(state, 200);
+        rbits3 = 2 + n_randint(state, 200);
+
+        n = n_randint(state, 10);
+
+        fmpq_poly_init(F);
+        X = _fmpq_vec_init(n);
+        Y = _fmpq_vec_init(n);
+
+        fmpcb_poly_init(f);
+        x = _fmpcb_vec_init(n);
+        y = _fmpcb_vec_init(n);
+
+        fmpq_poly_randtest(F, state, 1 + n_randint(state, 20), qbits1);
+        for (i = 0; i < n; i++)
+            fmpq_randtest(X + i, state, qbits2);
+        for (i = 0; i < n; i++)
+            fmpq_poly_evaluate_fmpq(Y + i, F, X + i);
+
+        fmpcb_poly_set_fmpq_poly(f, F, rbits1);
+        for (i = 0; i < n; i++)
+            fmpcb_set_fmpq(x + i, X + i, rbits2);
+        fmpcb_poly_evaluate_vec_fast(y, f, x, n, rbits3);
+
+        for (i = 0; i < n; i++)
+        {
+            if (!fmpcb_contains_fmpq(y + i, Y + i))
+            {
+                printf("FAIL (%ld of %ld)\n\n", i, n);
+
+                printf("F = "); fmpq_poly_print(F); printf("\n\n");
+                printf("X = "); fmpq_print(X + i); printf("\n\n");
+                printf("Y = "); fmpq_print(Y + i); printf("\n\n");
+
+                printf("f = "); fmpcb_poly_printd(f, 15); printf("\n\n");
+                printf("x = "); fmpcb_printd(x + i, 15); printf("\n\n");
+                printf("y = "); fmpcb_printd(y + i, 15); printf("\n\n");
+
+                abort();
+            }
+        }
+
+        fmpq_poly_clear(F);
+        _fmpq_vec_clear(X, n);
+        _fmpq_vec_clear(Y, n);
+
+        fmpcb_poly_clear(f);
+        _fmpcb_vec_clear(x, n);
+        _fmpcb_vec_clear(y, n);
+    }
+
+    flint_randclear(state);
+    _fmpz_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}

fmpcb_poly/test/t-evaluate_vec_iter.c

@@ -0,0 +1,104 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmpcb_poly.h"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("evaluate_vec_iter....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 10000; iter++)
+    {
+        long i, n, qbits1, qbits2, rbits1, rbits2, rbits3;
+        fmpq_poly_t F;
+        fmpq * X, * Y;
+        fmpcb_poly_t f;
+        fmpcb_struct * x, * y;
+
+        qbits1 = 2 + n_randint(state, 100);
+        qbits2 = 2 + n_randint(state, 100);
+        rbits1 = 2 + n_randint(state, 200);
+        rbits2 = 2 + n_randint(state, 200);
+        rbits3 = 2 + n_randint(state, 200);
+
+        n = n_randint(state, 10);
+
+        fmpq_poly_init(F);
+        X = _fmpq_vec_init(n);
+        Y = _fmpq_vec_init(n);
+
+        fmpcb_poly_init(f);
+        x = _fmpcb_vec_init(n);
+        y = _fmpcb_vec_init(n);
+
+        fmpq_poly_randtest(F, state, 1 + n_randint(state, 20), qbits1);
+        for (i = 0; i < n; i++)
+            fmpq_randtest(X + i, state, qbits2);
+        for (i = 0; i < n; i++)
+            fmpq_poly_evaluate_fmpq(Y + i, F, X + i);
+
+        fmpcb_poly_set_fmpq_poly(f, F, rbits1);
+        for (i = 0; i < n; i++)
+            fmpcb_set_fmpq(x + i, X + i, rbits2);
+        fmpcb_poly_evaluate_vec_iter(y, f, x, n, rbits3);
+
+        for (i = 0; i < n; i++)
+        {
+            if (!fmpcb_contains_fmpq(y + i, Y + i))
+            {
+                printf("FAIL (%ld of %ld)\n\n", i, n);
+
+                printf("F = "); fmpq_poly_print(F); printf("\n\n");
+                printf("X = "); fmpq_print(X + i); printf("\n\n");
+                printf("Y = "); fmpq_print(Y + i); printf("\n\n");
+
+                printf("f = "); fmpcb_poly_printd(f, 15); printf("\n\n");
+                printf("x = "); fmpcb_printd(x + i, 15); printf("\n\n");
+                printf("y = "); fmpcb_printd(y + i, 15); printf("\n\n");
+
+                abort();
+            }
+        }
+
+        fmpq_poly_clear(F);
+        _fmpq_vec_clear(X, n);
+        _fmpq_vec_clear(Y, n);
+
+        fmpcb_poly_clear(f);
+        _fmpcb_vec_clear(x, n);
+        _fmpcb_vec_clear(y, n);
+    }
+
+    flint_randclear(state);
+    _fmpz_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}

fmpcb_poly/test/t-interpolate_barycentric.c

@@ -0,0 +1,107 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmpcb_poly.h"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("interpolate_barycentric....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 10000; iter++)
+    {
+        long i, n, qbits1, qbits2, rbits1, rbits2, rbits3;
+        fmpq_poly_t P;
+        fmpcb_poly_t R, S;
+        fmpq_t t, u;
+        fmpcb_struct * xs, * ys;
+
+        fmpq_poly_init(P);
+        fmpcb_poly_init(R);
+        fmpcb_poly_init(S);
+        fmpq_init(t);
+        fmpq_init(u);
+
+        qbits1 = 2 + n_randint(state, 200);
+        qbits2 = 2 + n_randint(state, 5);
+        rbits1 = 2 + n_randint(state, 200);
+        rbits2 = 2 + n_randint(state, 200);
+        rbits3 = 2 + n_randint(state, 200);
+
+        fmpq_poly_randtest(P, state, 1 + n_randint(state, 20), qbits1);
+        n = P->length;
+
+        xs = _fmpcb_vec_init(n);
+        ys = _fmpcb_vec_init(n);
+
+        fmpcb_poly_set_fmpq_poly(R, P, rbits1);
+
+        if (n > 0)
+        {
+            fmpq_randtest(t, state, qbits2);
+            fmpcb_set_fmpq(xs, t, rbits2);
+
+            for (i = 1; i < n; i++)
+            {
+                fmpq_randtest_not_zero(u, state, qbits2);
+                fmpq_abs(u, u);
+                fmpq_add(t, t, u);
+                fmpcb_set_fmpq(xs + i, t, rbits2);
+            }
+        }
+
+        for (i = 0; i < n; i++)
+            fmpcb_poly_evaluate(ys + i, R, xs + i, rbits2);
+
+        fmpcb_poly_interpolate_barycentric(S, xs, ys, n, rbits3);
+
+        if (!fmpcb_poly_contains_fmpq_poly(S, P))
+        {
+            printf("FAIL:\n");
+            printf("P = "); fmpq_poly_print(P); printf("\n\n");
+            printf("R = "); fmpcb_poly_printd(R, 15); printf("\n\n");
+            printf("S = "); fmpcb_poly_printd(S, 15); printf("\n\n");
+            abort();
+        }
+
+        fmpq_poly_clear(P);
+        fmpcb_poly_clear(R);
+        fmpcb_poly_clear(S);
+        fmpq_clear(t);
+        fmpq_clear(u);
+        _fmpcb_vec_clear(xs, n);
+        _fmpcb_vec_clear(ys, n);
+    }
+
+    flint_randclear(state);
+    _fmpz_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}

fmpcb_poly/test/t-interpolate_fast.c

@@ -0,0 +1,107 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmpcb_poly.h"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("interpolate_fast....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 10000; iter++)
+    {
+        long i, n, qbits1, qbits2, rbits1, rbits2, rbits3;
+        fmpq_poly_t P;
+        fmpcb_poly_t R, S;
+        fmpq_t t, u;
+        fmpcb_struct * xs, * ys;
+
+        fmpq_poly_init(P);
+        fmpcb_poly_init(R);
+        fmpcb_poly_init(S);
+        fmpq_init(t);
+        fmpq_init(u);
+
+        qbits1 = 2 + n_randint(state, 200);
+        qbits2 = 2 + n_randint(state, 5);
+        rbits1 = 2 + n_randint(state, 200);
+        rbits2 = 2 + n_randint(state, 200);
+        rbits3 = 2 + n_randint(state, 200);
+
+        fmpq_poly_randtest(P, state, 1 + n_randint(state, 20), qbits1);
+        n = P->length;
+
+        xs = _fmpcb_vec_init(n);
+        ys = _fmpcb_vec_init(n);
+
+        fmpcb_poly_set_fmpq_poly(R, P, rbits1);
+
+        if (n > 0)
+        {
+            fmpq_randtest(t, state, qbits2);
+            fmpcb_set_fmpq(xs, t, rbits2);
+
+            for (i = 1; i < n; i++)
+            {
+                fmpq_randtest_not_zero(u, state, qbits2);
+                fmpq_abs(u, u);
+                fmpq_add(t, t, u);
+                fmpcb_set_fmpq(xs + i, t, rbits2);
+            }
+        }
+
+        for (i = 0; i < n; i++)
+            fmpcb_poly_evaluate(ys + i, R, xs + i, rbits2);
+
+        fmpcb_poly_interpolate_fast(S, xs, ys, n, rbits3);
+
+        if (!fmpcb_poly_contains_fmpq_poly(S, P))
+        {
+            printf("FAIL:\n");
+            printf("P = "); fmpq_poly_print(P); printf("\n\n");
+            printf("R = "); fmpcb_poly_printd(R, 15); printf("\n\n");
+            printf("S = "); fmpcb_poly_printd(S, 15); printf("\n\n");
+            abort();
+        }
+
+        fmpq_poly_clear(P);
+        fmpcb_poly_clear(R);
+        fmpcb_poly_clear(S);
+        fmpq_clear(t);
+        fmpq_clear(u);
+        _fmpcb_vec_clear(xs, n);
+        _fmpcb_vec_clear(ys, n);
+    }
+
+    flint_randclear(state);
+    _fmpz_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}

fmpcb_poly/test/t-interpolate_newton.c

@@ -0,0 +1,107 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmpcb_poly.h"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("interpolate_newton....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 10000; iter++)
+    {
+        long i, n, qbits1, qbits2, rbits1, rbits2, rbits3;
+        fmpq_poly_t P;
+        fmpcb_poly_t R, S;
+        fmpq_t t, u;
+        fmpcb_struct * xs, * ys;
+
+        fmpq_poly_init(P);
+        fmpcb_poly_init(R);
+        fmpcb_poly_init(S);
+        fmpq_init(t);
+        fmpq_init(u);
+
+        qbits1 = 2 + n_randint(state, 200);
+        qbits2 = 2 + n_randint(state, 5);
+        rbits1 = 2 + n_randint(state, 200);
+        rbits2 = 2 + n_randint(state, 200);
+        rbits3 = 2 + n_randint(state, 200);
+
+        fmpq_poly_randtest(P, state, 1 + n_randint(state, 20), qbits1);
+        n = P->length;
+
+        xs = _fmpcb_vec_init(n);
+        ys = _fmpcb_vec_init(n);
+
+        fmpcb_poly_set_fmpq_poly(R, P, rbits1);
+
+        if (n > 0)
+        {
+            fmpq_randtest(t, state, qbits2);
+            fmpcb_set_fmpq(xs, t, rbits2);
+
+            for (i = 1; i < n; i++)
+            {
+                fmpq_randtest_not_zero(u, state, qbits2);
+                fmpq_abs(u, u);
+                fmpq_add(t, t, u);
+                fmpcb_set_fmpq(xs + i, t, rbits2);
+            }
+        }
+
+        for (i = 0; i < n; i++)
+            fmpcb_poly_evaluate(ys + i, R, xs + i, rbits2);
+
+        fmpcb_poly_interpolate_newton(S, xs, ys, n, rbits3);
+
+        if (!fmpcb_poly_contains_fmpq_poly(S, P))
+        {
+            printf("FAIL:\n");
+            printf("P = "); fmpq_poly_print(P); printf("\n\n");
+            printf("R = "); fmpcb_poly_printd(R, 15); printf("\n\n");
+            printf("S = "); fmpcb_poly_printd(S, 15); printf("\n\n");
+            abort();
+        }
+
+        fmpq_poly_clear(P);
+        fmpcb_poly_clear(R);
+        fmpcb_poly_clear(S);
+        fmpq_clear(t);
+        fmpq_clear(u);
+        _fmpcb_vec_clear(xs, n);
+        _fmpcb_vec_clear(ys, n);
+    }
+
+    flint_randclear(state);
+    _fmpz_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}

fmpcb_poly/tree.c

@@ -0,0 +1,120 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmpcb_poly.h"
+
+fmpcb_struct ** _fmpcb_poly_tree_alloc(long len)
+{
+    fmpcb_struct ** tree = NULL;
+
+    if (len)
+    {
+        long i, height = FLINT_CLOG2(len);
+
+        tree = flint_malloc(sizeof(fmpcb_struct *) * (height + 1));
+        for (i = 0; i <= height; i++)
+            tree[i] = _fmpcb_vec_init(len + (len >> i) + 1);
+    }
+
+    return tree;
+}
+
+void _fmpcb_poly_tree_free(fmpcb_struct ** tree, long len)
+{
+    if (len)
+    {
+        long i, height = FLINT_CLOG2(len);
+
+        for (i = 0; i <= height; i++)
+            _fmpcb_vec_clear(tree[i], len + (len >> i) + 1);
+
+        flint_free(tree);
+    }
+}
+
+void
+_fmpcb_poly_tree_build(fmpcb_struct ** tree, const fmpcb_struct * roots, long len, long prec)
+{
+    long height, pow, left, i;
+    fmpcb_struct *pa, *pb, *a, *b;
+
+    if (len == 0)
+        return;
+
+    height = FLINT_CLOG2(len);
+
+    /* zeroth level, (x-a) */
+    for (i = 0; i < len; i++)
+    {
+        fmpcb_one(tree[0] + (2 * i + 1));
+        fmpcb_neg(tree[0] + (2 * i), roots + i);
+    }
+
+    /* first level, (x-a)(x-b) = x^2 + (-a-b)*x + a*b */
+    if (height > 1)
+    {
+        pa = tree[1];
+
+        for (i = 0; i < len / 2; i++)
+        {
+            a = (fmpcb_struct *) (roots + (2 * i));
+            b = (fmpcb_struct *) (roots + (2 * i + 1));
+
+            fmpcb_mul(pa + (3 * i), a, b, prec);
+            fmpcb_add(pa + (3 * i + 1), a, b, prec);
+            fmpcb_neg(pa + (3 * i + 1), pa + (3 * i + 1));
+            fmpcb_one(pa + (3 * i + 2));
+        }
+
+        if (len & 1)
+        {
+            fmpcb_neg(pa + (3 * (len / 2)), roots + len - 1);
+            fmpcb_one(pa + (3 * (len / 2) + 1));
+        }
+    }
+
+    for (i = 1; i < height - 1; i++)
+    {
+        left = len;
+        pow = 1L << i;
+        pa = tree[i];
+        pb = tree[i + 1];
+
+        while (left >= 2 * pow)
+        {
+            _fmpcb_poly_mul_monic(pb, pa, pow + 1, pa + pow + 1, pow + 1, prec);
+            left -= 2 * pow;
+            pa += 2 * pow + 2;
+            pb += 2 * pow + 1;
+        }
+
+        if (left > pow)
+        {
+            _fmpcb_poly_mul_monic(pb, pa, pow + 1, pa + pow + 1, left - pow + 1, prec);
+        }
+        else if (left > 0)
+            _fmpcb_vec_set(pb, pa, left + 1);
+    }
+}