Add _arb_poly_to_taylor_model, for creating a Taylor model of a polynomial

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);
 

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 <http://www.gnu.org/licenses/>.
+*/
+
+#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;
+}

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 <http://www.gnu.org/licenses/>.
+*/
+
+#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);
+}

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
 -------------------------------------------------------------------------------