implement power series composition

doc/fmprb_poly.txt

@@ -249,6 +249,38 @@ void fmprb_poly_compose(fmprb_poly_t res,
     The underscore methods do not support aliasing of the output
     with either input polynomial.
 
+void _fmprb_poly_compose_series_horner(fmprb_struct * res,
+        const fmprb_struct * poly1, long len1,
+        const fmprb_struct * poly2, long len2, long n, long prec)
+
+void fmprb_poly_compose_series_horner(fmprb_poly_t res,
+                    const fmprb_poly_t poly1,
+                    const fmprb_poly_t poly2, long n, long prec)
+
+void _fmprb_poly_compose_series_brent_kung(fmprb_struct * res,
+            const fmprb_struct * poly1, long len1,
+            const fmprb_struct * poly2, long len2, long n, long prec)
+
+void fmprb_poly_compose_series_brent_kung(fmprb_poly_t res,
+                    const fmprb_poly_t poly1,
+                    const fmprb_poly_t poly2, long n, long prec)
+
+void _fmprb_poly_compose_series(fmprb_struct * res,
+            const fmprb_struct * poly1, long len1,
+            const fmprb_struct * poly2, long len2, long n, long prec)
+
+void fmprb_poly_compose_series(fmprb_poly_t res,
+                    const fmprb_poly_t poly1,
+                    const fmprb_poly_t poly2, long n, long prec)
+
+    Sets res to the power series composition $h(x) = f(g(x))$ truncated
+    to order $O(x^n)$ where $f$ is given by poly1 and $g$ is given by poly2,
+    respectively using Horner's rule, the Brent-Kung baby step-giant step
+    algorithm, and an automatic choice between the two algorithms.
+    We require that the constant term in $g(x)$ is exactly zero.
+    The underscore methods do not support aliasing of the output
+    with either input polynomial.
+
 *******************************************************************************
 
     Evaluation and interpolation

fmprb_poly.h

@@ -217,6 +217,27 @@ void _fmprb_poly_compose_divconquer(fmprb_struct * res,
 void fmprb_poly_compose_divconquer(fmprb_poly_t res,
               const fmprb_poly_t poly1, const fmprb_poly_t poly2, long prec);
 
+void _fmprb_poly_compose_series_horner(fmprb_struct * res, const fmprb_struct * poly1, long len1,
+                            const fmprb_struct * poly2, long len2, long n, long prec);
+
+void fmprb_poly_compose_series_horner(fmprb_poly_t res,
+                    const fmprb_poly_t poly1,
+                    const fmprb_poly_t poly2, long n, long prec);
+
+void _fmprb_poly_compose_series_brent_kung(fmprb_struct * res, const fmprb_struct * poly1, long len1,
+                            const fmprb_struct * poly2, long len2, long n, long prec);
+
+void fmprb_poly_compose_series_brent_kung(fmprb_poly_t res,
+                    const fmprb_poly_t poly1,
+                    const fmprb_poly_t poly2, long n, long prec);
+
+void _fmprb_poly_compose_series(fmprb_struct * res, const fmprb_struct * poly1, long len1,
+                            const fmprb_struct * poly2, long len2, long n, long prec);
+
+void fmprb_poly_compose_series(fmprb_poly_t res,
+                    const fmprb_poly_t poly1,
+                    const fmprb_poly_t poly2, long n, long prec);
+
 /* Evaluation and interpolation */
 
 void

fmprb_poly/compose_series.c

@@ -0,0 +1,89 @@
+/*=============================================================================
+
+    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 "fmprb_poly.h"
+
+void
+_fmprb_poly_compose_series(fmprb_struct * res, const fmprb_struct * poly1, long len1,
+                            const fmprb_struct * poly2, long len2, long n, long prec)
+{
+    if (len1 < 6 || n < 6)
+        _fmprb_poly_compose_series_horner(res, poly1, len1, poly2, len2, n, prec);
+    else
+        _fmprb_poly_compose_series_brent_kung(res, poly1, len1, poly2, len2, n, prec);
+}
+
+void
+fmprb_poly_compose_series(fmprb_poly_t res,
+                    const fmprb_poly_t poly1,
+                    const fmprb_poly_t poly2, long n, long prec)
+{
+    long len1 = poly1->length;
+    long len2 = poly2->length;
+    long lenr;
+
+    if (len2 != 0 && !fmprb_is_zero(poly2->coeffs))
+    {
+        printf("exception: compose_series: inner "
+                "polynomial must have zero constant term\n");
+        abort();
+    }
+
+    if (len1 == 0 || n == 0)
+    {
+        fmprb_poly_zero(res);
+        return;
+    }
+
+    if (len2 == 0 || len1 == 1)
+    {
+        fmprb_poly_set_fmprb(res, poly1->coeffs);
+        return;
+    }
+
+    lenr = FLINT_MIN((len1 - 1) * (len2 - 1) + 1, n);
+    len1 = FLINT_MIN(len1, lenr);
+    len2 = FLINT_MIN(len2, lenr);
+
+    if ((res != poly1) && (res != poly2))
+    {
+        fmprb_poly_fit_length(res, lenr);
+        _fmprb_poly_compose_series(res->coeffs, poly1->coeffs, len1,
+                                        poly2->coeffs, len2, lenr, prec);
+        _fmprb_poly_set_length(res, lenr);
+        _fmprb_poly_normalise(res);
+    }
+    else
+    {
+        fmprb_poly_t t;
+        fmprb_poly_init2(t, lenr);
+        _fmprb_poly_compose_series(t->coeffs, poly1->coeffs, len1,
+                                        poly2->coeffs, len2, lenr, prec);
+        _fmprb_poly_set_length(t, lenr);
+        _fmprb_poly_normalise(t);
+        fmprb_poly_swap(res, t);
+        fmprb_poly_clear(t);
+    }
+}

fmprb_poly/compose_series_brent_kung.c

@@ -0,0 +1,135 @@
+/*=============================================================================
+
+    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 "fmprb_poly.h"
+#include "fmprb_mat.h"
+
+void
+_fmprb_poly_compose_series_brent_kung(fmprb_struct * res,
+    const fmprb_struct * poly1, long len1,
+    const fmprb_struct * poly2, long len2, long n, long prec)
+{
+    fmprb_mat_t A, B, C;
+    fmprb_struct *t, *h;
+    long i, m;
+
+    if (n == 1)
+    {
+        fmprb_set(res, poly1);
+        return;
+    }
+
+    m = n_sqrt(n) + 1;
+
+    fmprb_mat_init(A, m, n);
+    fmprb_mat_init(B, m, m);
+    fmprb_mat_init(C, m, n);
+
+    h = _fmprb_vec_init(n);
+    t = _fmprb_vec_init(n);
+
+    /* Set rows of B to the segments of poly1 */
+    for (i = 0; i < len1 / m; i++)
+        _fmprb_vec_set(B->rows[i], poly1 + i*m, m);
+    _fmprb_vec_set(B->rows[i], poly1 + i*m, len1 % m);
+
+    /* Set rows of A to powers of poly2 */
+    fmprb_set_ui(A->rows[0] + 0, 1UL);
+    _fmprb_vec_set(A->rows[1], poly2, len2);
+    for (i = 2; i < m; i++)
+        _fmprb_poly_mullow(A->rows[i], A->rows[i-1], n, poly2, len2, n, prec);
+
+    fmprb_mat_mul(C, B, A, prec);
+
+    /* Evaluate block composition using the Horner scheme */
+    _fmprb_vec_set(res, C->rows[m - 1], n);
+    _fmprb_poly_mullow(h, A->rows[m - 1], n, poly2, len2, n, prec);
+
+    for (i = m - 2; i >= 0; i--)
+    {
+        _fmprb_poly_mullow(t, res, n, h, n, n, prec);
+        _fmprb_poly_add(res, t, n, C->rows[i], n, prec);
+    }
+
+    _fmprb_vec_clear(h, n);
+    _fmprb_vec_clear(t, n);
+
+    fmprb_mat_clear(A);
+    fmprb_mat_clear(B);
+    fmprb_mat_clear(C);
+}
+
+void
+fmprb_poly_compose_series_brent_kung(fmprb_poly_t res,
+                    const fmprb_poly_t poly1,
+                    const fmprb_poly_t poly2, long n, long prec)
+{
+    long len1 = poly1->length;
+    long len2 = poly2->length;
+    long lenr;
+
+    if (len2 != 0 && !fmprb_is_zero(poly2->coeffs))
+    {
+        printf("exception: compose_series: inner "
+                "polynomial must have zero constant term\n");
+        abort();
+    }
+
+    if (len1 == 0 || n == 0)
+    {
+        fmprb_poly_zero(res);
+        return;
+    }
+
+    if (len2 == 0 || len1 == 1)
+    {
+        fmprb_poly_set_fmprb(res, poly1->coeffs);
+        return;
+    }
+
+    lenr = FLINT_MIN((len1 - 1) * (len2 - 1) + 1, n);
+    len1 = FLINT_MIN(len1, lenr);
+    len2 = FLINT_MIN(len2, lenr);
+
+    if ((res != poly1) && (res != poly2))
+    {
+        fmprb_poly_fit_length(res, lenr);
+        _fmprb_poly_compose_series_brent_kung(res->coeffs, poly1->coeffs, len1,
+                                        poly2->coeffs, len2, lenr, prec);
+        _fmprb_poly_set_length(res, lenr);
+        _fmprb_poly_normalise(res);
+    }
+    else
+    {
+        fmprb_poly_t t;
+        fmprb_poly_init2(t, lenr);
+        _fmprb_poly_compose_series_brent_kung(t->coeffs, poly1->coeffs, len1,
+                                        poly2->coeffs, len2, lenr, prec);
+        _fmprb_poly_set_length(t, lenr);
+        _fmprb_poly_normalise(t);
+        fmprb_poly_swap(res, t);
+        fmprb_poly_clear(t);
+    }
+}

fmprb_poly/compose_series_horner.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 "fmprb_poly.h"
+
+void
+_fmprb_poly_compose_series_horner(fmprb_struct * res, const fmprb_struct * poly1, long len1,
+                            const fmprb_struct * poly2, long len2, long n, long prec)
+{
+    if (n == 1)
+    {
+        fmprb_set(res, poly1);
+    }
+    else
+    {
+        long i = len1 - 1;
+        long lenr;
+
+        fmprb_struct * t = _fmprb_vec_init(n);
+
+        lenr = len2;
+        _fmprb_vec_scalar_mul(res, poly2, len2, poly1 + i, prec);
+        i--;
+        fmprb_add(res, res, poly1 + i, prec);
+
+        while (i > 0)
+        {
+            i--;
+            if (lenr + len2 - 1 < n)
+            {
+                _fmprb_poly_mul(t, res, lenr, poly2, len2, prec);
+                lenr = lenr + len2 - 1;
+            }
+            else
+            {
+                _fmprb_poly_mullow(t, res, lenr, poly2, len2, n, prec);
+                lenr = n;
+            }
+            _fmprb_poly_add(res, t, lenr, poly1 + i, 1, prec);
+        }
+
+        _fmprb_vec_zero(res + lenr, n - lenr);
+        _fmprb_vec_clear(t, n);
+    }
+}
+
+void
+fmprb_poly_compose_series_horner(fmprb_poly_t res,
+                    const fmprb_poly_t poly1,
+                    const fmprb_poly_t poly2, long n, long prec)
+{
+    long len1 = poly1->length;
+    long len2 = poly2->length;
+    long lenr;
+
+    if (len2 != 0 && !fmprb_is_zero(poly2->coeffs))
+    {
+        printf("exception: compose_series: inner "
+                "polynomial must have zero constant term\n");
+        abort();
+    }
+
+    if (len1 == 0 || n == 0)
+    {
+        fmprb_poly_zero(res);
+        return;
+    }
+
+    if (len2 == 0 || len1 == 1)
+    {
+        fmprb_poly_set_fmprb(res, poly1->coeffs);
+        return;
+    }
+
+    lenr = FLINT_MIN((len1 - 1) * (len2 - 1) + 1, n);
+    len1 = FLINT_MIN(len1, lenr);
+    len2 = FLINT_MIN(len2, lenr);
+
+    if ((res != poly1) && (res != poly2))
+    {
+        fmprb_poly_fit_length(res, lenr);
+        _fmprb_poly_compose_series_horner(res->coeffs, poly1->coeffs, len1,
+                                        poly2->coeffs, len2, lenr, prec);
+        _fmprb_poly_set_length(res, lenr);
+        _fmprb_poly_normalise(res);
+    }
+    else
+    {
+        fmprb_poly_t t;
+        fmprb_poly_init2(t, lenr);
+        _fmprb_poly_compose_series_horner(t->coeffs, poly1->coeffs, len1,
+                                        poly2->coeffs, len2, lenr, prec);
+        _fmprb_poly_set_length(t, lenr);
+        _fmprb_poly_normalise(t);
+        fmprb_poly_swap(res, t);
+        fmprb_poly_clear(t);
+    }
+}

fmprb_poly/test/t-compose_series.c

@@ -0,0 +1,116 @@
+/*=============================================================================
+
+    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 "fmprb_poly.h"
+
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("compose_series....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 3000; iter++)
+    {
+        long qbits1, qbits2, rbits1, rbits2, rbits3, n;
+        fmpq_poly_t A, B, C;
+        fmprb_poly_t a, b, c, d;
+
+        qbits1 = 2 + n_randint(state, 200);
+        qbits2 = 2 + n_randint(state, 200);
+        rbits1 = 2 + n_randint(state, 200);
+        rbits2 = 2 + n_randint(state, 200);
+        rbits3 = 2 + n_randint(state, 200);
+        n = 2 + n_randint(state, 25);
+
+        fmpq_poly_init(A);
+        fmpq_poly_init(B);
+        fmpq_poly_init(C);
+
+        fmprb_poly_init(a);
+        fmprb_poly_init(b);
+        fmprb_poly_init(c);
+        fmprb_poly_init(d);
+
+        fmpq_poly_randtest(A, state, 1 + n_randint(state, 25), qbits1);
+        fmpq_poly_randtest(B, state, 1 + n_randint(state, 25), qbits2);
+        fmpq_poly_set_coeff_ui(B, 0, 0);
+        fmpq_poly_compose_series(C, A, B, n);
+
+        fmprb_poly_set_fmpq_poly(a, A, rbits1);
+        fmprb_poly_set_fmpq_poly(b, B, rbits2);
+        fmprb_poly_compose_series(c, a, b, n, rbits3);
+
+        if (!fmprb_poly_contains_fmpq_poly(c, C))
+        {
+            printf("FAIL\n\n");
+            printf("n = %ld, bits3 = %ld\n", n, rbits3);
+
+            printf("A = "); fmpq_poly_print(A); printf("\n\n");
+            printf("B = "); fmpq_poly_print(B); printf("\n\n");
+            printf("C = "); fmpq_poly_print(C); printf("\n\n");
+
+            printf("a = "); fmprb_poly_printd(a, 15); printf("\n\n");
+            printf("b = "); fmprb_poly_printd(b, 15); printf("\n\n");
+            printf("c = "); fmprb_poly_printd(c, 15); printf("\n\n");
+
+            abort();
+        }
+
+        fmprb_poly_set(d, a);
+        fmprb_poly_compose_series(d, d, b, n, rbits3);
+        if (!fmprb_poly_equal(d, c))
+        {
+            printf("FAIL (aliasing 1)\n\n");
+            abort();
+        }
+
+        fmprb_poly_set(d, b);
+        fmprb_poly_compose_series(d, a, d, n, rbits3);
+        if (!fmprb_poly_equal(d, c))
+        {
+            printf("FAIL (aliasing 2)\n\n");
+            abort();
+        }
+
+        fmpq_poly_clear(A);
+        fmpq_poly_clear(B);
+        fmpq_poly_clear(C);
+
+        fmprb_poly_clear(a);
+        fmprb_poly_clear(b);
+        fmprb_poly_clear(c);
+        fmprb_poly_clear(d);
+    }
+
+    flint_randclear(state);
+    _fmpz_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}

fmprb_poly/test/t-compose_series_brent_kung.c

@@ -0,0 +1,116 @@
+/*=============================================================================
+
+    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 "fmprb_poly.h"
+
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("compose_series_brent_kung....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 3000; iter++)
+    {
+        long qbits1, qbits2, rbits1, rbits2, rbits3, n;
+        fmpq_poly_t A, B, C;
+        fmprb_poly_t a, b, c, d;
+
+        qbits1 = 2 + n_randint(state, 200);
+        qbits2 = 2 + n_randint(state, 200);
+        rbits1 = 2 + n_randint(state, 200);
+        rbits2 = 2 + n_randint(state, 200);
+        rbits3 = 2 + n_randint(state, 200);
+        n = 2 + n_randint(state, 25);
+
+        fmpq_poly_init(A);
+        fmpq_poly_init(B);
+        fmpq_poly_init(C);
+
+        fmprb_poly_init(a);
+        fmprb_poly_init(b);
+        fmprb_poly_init(c);
+        fmprb_poly_init(d);
+
+        fmpq_poly_randtest(A, state, 1 + n_randint(state, 25), qbits1);
+        fmpq_poly_randtest(B, state, 1 + n_randint(state, 25), qbits2);
+        fmpq_poly_set_coeff_ui(B, 0, 0);
+        fmpq_poly_compose_series(C, A, B, n);
+
+        fmprb_poly_set_fmpq_poly(a, A, rbits1);
+        fmprb_poly_set_fmpq_poly(b, B, rbits2);
+        fmprb_poly_compose_series_brent_kung(c, a, b, n, rbits3);
+
+        if (!fmprb_poly_contains_fmpq_poly(c, C))
+        {
+            printf("FAIL\n\n");
+            printf("n = %ld, bits3 = %ld\n", n, rbits3);
+
+            printf("A = "); fmpq_poly_print(A); printf("\n\n");
+            printf("B = "); fmpq_poly_print(B); printf("\n\n");
+            printf("C = "); fmpq_poly_print(C); printf("\n\n");
+
+            printf("a = "); fmprb_poly_printd(a, 15); printf("\n\n");
+            printf("b = "); fmprb_poly_printd(b, 15); printf("\n\n");
+            printf("c = "); fmprb_poly_printd(c, 15); printf("\n\n");
+
+            abort();
+        }
+
+        fmprb_poly_set(d, a);
+        fmprb_poly_compose_series_brent_kung(d, d, b, n, rbits3);
+        if (!fmprb_poly_equal(d, c))
+        {
+            printf("FAIL (aliasing 1)\n\n");
+            abort();
+        }
+
+        fmprb_poly_set(d, b);
+        fmprb_poly_compose_series_brent_kung(d, a, d, n, rbits3);
+        if (!fmprb_poly_equal(d, c))
+        {
+            printf("FAIL (aliasing 2)\n\n");
+            abort();
+        }
+
+        fmpq_poly_clear(A);
+        fmpq_poly_clear(B);
+        fmpq_poly_clear(C);
+
+        fmprb_poly_clear(a);
+        fmprb_poly_clear(b);
+        fmprb_poly_clear(c);
+        fmprb_poly_clear(d);
+    }
+
+    flint_randclear(state);
+    _fmpz_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}

fmprb_poly/test/t-compose_series_horner.c

@@ -0,0 +1,116 @@
+/*=============================================================================
+
+    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 "fmprb_poly.h"
+
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("compose_series_horner....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 3000; iter++)
+    {
+        long qbits1, qbits2, rbits1, rbits2, rbits3, n;
+        fmpq_poly_t A, B, C;
+        fmprb_poly_t a, b, c, d;
+
+        qbits1 = 2 + n_randint(state, 200);
+        qbits2 = 2 + n_randint(state, 200);
+        rbits1 = 2 + n_randint(state, 200);
+        rbits2 = 2 + n_randint(state, 200);
+        rbits3 = 2 + n_randint(state, 200);
+        n = 2 + n_randint(state, 25);
+
+        fmpq_poly_init(A);
+        fmpq_poly_init(B);
+        fmpq_poly_init(C);
+
+        fmprb_poly_init(a);
+        fmprb_poly_init(b);
+        fmprb_poly_init(c);
+        fmprb_poly_init(d);
+
+        fmpq_poly_randtest(A, state, 1 + n_randint(state, 25), qbits1);
+        fmpq_poly_randtest(B, state, 1 + n_randint(state, 25), qbits2);
+        fmpq_poly_set_coeff_ui(B, 0, 0);
+        fmpq_poly_compose_series(C, A, B, n);
+
+        fmprb_poly_set_fmpq_poly(a, A, rbits1);
+        fmprb_poly_set_fmpq_poly(b, B, rbits2);
+        fmprb_poly_compose_series_horner(c, a, b, n, rbits3);
+
+        if (!fmprb_poly_contains_fmpq_poly(c, C))
+        {
+            printf("FAIL\n\n");
+            printf("n = %ld, bits3 = %ld\n", n, rbits3);
+
+            printf("A = "); fmpq_poly_print(A); printf("\n\n");
+            printf("B = "); fmpq_poly_print(B); printf("\n\n");
+            printf("C = "); fmpq_poly_print(C); printf("\n\n");
+
+            printf("a = "); fmprb_poly_printd(a, 15); printf("\n\n");
+            printf("b = "); fmprb_poly_printd(b, 15); printf("\n\n");
+            printf("c = "); fmprb_poly_printd(c, 15); printf("\n\n");
+
+            abort();
+        }
+
+        fmprb_poly_set(d, a);
+        fmprb_poly_compose_series_horner(d, d, b, n, rbits3);
+        if (!fmprb_poly_equal(d, c))
+        {
+            printf("FAIL (aliasing 1)\n\n");
+            abort();
+        }
+
+        fmprb_poly_set(d, b);
+        fmprb_poly_compose_series_horner(d, a, d, n, rbits3);
+        if (!fmprb_poly_equal(d, c))
+        {
+            printf("FAIL (aliasing 2)\n\n");
+            abort();
+        }
+
+        fmpq_poly_clear(A);
+        fmpq_poly_clear(B);
+        fmpq_poly_clear(C);
+
+        fmprb_poly_clear(a);
+        fmprb_poly_clear(b);
+        fmprb_poly_clear(c);
+        fmprb_poly_clear(d);
+    }
+
+    flint_randclear(state);
+    _fmpz_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}