implement fast Taylor shift and use in special-form polynomial composition

acb_poly.h

@@ -322,6 +322,18 @@ void _acb_poly_div_root(acb_ptr Q, acb_t R, acb_srcptr A,
 
 /* Composition */
 
+void _acb_poly_taylor_shift_horner(acb_ptr poly, const acb_t c, slong n, slong prec);
+
+void acb_poly_taylor_shift_horner(acb_poly_t g, const acb_poly_t f, const acb_t c, slong prec);
+
+void _acb_poly_taylor_shift_divconquer(acb_ptr poly, const acb_t c, slong n, slong prec);
+
+void acb_poly_taylor_shift_divconquer(acb_poly_t g, const acb_poly_t f, const acb_t c, slong prec);
+
+void _acb_poly_taylor_shift(acb_ptr poly, const acb_t c, slong n, slong prec);
+
+void acb_poly_taylor_shift(acb_poly_t g, const acb_poly_t f, const acb_t c, slong prec);
+
 void _acb_poly_compose(acb_ptr res,
     acb_srcptr poly1, slong len1,
     acb_srcptr poly2, slong len2, slong prec);

acb_poly/compose.c

@@ -21,40 +21,76 @@
 
     Copyright (C) 2010 William Hart
     Copyright (C) 2012 Sebastian Pancratz
-    Copyright (C) 2012 Fredrik Johansson
+    Copyright (C) 2012, 2016 Fredrik Johansson
 
 ******************************************************************************/
 
 #include "acb_poly.h"
 
+/* compose by poly2 = a*x^n + c, no aliasing; n >= 1 */
+void
+_acb_poly_compose_axnc(acb_ptr res, acb_srcptr poly1, slong len1,
+    const acb_t c, const acb_t a, slong n, slong prec)
+{
+    slong i;
+
+    _acb_vec_set_round(res, poly1, len1, prec);
+    /* shift by c (c = 0 case will be fast) */
+    _acb_poly_taylor_shift(res, c, len1, prec);
+
+    /* multiply by powers of a */
+    if (!acb_is_one(a))
+    {
+        if (acb_equal_si(a, -1))
+        {
+            for (i = 1; i < len1; i += 2)
+                acb_neg(res + i, res + i);
+        }
+        else if (len1 == 2)
+        {
+            acb_mul(res + 1, res + 1, a, prec);
+        }
+        else
+        {
+            acb_t t;
+            acb_init(t);
+            acb_set(t, a);
+
+            for (i = 1; i < len1; i++)
+            {
+                acb_mul(res + i, res + i, t, prec);
+                if (i + 1 < len1)
+                    acb_mul(t, t, a, prec);
+            }
+
+            acb_clear(t);
+        }
+    }
+
+    /* stretch */
+    for (i = len1 - 1; i >= 1 && n > 1; i--)
+    {
+        acb_swap(res + i * n, res + i);
+        _acb_vec_zero(res + (i - 1) * n + 1, n - 1);
+    }
+}
+
 void
 _acb_poly_compose(acb_ptr res,
     acb_srcptr poly1, slong len1,
     acb_srcptr poly2, slong len2, slong prec)
 {
-    if (len2 == 1)
+    if (len1 == 1)
+    {
+        acb_set_round(res, poly1, prec);
+    }
+    else if (len2 == 1)
     {
         _acb_poly_evaluate(res, poly1, len1, poly2, prec);
     }
-    else if (_acb_vec_is_zero(poly2, len2 - 1)) /* poly2 is a monomial */
+    else if (_acb_vec_is_zero(poly2 + 1, len2 - 2))
     {
-        slong i;
-        acb_t t;
-
-        acb_init(t);
-        acb_set(t, poly2 + len2 - 1);
-        acb_set_round(res, poly1, prec);
-
-        for (i = 1; i < len1; i++)
-        {
-            acb_mul(res + i * (len2 - 1), poly1 + i, t, prec);
-            if (i + 1 < len1)
-                acb_mul(t, t, poly2 + len2 - 1, prec);
-
-            _acb_vec_zero(res + (i - 1) * (len2 - 1) + 1, len2 - 2);
-        }
-
-        acb_clear(t);
+        _acb_poly_compose_axnc(res, poly1, len1, poly2, poly2 + len2 - 1, len2 - 1, prec);
     }
     else if (len1 <= 7)
     {
@@ -104,3 +140,4 @@ void acb_poly_compose(acb_poly_t res,
         _acb_poly_normalise(res);
     }
 }
+

acb_poly/taylor_shift.c

@@ -0,0 +1,47 @@
+/*=============================================================================
+
+    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) 2016 Fredrik Johansson
+
+******************************************************************************/
+
+#include "acb_poly.h"
+
+void
+_acb_poly_taylor_shift(acb_ptr poly, const acb_t c, slong n, slong prec)
+{
+    if (n <= 30 || (n <= 1000 && acb_bits(c) == 1 && n < 30 + 3 * sqrt(prec))
+            || (n <= 100 && acb_bits(c) < 0.01 * prec))
+        _acb_poly_taylor_shift_horner(poly, c, n, prec);
+    else
+        _acb_poly_taylor_shift_divconquer(poly, c, n, prec);
+}
+
+void
+acb_poly_taylor_shift(acb_poly_t g, const acb_poly_t f,
+    const acb_t c, slong prec)
+{
+    if (f != g)
+        acb_poly_set_round(g, f, prec);
+
+    _acb_poly_taylor_shift(g->coeffs, c, g->length, prec);
+}
+

acb_poly/taylor_shift_divconquer.c

@@ -0,0 +1,59 @@
+/*=============================================================================
+
+    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) 2016 Fredrik Johansson
+
+******************************************************************************/
+
+#include "acb_poly.h"
+
+void
+_acb_poly_taylor_shift_divconquer(acb_ptr poly, const acb_t c, slong len, slong prec)
+{
+    acb_struct d[2];
+
+    if (len <= 1 || acb_is_zero(c))
+        return;
+
+    if (len == 2)
+    {
+        acb_addmul(poly, poly + 1, c, prec);
+        return;
+    }
+
+    d[0] = *c;
+
+    acb_init(d + 1);
+    acb_one(d + 1); /* no need to free */
+
+    _acb_poly_compose_divconquer(poly, poly, len, d, 2, prec);
+}
+
+void
+acb_poly_taylor_shift_divconquer(acb_poly_t g, const acb_poly_t f,
+    const acb_t c, slong prec)
+{
+    if (f != g)
+        acb_poly_set_round(g, f, prec);
+
+    _acb_poly_taylor_shift_divconquer(g->coeffs, c, g->length, prec);
+}
+

acb_poly/taylor_shift_horner.c

@@ -0,0 +1,62 @@
+/*=============================================================================
+
+    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) 2016 Fredrik Johansson
+
+******************************************************************************/
+
+#include "acb_poly.h"
+
+void
+_acb_poly_taylor_shift_horner(acb_ptr poly, const acb_t c, slong n, slong prec)
+{
+    slong i, j;
+
+    if (acb_is_one(c))
+    {
+        for (i = n - 2; i >= 0; i--)
+            for (j = i; j < n - 1; j++)
+                acb_add(poly + j, poly + j, poly + j + 1, prec);
+    }
+    else if (acb_equal_si(c, -1))
+    {
+        for (i = n - 2; i >= 0; i--)
+            for (j = i; j < n - 1; j++)
+                acb_sub(poly + j, poly + j, poly + j + 1, prec);
+    }
+    else if (!acb_is_zero(c))
+    {
+        for (i = n - 2; i >= 0; i--)
+            for (j = i; j < n - 1; j++)
+                acb_addmul(poly + j, poly + j + 1, c, prec);
+    }
+}
+
+void
+acb_poly_taylor_shift_horner(acb_poly_t g, const acb_poly_t f,
+    const acb_t c, slong prec)
+{
+    if (f != g)
+        acb_poly_set_round(g, f, prec);
+
+    _acb_poly_taylor_shift_horner(g->coeffs, c, g->length, prec);
+}
+

acb_poly/test/t-taylor_shift.c

@@ -0,0 +1,99 @@
+/*=============================================================================
+
+    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) 2016 Fredrik Johansson
+
+******************************************************************************/
+
+#include "acb_poly.h"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("taylor_shift....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 1000; iter++)
+    {
+        slong prec1, prec2;
+        acb_poly_t f, g;
+        acb_t c, d, e;
+
+        prec1 = 2 + n_randint(state, 500);
+        prec2 = 2 + n_randint(state, 500);
+
+        acb_poly_init(f);
+        acb_poly_init(g);
+
+        acb_init(c);
+        acb_init(d);
+        acb_init(e);
+
+        acb_poly_randtest(f, state, 1 + n_randint(state, 40), 1 + n_randint(state, 500), 10);
+        acb_poly_randtest(g, state, 1 + n_randint(state, 20), 1 + n_randint(state, 500), 10);
+
+        if (n_randint(state, 2))
+            acb_set_si(c, n_randint(state, 5) - 2);
+        else
+            acb_randtest(c, state, 1 + n_randint(state, 500), 1 + n_randint(state, 100));
+
+        if (n_randint(state, 2))
+            acb_set_si(d, n_randint(state, 5) - 2);
+        else
+            acb_randtest(d, state, 1 + n_randint(state, 500), 1 + n_randint(state, 100));
+
+        acb_add(e, c, d, prec1);
+
+        /* check f(x+c)(x+d) = f(x+c+d) */
+        acb_poly_taylor_shift(g, f, e, prec2);
+        acb_poly_taylor_shift(f, f, c, prec1);
+        acb_poly_taylor_shift(f, f, d, prec1);
+
+        if (!acb_poly_overlaps(f, g))
+        {
+            flint_printf("FAIL\n\n");
+
+            flint_printf("c = "); acb_printd(c, 15); flint_printf("\n\n");
+            flint_printf("d = "); acb_printd(d, 15); flint_printf("\n\n");
+
+            flint_printf("f = "); acb_poly_printd(f, 15); flint_printf("\n\n");
+            flint_printf("g = "); acb_poly_printd(g, 15); flint_printf("\n\n");
+
+            abort();
+        }
+
+        acb_poly_clear(f);
+        acb_poly_clear(g);
+
+        acb_clear(c);
+        acb_clear(d);
+        acb_clear(e);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}

acb_poly/test/t-taylor_shift_divconquer.c

@@ -0,0 +1,99 @@
+/*=============================================================================
+
+    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) 2016 Fredrik Johansson
+
+******************************************************************************/
+
+#include "acb_poly.h"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("taylor_shift_divconquer....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 1000; iter++)
+    {
+        slong prec1, prec2;
+        acb_poly_t f, g;
+        acb_t c, d, e;
+
+        prec1 = 2 + n_randint(state, 500);
+        prec2 = 2 + n_randint(state, 500);
+
+        acb_poly_init(f);
+        acb_poly_init(g);
+
+        acb_init(c);
+        acb_init(d);
+        acb_init(e);
+
+        acb_poly_randtest(f, state, 1 + n_randint(state, 40), 1 + n_randint(state, 500), 10);
+        acb_poly_randtest(g, state, 1 + n_randint(state, 20), 1 + n_randint(state, 500), 10);
+
+        if (n_randint(state, 2))
+            acb_set_si(c, n_randint(state, 5) - 2);
+        else
+            acb_randtest(c, state, 1 + n_randint(state, 500), 1 + n_randint(state, 100));
+
+        if (n_randint(state, 2))
+            acb_set_si(d, n_randint(state, 5) - 2);
+        else
+            acb_randtest(d, state, 1 + n_randint(state, 500), 1 + n_randint(state, 100));
+
+        acb_add(e, c, d, prec1);
+
+        /* check f(x+c)(x+d) = f(x+c+d) */
+        acb_poly_taylor_shift_divconquer(g, f, e, prec2);
+        acb_poly_taylor_shift_divconquer(f, f, c, prec1);
+        acb_poly_taylor_shift_divconquer(f, f, d, prec1);
+
+        if (!acb_poly_overlaps(f, g))
+        {
+            flint_printf("FAIL\n\n");
+
+            flint_printf("c = "); acb_printd(c, 15); flint_printf("\n\n");
+            flint_printf("d = "); acb_printd(d, 15); flint_printf("\n\n");
+
+            flint_printf("f = "); acb_poly_printd(f, 15); flint_printf("\n\n");
+            flint_printf("g = "); acb_poly_printd(g, 15); flint_printf("\n\n");
+
+            abort();
+        }
+
+        acb_poly_clear(f);
+        acb_poly_clear(g);
+
+        acb_clear(c);
+        acb_clear(d);
+        acb_clear(e);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}

acb_poly/test/t-taylor_shift_horner.c

@@ -0,0 +1,99 @@
+/*=============================================================================
+
+    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) 2016 Fredrik Johansson
+
+******************************************************************************/
+
+#include "acb_poly.h"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("taylor_shift_horner....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 1000; iter++)
+    {
+        slong prec1, prec2;
+        acb_poly_t f, g;
+        acb_t c, d, e;
+
+        prec1 = 2 + n_randint(state, 500);
+        prec2 = 2 + n_randint(state, 500);
+
+        acb_poly_init(f);
+        acb_poly_init(g);
+
+        acb_init(c);
+        acb_init(d);
+        acb_init(e);
+
+        acb_poly_randtest(f, state, 1 + n_randint(state, 40), 1 + n_randint(state, 500), 10);
+        acb_poly_randtest(g, state, 1 + n_randint(state, 20), 1 + n_randint(state, 500), 10);
+
+        if (n_randint(state, 2))
+            acb_set_si(c, n_randint(state, 5) - 2);
+        else
+            acb_randtest(c, state, 1 + n_randint(state, 500), 1 + n_randint(state, 100));
+
+        if (n_randint(state, 2))
+            acb_set_si(d, n_randint(state, 5) - 2);
+        else
+            acb_randtest(d, state, 1 + n_randint(state, 500), 1 + n_randint(state, 100));
+
+        acb_add(e, c, d, prec1);
+
+        /* check f(x+c)(x+d) = f(x+c+d) */
+        acb_poly_taylor_shift_horner(g, f, e, prec2);
+        acb_poly_taylor_shift_horner(f, f, c, prec1);
+        acb_poly_taylor_shift_horner(f, f, d, prec1);
+
+        if (!acb_poly_overlaps(f, g))
+        {
+            flint_printf("FAIL\n\n");
+
+            flint_printf("c = "); acb_printd(c, 15); flint_printf("\n\n");
+            flint_printf("d = "); acb_printd(d, 15); flint_printf("\n\n");
+
+            flint_printf("f = "); acb_poly_printd(f, 15); flint_printf("\n\n");
+            flint_printf("g = "); acb_poly_printd(g, 15); flint_printf("\n\n");
+
+            abort();
+        }
+
+        acb_poly_clear(f);
+        acb_poly_clear(g);
+
+        acb_clear(c);
+        acb_clear(d);
+        acb_clear(e);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}

arb_poly.h

@@ -310,6 +310,18 @@ void _arb_poly_tree_build(arb_ptr * tree, arb_srcptr roots, slong len, slong pre
 
 /* Composition */
 
+void _arb_poly_taylor_shift_horner(arb_ptr poly, const arb_t c, slong n, slong prec);
+
+void arb_poly_taylor_shift_horner(arb_poly_t g, const arb_poly_t f, const arb_t c, slong prec);
+
+void _arb_poly_taylor_shift_divconquer(arb_ptr poly, const arb_t c, slong n, slong prec);
+
+void arb_poly_taylor_shift_divconquer(arb_poly_t g, const arb_poly_t f, const arb_t c, slong prec);
+
+void _arb_poly_taylor_shift(arb_ptr poly, const arb_t c, slong n, slong prec);
+
+void arb_poly_taylor_shift(arb_poly_t g, const arb_poly_t f, const arb_t c, slong prec);
+
 void _arb_poly_compose(arb_ptr res,
     arb_srcptr poly1, slong len1,
     arb_srcptr poly2, slong len2, slong prec);

arb_poly/compose.c

@@ -21,40 +21,76 @@
 
     Copyright (C) 2010 William Hart
     Copyright (C) 2012 Sebastian Pancratz
-    Copyright (C) 2012 Fredrik Johansson
+    Copyright (C) 2012, 2016 Fredrik Johansson
 
 ******************************************************************************/
 
 #include "arb_poly.h"
 
+/* compose by poly2 = a*x^n + c, no aliasing; n >= 1 */
+void
+_arb_poly_compose_axnc(arb_ptr res, arb_srcptr poly1, slong len1,
+    const arb_t c, const arb_t a, slong n, slong prec)
+{
+    slong i;
+
+    _arb_vec_set_round(res, poly1, len1, prec);
+    /* shift by c (c = 0 case will be fast) */
+    _arb_poly_taylor_shift(res, c, len1, prec);
+
+    /* multiply by powers of a */
+    if (!arb_is_one(a))
+    {
+        if (arb_equal_si(a, -1))
+        {
+            for (i = 1; i < len1; i += 2)
+                arb_neg(res + i, res + i);
+        }
+        else if (len1 == 2)
+        {
+            arb_mul(res + 1, res + 1, a, prec);
+        }
+        else
+        {
+            arb_t t;
+            arb_init(t);
+            arb_set(t, a);
+
+            for (i = 1; i < len1; i++)
+            {
+                arb_mul(res + i, res + i, t, prec);
+                if (i + 1 < len1)
+                    arb_mul(t, t, a, prec);
+            }
+
+            arb_clear(t);
+        }
+    }
+
+    /* stretch */
+    for (i = len1 - 1; i >= 1 && n > 1; i--)
+    {
+        arb_swap(res + i * n, res + i);
+        _arb_vec_zero(res + (i - 1) * n + 1, n - 1);
+    }
+}
+
 void
 _arb_poly_compose(arb_ptr res,
     arb_srcptr poly1, slong len1,
     arb_srcptr poly2, slong len2, slong prec)
 {
-    if (len2 == 1)
+    if (len1 == 1)
+    {
+        arb_set_round(res, poly1, prec);
+    }
+    else if (len2 == 1)
     {
         _arb_poly_evaluate(res, poly1, len1, poly2, prec);
     }
-    else if (_arb_vec_is_zero(poly2, len2 - 1)) /* poly2 is a monomial */
+    else if (_arb_vec_is_zero(poly2 + 1, len2 - 2))
     {
-        slong i;
-        arb_t t;
-
-        arb_init(t);
-        arb_set(t, poly2 + len2 - 1);
-        arb_set_round(res, poly1, prec);
-
-        for (i = 1; i < len1; i++)
-        {
-            arb_mul(res + i * (len2 - 1), poly1 + i, t, prec);
-            if (i + 1 < len1)
-                arb_mul(t, t, poly2 + len2 - 1, prec);
-
-            _arb_vec_zero(res + (i - 1) * (len2 - 1) + 1, len2 - 2);
-        }
-
-        arb_clear(t);
+        _arb_poly_compose_axnc(res, poly1, len1, poly2, poly2 + len2 - 1, len2 - 1, prec);
     }
     else if (len1 <= 7)
     {
@@ -104,3 +140,4 @@ void arb_poly_compose(arb_poly_t res,
         _arb_poly_normalise(res);
     }
 }
+

arb_poly/taylor_shift.c

@@ -0,0 +1,47 @@
+/*=============================================================================
+
+    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) 2016 Fredrik Johansson
+
+******************************************************************************/
+
+#include "arb_poly.h"
+
+void
+_arb_poly_taylor_shift(arb_ptr poly, const arb_t c, slong n, slong prec)
+{
+    if (n <= 30 || (n <= 1000 && arb_bits(c) == 1 && n < 30 + 3 * sqrt(prec))
+            || (n <= 100 && arb_bits(c) < 0.01 * prec))
+        _arb_poly_taylor_shift_horner(poly, c, n, prec);
+    else
+        _arb_poly_taylor_shift_divconquer(poly, c, n, prec);
+}
+
+void
+arb_poly_taylor_shift(arb_poly_t g, const arb_poly_t f,
+    const arb_t c, slong prec)
+{
+    if (f != g)
+        arb_poly_set_round(g, f, prec);
+
+    _arb_poly_taylor_shift(g->coeffs, c, g->length, prec);
+}
+

arb_poly/taylor_shift_divconquer.c

@@ -0,0 +1,59 @@
+/*=============================================================================
+
+    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) 2016 Fredrik Johansson
+
+******************************************************************************/
+
+#include "arb_poly.h"
+
+void
+_arb_poly_taylor_shift_divconquer(arb_ptr poly, const arb_t c, slong len, slong prec)
+{
+    arb_struct d[2];
+
+    if (len <= 1 || arb_is_zero(c))
+        return;
+
+    if (len == 2)
+    {
+        arb_addmul(poly, poly + 1, c, prec);
+        return;
+    }
+
+    d[0] = *c;
+
+    arb_init(d + 1);
+    arb_one(d + 1); /* no need to free */
+
+    _arb_poly_compose_divconquer(poly, poly, len, d, 2, prec);
+}
+
+void
+arb_poly_taylor_shift_divconquer(arb_poly_t g, const arb_poly_t f,
+    const arb_t c, slong prec)
+{
+    if (f != g)
+        arb_poly_set_round(g, f, prec);
+
+    _arb_poly_taylor_shift_divconquer(g->coeffs, c, g->length, prec);
+}
+

arb_poly/taylor_shift_horner.c

@@ -0,0 +1,62 @@
+/*=============================================================================
+
+    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) 2016 Fredrik Johansson
+
+******************************************************************************/
+
+#include "arb_poly.h"
+
+void
+_arb_poly_taylor_shift_horner(arb_ptr poly, const arb_t c, slong n, slong prec)
+{
+    slong i, j;
+
+    if (arb_is_one(c))
+    {
+        for (i = n - 2; i >= 0; i--)
+            for (j = i; j < n - 1; j++)
+                arb_add(poly + j, poly + j, poly + j + 1, prec);
+    }
+    else if (arb_equal_si(c, -1))
+    {
+        for (i = n - 2; i >= 0; i--)
+            for (j = i; j < n - 1; j++)
+                arb_sub(poly + j, poly + j, poly + j + 1, prec);
+    }
+    else if (!arb_is_zero(c))
+    {
+        for (i = n - 2; i >= 0; i--)
+            for (j = i; j < n - 1; j++)
+                arb_addmul(poly + j, poly + j + 1, c, prec);
+    }
+}
+
+void
+arb_poly_taylor_shift_horner(arb_poly_t g, const arb_poly_t f,
+    const arb_t c, slong prec)
+{
+    if (f != g)
+        arb_poly_set_round(g, f, prec);
+
+    _arb_poly_taylor_shift_horner(g->coeffs, c, g->length, prec);
+}
+

arb_poly/test/t-taylor_shift.c

@@ -0,0 +1,99 @@
+/*=============================================================================
+
+    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) 2016 Fredrik Johansson
+
+******************************************************************************/
+
+#include "arb_poly.h"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("taylor_shift....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 2000; iter++)
+    {
+        slong prec1, prec2;
+        arb_poly_t f, g;
+        arb_t c, d, e;
+
+        prec1 = 2 + n_randint(state, 500);
+        prec2 = 2 + n_randint(state, 500);
+
+        arb_poly_init(f);
+        arb_poly_init(g);
+
+        arb_init(c);
+        arb_init(d);
+        arb_init(e);
+
+        arb_poly_randtest(f, state, 1 + n_randint(state, 40), 1 + n_randint(state, 500), 10);
+        arb_poly_randtest(g, state, 1 + n_randint(state, 20), 1 + n_randint(state, 500), 10);
+
+        if (n_randint(state, 2))
+            arb_set_si(c, n_randint(state, 5) - 2);
+        else
+            arb_randtest(c, state, 1 + n_randint(state, 500), 1 + n_randint(state, 100));
+
+        if (n_randint(state, 2))
+            arb_set_si(d, n_randint(state, 5) - 2);
+        else
+            arb_randtest(d, state, 1 + n_randint(state, 500), 1 + n_randint(state, 100));
+
+        arb_add(e, c, d, prec1);
+
+        /* check f(x+c)(x+d) = f(x+c+d) */
+        arb_poly_taylor_shift(g, f, e, prec2);
+        arb_poly_taylor_shift(f, f, c, prec1);
+        arb_poly_taylor_shift(f, f, d, prec1);
+
+        if (!arb_poly_overlaps(f, g))
+        {
+            flint_printf("FAIL\n\n");
+
+            flint_printf("c = "); arb_printd(c, 15); flint_printf("\n\n");
+            flint_printf("d = "); arb_printd(d, 15); flint_printf("\n\n");
+
+            flint_printf("f = "); arb_poly_printd(f, 15); flint_printf("\n\n");
+            flint_printf("g = "); arb_poly_printd(g, 15); flint_printf("\n\n");
+
+            abort();
+        }
+
+        arb_poly_clear(f);
+        arb_poly_clear(g);
+
+        arb_clear(c);
+        arb_clear(d);
+        arb_clear(e);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}

arb_poly/test/t-taylor_shift_divconquer.c

@@ -0,0 +1,99 @@
+/*=============================================================================
+
+    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) 2016 Fredrik Johansson
+
+******************************************************************************/
+
+#include "arb_poly.h"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("taylor_shift_divconquer....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 2000; iter++)
+    {
+        slong prec1, prec2;
+        arb_poly_t f, g;
+        arb_t c, d, e;
+
+        prec1 = 2 + n_randint(state, 500);
+        prec2 = 2 + n_randint(state, 500);
+
+        arb_poly_init(f);
+        arb_poly_init(g);
+
+        arb_init(c);
+        arb_init(d);
+        arb_init(e);
+
+        arb_poly_randtest(f, state, 1 + n_randint(state, 40), 1 + n_randint(state, 500), 10);
+        arb_poly_randtest(g, state, 1 + n_randint(state, 20), 1 + n_randint(state, 500), 10);
+
+        if (n_randint(state, 2))
+            arb_set_si(c, n_randint(state, 5) - 2);
+        else
+            arb_randtest(c, state, 1 + n_randint(state, 500), 1 + n_randint(state, 100));
+
+        if (n_randint(state, 2))
+            arb_set_si(d, n_randint(state, 5) - 2);
+        else
+            arb_randtest(d, state, 1 + n_randint(state, 500), 1 + n_randint(state, 100));
+
+        arb_add(e, c, d, prec1);
+
+        /* check f(x+c)(x+d) = f(x+c+d) */
+        arb_poly_taylor_shift_divconquer(g, f, e, prec2);
+        arb_poly_taylor_shift_divconquer(f, f, c, prec1);
+        arb_poly_taylor_shift_divconquer(f, f, d, prec1);
+
+        if (!arb_poly_overlaps(f, g))
+        {
+            flint_printf("FAIL\n\n");
+
+            flint_printf("c = "); arb_printd(c, 15); flint_printf("\n\n");
+            flint_printf("d = "); arb_printd(d, 15); flint_printf("\n\n");
+
+            flint_printf("f = "); arb_poly_printd(f, 15); flint_printf("\n\n");
+            flint_printf("g = "); arb_poly_printd(g, 15); flint_printf("\n\n");
+
+            abort();
+        }
+
+        arb_poly_clear(f);
+        arb_poly_clear(g);
+
+        arb_clear(c);
+        arb_clear(d);
+        arb_clear(e);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}

arb_poly/test/t-taylor_shift_horner.c

@@ -0,0 +1,99 @@
+/*=============================================================================
+
+    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) 2016 Fredrik Johansson
+
+******************************************************************************/
+
+#include "arb_poly.h"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("taylor_shift_horner....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 2000; iter++)
+    {
+        slong prec1, prec2;
+        arb_poly_t f, g;
+        arb_t c, d, e;
+
+        prec1 = 2 + n_randint(state, 500);
+        prec2 = 2 + n_randint(state, 500);
+
+        arb_poly_init(f);
+        arb_poly_init(g);
+
+        arb_init(c);
+        arb_init(d);
+        arb_init(e);
+
+        arb_poly_randtest(f, state, 1 + n_randint(state, 40), 1 + n_randint(state, 500), 10);
+        arb_poly_randtest(g, state, 1 + n_randint(state, 20), 1 + n_randint(state, 500), 10);
+
+        if (n_randint(state, 2))
+            arb_set_si(c, n_randint(state, 5) - 2);
+        else
+            arb_randtest(c, state, 1 + n_randint(state, 500), 1 + n_randint(state, 100));
+
+        if (n_randint(state, 2))
+            arb_set_si(d, n_randint(state, 5) - 2);
+        else
+            arb_randtest(d, state, 1 + n_randint(state, 500), 1 + n_randint(state, 100));
+
+        arb_add(e, c, d, prec1);
+
+        /* check f(x+c)(x+d) = f(x+c+d) */
+        arb_poly_taylor_shift_horner(g, f, e, prec2);
+        arb_poly_taylor_shift_horner(f, f, c, prec1);
+        arb_poly_taylor_shift_horner(f, f, d, prec1);
+
+        if (!arb_poly_overlaps(f, g))
+        {
+            flint_printf("FAIL\n\n");
+
+            flint_printf("c = "); arb_printd(c, 15); flint_printf("\n\n");
+            flint_printf("d = "); arb_printd(d, 15); flint_printf("\n\n");
+
+            flint_printf("f = "); arb_poly_printd(f, 15); flint_printf("\n\n");
+            flint_printf("g = "); arb_poly_printd(g, 15); flint_printf("\n\n");
+
+            abort();
+        }
+
+        arb_poly_clear(f);
+        arb_poly_clear(g);
+
+        arb_clear(c);
+        arb_clear(d);
+        arb_clear(e);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}

doc/source/acb_poly.rst

@@ -346,6 +346,24 @@ Arithmetic
 Composition
 -------------------------------------------------------------------------------
 
+.. function:: void _acb_poly_taylor_shift_horner(acb_ptr g, const acb_t c, slong n, slong prec)
+
+.. function:: void acb_poly_taylor_shift_horner(acb_poly_t g, const acb_poly_t f, const acb_t c, slong prec)
+
+.. function:: void _acb_poly_taylor_shift_divconquer(acb_ptr g, const acb_t c, slong n, slong prec)
+
+.. function:: void acb_poly_taylor_shift_divconquer(acb_poly_t g, const acb_poly_t f, const acb_t c, slong prec)
+
+.. function:: void _acb_poly_taylor_shift(acb_ptr g, const acb_t c, slong n, slong prec)
+
+.. function:: void acb_poly_taylor_shift(acb_poly_t g, const acb_poly_t f, const acb_t c, slong prec)
+
+    Sets *g* to the Taylor shift `f(x+c)`, computed respectively using
+    an optimized form of Horner's rule, divide-and-conquer, and an automatic
+    choice between the two algorithms.
+
+    The underscore methods act in-place on *g* = *f* which has length *n*.
+
 .. function:: void _acb_poly_compose_horner(acb_ptr res, acb_srcptr poly1, slong len1, acb_srcptr poly2, slong len2, slong prec)
 
 .. function:: void acb_poly_compose_horner(acb_poly_t res, const acb_poly_t poly1, const acb_poly_t poly2, slong prec)
@@ -361,6 +379,10 @@ Composition
     Sets *res* to the composition `h(x) = f(g(x))` where `f` is given by
     *poly1* and `g` is given by *poly2*, respectively using Horner's rule,
     divide-and-conquer, and an automatic choice between the two algorithms.
+
+    The default algorithm also handles special-form input `g = ax^n + c`
+    efficiently by performing a Taylor shift followed by a rescaling.
+
     The underscore methods do not support aliasing of the output
     with either input polynomial.
 
@@ -380,11 +402,13 @@ Composition
     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.
+
+    The default algorithm also handles special-form input `g = ax^n` efficiently.
+
     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.
 
-
 .. function:: void _acb_poly_revert_series_lagrange(acb_ptr h, acb_srcptr f, slong flen, slong n, slong prec)
 
 .. function:: void acb_poly_revert_series_lagrange(acb_poly_t h, const acb_poly_t f, slong n, slong prec)

doc/source/arb_poly.rst

@@ -345,6 +345,24 @@ Arithmetic
 Composition
 -------------------------------------------------------------------------------
 
+.. function:: void _arb_poly_taylor_shift_horner(arb_ptr g, const arb_t c, slong n, slong prec)
+
+.. function:: void arb_poly_taylor_shift_horner(arb_poly_t g, const arb_poly_t f, const arb_t c, slong prec)
+
+.. function:: void _arb_poly_taylor_shift_divconquer(arb_ptr g, const arb_t c, slong n, slong prec)
+
+.. function:: void arb_poly_taylor_shift_divconquer(arb_poly_t g, const arb_poly_t f, const arb_t c, slong prec)
+
+.. function:: void _arb_poly_taylor_shift(arb_ptr g, const arb_t c, slong n, slong prec)
+
+.. function:: void arb_poly_taylor_shift(arb_poly_t g, const arb_poly_t f, const arb_t c, slong prec)
+
+    Sets *g* to the Taylor shift `f(x+c)`, computed respectively using
+    an optimized form of Horner's rule, divide-and-conquer, and an automatic
+    choice between the two algorithms.
+
+    The underscore methods act in-place on *g* = *f* which has length *n*.
+
 .. function:: void _arb_poly_compose_horner(arb_ptr res, arb_srcptr poly1, slong len1, arb_srcptr poly2, slong len2, slong prec)
 
 .. function:: void arb_poly_compose_horner(arb_poly_t res, const arb_poly_t poly1, const arb_poly_t poly2, slong prec)
@@ -360,6 +378,10 @@ Composition
     Sets *res* to the composition `h(x) = f(g(x))` where `f` is given by
     *poly1* and `g` is given by *poly2*, respectively using Horner's rule,
     divide-and-conquer, and an automatic choice between the two algorithms.
+
+    The default algorithm also handles special-form input `g = ax^n + c`
+    efficiently by performing a Taylor shift followed by a rescaling.
+
     The underscore methods do not support aliasing of the output
     with either input polynomial.
 
@@ -379,11 +401,13 @@ Composition
     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.
+
+    The default algorithm also handles special-form input `g = ax^n` efficiently.
+
     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.
 
-
 .. function:: void _arb_poly_revert_series_lagrange(arb_ptr h, arb_srcptr f, slong flen, slong n, slong prec)
 
 .. function:: void arb_poly_revert_series_lagrange(arb_poly_t h, const arb_poly_t f, slong n, slong prec)