diff --git a/arb_poly.h b/arb_poly.h
index 1fa2940b2..7569b1fb6 100644
--- a/arb_poly.h
+++ b/arb_poly.h
@@ -454,6 +454,10 @@ void _arb_poly_evaluate_vec_fast(arb_ptr ys, arb_srcptr poly, slong plen,
void arb_poly_evaluate_vec_fast(arb_ptr ys,
const arb_poly_t poly, arb_srcptr xs, slong n, slong prec);
+void _arb_poly_to_taylor_model(arb_ptr g, mag_t radius, arb_srcptr f, slong len, const arb_t x, slong glen,
+ slong prec);
+
+
void _arb_poly_interpolate_newton(arb_ptr poly, arb_srcptr xs,
arb_srcptr ys, slong n, slong prec);
diff --git a/arb_poly/test/t-to_taylor_model.c b/arb_poly/test/t-to_taylor_model.c
new file mode 100644
index 000000000..da172b62b
--- /dev/null
+++ b/arb_poly/test/t-to_taylor_model.c
@@ -0,0 +1,141 @@
+/*
+ Copyright (C) 2022 Erik Postma
+
+ This file is part of Arb.
+
+ Arb is free software: you can redistribute it and/or modify it under
+ the terms of the GNU Lesser General Public License (LGPL) as published
+ by the Free Software Foundation; either version 2.1 of the License, or
+ (at your option) any later version. See .
+*/
+
+#include "arb_poly.h"
+
+int main()
+{
+ slong iter;
+ flint_rand_t state;
+
+ flint_printf("to_taylor_model....");
+ fflush(stdout);
+
+ flint_randinit(state);
+
+ for (iter = 0; iter < 50 * arb_test_multiplier(); ++iter)
+ {
+ slong prec, glen, n_iterations, y_prec, i;
+ arb_poly_t fp;
+ arb_t x;
+ arb_ptr g;
+ mag_t radius;
+
+ prec = 2 + n_randint(state, 500);
+
+ arb_poly_init(fp);
+ arb_init(x);
+ mag_init(radius);
+
+ arb_poly_randtest(fp, state, 1 + n_randint(state, 10), 1 + n_randint(state, 500), 10);
+ arb_randtest(x, state, 1 + n_randint(state, 100), 10);
+
+ glen = 1 + n_randint(state, 10); /* We can have glen >= len. That's not very useful, but it
+ * should work. */
+ g = _arb_vec_init(glen);
+
+ _arb_poly_to_taylor_model(g, radius, fp->coeffs, fp->length, x, glen, prec);
+
+ if(arb_is_exact(x))
+ {
+ n_iterations = 1;
+ y_prec = prec + 30;
+ }
+ else
+ {
+ n_iterations = 25;
+ y_prec = prec;
+ }
+
+ for(i = 0; i < n_iterations; ++i)
+ {
+ /* Generate a random point y = arb_midref(x) + yy inside the interval represented by
+ * x. Verify that f(y) is contained in g(y) +/- radius. */
+ arb_t y, yy, fy, gy;
+
+ arb_init(y);
+ arb_init(yy);
+ arb_init(fy);
+ arb_init(gy);
+
+ if(arb_is_exact(x))
+ {
+ arb_zero(yy);
+ arb_set(y, x);
+ }
+ else
+ {
+ arf_t rad_x_as_arf_t;
+ arf_init(rad_x_as_arf_t);
+
+ arb_urandom(yy, state, y_prec);
+
+ arf_set_mag(rad_x_as_arf_t, arb_radref(x));
+ arb_mul_arf(yy, yy, rad_x_as_arf_t, y_prec);
+ arb_add_arf(y, yy, arb_midref(x), y_prec);
+
+ arf_clear(rad_x_as_arf_t);
+ }
+
+ arb_poly_evaluate(fy, fp, y, y_prec + 30);
+
+ _arb_poly_evaluate(gy, g, glen, yy, y_prec);
+ flint_printf("gy before = "); arb_printd(gy, 15); flint_printf("\n");
+ arb_add_error_mag(gy, radius);
+ flint_printf("gy after = "); arb_printd(gy, 15); flint_printf("\n");
+
+ if(! (glen >= fp->length ? arb_overlaps(gy, fy) : arb_contains(gy, fy)))
+ {
+ flint_printf("FAIL\n\n");
+
+ flint_printf("prec = %d\n", prec);
+ flint_printf("len = %d\n", fp->length);
+ flint_printf("glen = %d\n\n", glen);
+
+ flint_printf("f = "); arb_poly_printd(fp, 15); flint_printf("\n");
+ flint_printf("x = "); arb_printd(x, 15); flint_printf("\n");
+ flint_printf("g = [");
+ for(i = 0; i < glen; ++i)
+ {
+ flint_printf("(");
+ arb_printd(g + i, 15);
+ if(i < glen-1)
+ flint_printf("), ");
+ else
+ flint_printf(")");
+ }
+ flint_printf("]\n");
+ flint_printf("radius = "); mag_printd(radius, 15); flint_printf("\n\n");
+
+ flint_printf("yy = "); arb_printd(yy, 15); flint_printf("\n");
+ flint_printf("y = "); arb_printd(y, 15); flint_printf("\n");
+ flint_printf("fy = "); arb_printd(fy, 50); flint_printf("\n");
+ flint_printf("gy = "); arb_printd(gy, 50); flint_printf("\n\n");
+
+ flint_abort();
+ }
+
+ arb_clear(gy);
+ arb_clear(fy);
+ arb_clear(y);
+ arb_clear(yy);
+ }
+
+ mag_clear(radius);
+ arb_clear(x);
+ arb_poly_clear(fp);
+ }
+
+ flint_randclear(state);
+ flint_cleanup();
+ flint_printf("PASS\n");
+ return EXIT_SUCCESS;
+}
diff --git a/arb_poly/to_taylor_model.c b/arb_poly/to_taylor_model.c
new file mode 100644
index 000000000..c668f1756
--- /dev/null
+++ b/arb_poly/to_taylor_model.c
@@ -0,0 +1,131 @@
+/*
+ Copyright (C) 2022 Erik Postma
+
+ This file is part of Arb.
+
+ Arb is free software: you can redistribute it and/or modify it under
+ the terms of the GNU Lesser General Public License (LGPL) as published
+ by the Free Software Foundation; either version 2.1 of the License, or
+ (at your option) any later version. See .
+*/
+
+#include "arb_poly.h"
+
+void _arb_poly_to_taylor_model(arb_ptr g, mag_t radius, arb_srcptr f, slong len, const arb_t x, slong glen,
+ slong prec)
+{
+ /*
+ * Consider a value x0 in x. Write x0 = xc + xr, with xc = arb_midref(x). Express f as an exact
+ * polynomial g in xr = x0 - xc, with a bound, radius, to enclose all values that f can assume.
+ *
+ * Define the stages of constructing the Horner form of f, expressed in an abstract variable y,
+ * as:
+ *
+ * HP[len](y) = 0,
+ * HP[d](y) = f[d] + y * HP[d+1](y), for 0 <= d < len.
+ *
+ * In addition to making sure that each g's radius is 0, we maintain this invariant in the loop
+ * below:
+ *
+ * abs(HP[i](x0) - sum(g[d] * xr^d, d = 0 .. glen-1)) <= radius.
+ *
+ * In order to establish this from the invariant for i+1, we use the definition of HP and write
+ *
+ * abs((f[i] + x0 * HP[i+1](x0)) - (f[i] + x0 * sum(g[d] * xr^d, d = 0 .. glen-1))) <=
+ * abs(f[i]) * radius.
+ *
+ * We first assume that f[i] is exact. The first two terms inside the absolute value are
+ * HP[i](x0); minus the rest can be expanded to
+ *
+ * f[i] + (xc + xr) * sum(g[d] * xr^d, d = 0 .. glen-1) =
+ * f[i] + sum(xc * g[d] * xr^d, d = 0 .. glen-1)
+ * + sum(g[d] * xr^(d+1), d = 0 .. glen-1) =
+ * f[i] + xc * g[0]
+ * + sum((g[d-1] + xc * g[d]) * xr^d, d = 1 .. glen-1)
+ * + g[glen-1] * xr^glen.
+ *
+ * So we need to:
+ *
+ * - multiply radius by abs(x);
+ * - set g[0] to f[i] + xc * g[0]
+ * - set g[d] to g[d-1] + xc * g[d] for d = 1 .. glen-1
+ * - account for g[glen-1] * xr^glen by bounding it appropriately:
+ * * if glen is odd, add abs(g[glen-1]) * arb_radref(x)^glen to radius;
+ * * if glen is even, add abs(g[glen-1]) * (arb_radref(x)^glen)/2 to radius and add g[glen-1]
+ * * (arb_radref(x)^glen)/2 to g[0].
+ *
+ * Finally, if f[i] has a nonzero radius, then we can just add it to radius.
+ */
+ arf_t g_glen_m_1; /* Stores a copy of g[glen - 1]. */
+ mag_t abs_x; /* Precompute abs(x). */
+ mag_t u; /* Scratch value. */
+ mag_t rad_x_glen_or_half; /* Precompute arb_radref(x)^glen, divided by 2 if glen is even. */
+ arf_t rad_x_glen_half_arf; /* If glen is even, same, but as an arf_t. This gets initialized
+ * *only* if glen is even. */
+ slong i, d;
+
+ arf_init(g_glen_m_1);
+
+ mag_init(abs_x);
+ arb_get_mag(abs_x, x);
+
+ mag_init(u);
+
+ mag_init(rad_x_glen_or_half);
+ mag_pow_ui(rad_x_glen_or_half, arb_radref(x), glen);
+ if(glen % 2 == 0)
+ {
+ mag_div_ui(rad_x_glen_or_half, rad_x_glen_or_half, 2);
+ arf_init(rad_x_glen_half_arf);
+ arf_set_mag(rad_x_glen_half_arf, rad_x_glen_or_half);
+ }
+
+ _arb_vec_zero(g, glen);
+ mag_zero(radius);
+
+ for(i = len - 1; i >= 0; --i)
+ {
+ /* Set radius to radius * abs(x) + abs(g[glen-1]) * arb_radref(x)^glen (maybe divided by 2)
+ * + arb_radref(f[i]). */
+ mag_mul(radius, radius, abs_x);
+ arf_get_mag(u, g_glen_m_1);
+ mag_addmul(radius, u, rad_x_glen_or_half);
+ mag_add(radius, radius, arb_radref(f + i));
+
+ /* We need to overwrite g[glen-1] before we can use it to update g[0]. */
+ arf_set(g_glen_m_1, arb_midref(g + (glen-1)));
+
+ for(d = glen - 1; d >= 1; --d)
+ {
+ /* Set g[d] to g[d-1] + xc * g[d]. */
+ arf_fma(arb_midref(g + d), arb_midref(x), arb_midref(g + d), arb_midref(g + (d-1)),
+ prec, ARF_RND_NEAR);
+ }
+
+ /* Set g[0] to f[i] + xc * g[0] (maybe plus g[glen-1] * arb_radref(x)^glen / 2). */
+ arf_fma(arb_midref(g + 0), arb_midref(x), arb_midref(g + 0), arb_midref(f + i),
+ prec, ARF_RND_NEAR);
+
+ if(glen % 2 == 0)
+ {
+ arf_addmul(arb_midref(g + 0), g_glen_m_1, rad_x_glen_half_arf, prec, ARF_RND_NEAR);
+ }
+
+ /*
+ for(d = 0; d < glen; ++d)
+ {
+ flint_printf("g[%d] = ", d); arb_printd(g + d, 15); flint_printf("\n");
+ }
+ flint_printf("radius = "); mag_printd(radius, 15); flint_printf("\n\n");
+ */
+ }
+
+ if(glen % 2 == 0)
+ {
+ arf_clear(rad_x_glen_half_arf);
+ }
+ mag_clear(rad_x_glen_or_half);
+ mag_clear(u);
+ mag_clear(abs_x);
+ arf_clear(g_glen_m_1);
+}
diff --git a/doc/source/arb_poly.rst b/doc/source/arb_poly.rst
index 8d74b87c3..bb827b87a 100644
--- a/doc/source/arb_poly.rst
+++ b/doc/source/arb_poly.rst
@@ -550,6 +550,11 @@ Evaluation
Sets `y = f(x), z = f'(x)`, evaluated respectively using Horner's rule,
rectangular splitting, and an automatic algorithm choice.
+.. function:: void _arb_poly_to_taylor_model(arb_ptr g, mag_t radius, arb_srcptr f, slong len, const
+ arb_t x, slong glen, slong prec)
+
+ Sets the `glen` entries of `g` to be exact coefficients of a polynomial, such that for values
+ `x_0` in `x`, we have `abs(g(x_0) - f(x_0)) <= radius`.
Product trees
-------------------------------------------------------------------------------