From f0b041228a9b9e2f8421fe8f846adfe7674638ce Mon Sep 17 00:00:00 2001 From: postmath Date: Thu, 7 Apr 2022 13:08:50 -0400 Subject: [PATCH] Add _arb_poly_to_taylor_model, for creating a Taylor model (in the sense of Makino and Berz) of a polynomial. --- arb_poly.h | 4 + arb_poly/test/t-to_taylor_model.c | 141 ++++++++++++++++++++++++++++++ arb_poly/to_taylor_model.c | 131 +++++++++++++++++++++++++++ doc/source/arb_poly.rst | 5 ++ 4 files changed, 281 insertions(+) create mode 100644 arb_poly/test/t-to_taylor_model.c create mode 100644 arb_poly/to_taylor_model.c 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 -------------------------------------------------------------------------------