From 2f942f1d226d54c24723e5f0ac6fdf55145acc2e Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Thu, 3 Oct 2013 13:03:55 +0200 Subject: [PATCH] first stab at some integration code; fix cauchy_bound --- doc/source/fmpcb_calc.rst | 60 +++++++- fmpcb_calc.h | 9 ++ fmpcb_calc/cauchy_bound.c | 1 - fmpcb_calc/integrate_taylor.c | 219 +++++++++++++++++++++++++++ fmpcb_calc/test/t-cauchy_bound.c | 8 +- fmpcb_calc/test/t-integrate_taylor.c | 109 +++++++++++++ 6 files changed, 399 insertions(+), 7 deletions(-) create mode 100644 fmpcb_calc/integrate_taylor.c create mode 100644 fmpcb_calc/test/t-integrate_taylor.c diff --git a/doc/source/fmpcb_calc.rst b/doc/source/fmpcb_calc.rst index 05b97f26..e1e87c24 100644 --- a/doc/source/fmpcb_calc.rst +++ b/doc/source/fmpcb_calc.rst @@ -34,13 +34,14 @@ Types, macros and constants Bounds ------------------------------------------------------------------------------- -.. function :: void fmpcb_calc_cauchy_bound(fmprb_t bound, fmpcb_calc_func_t func, void * param, const fmpcb_t x, const fmprb_t radius, long maxdepth, long prec) +.. function:: void fmpcb_calc_cauchy_bound(fmprb_t bound, fmpcb_calc_func_t func, void * param, const fmpcb_t x, const fmprb_t radius, long maxdepth, long prec) Sets *bound* to a ball containing the value of the integral .. math :: - b = \frac{1}{2 \pi r} \oint_{|t-x| = r} |f(t)| dt + C(x,r) = \frac{1}{2 \pi r} \oint_{|z-x| = r} |f(z)| dz + = \int_0^1 |f(x+re^{2\pi i t})| dt where *f* is specified by (*func*, *param*) and *r* is given by *radius*. The integral is computed using a simple step sum. @@ -51,3 +52,58 @@ Bounds repeatedly subdivides the whole integration range instead of performing adaptive subdivisions. +Integration +------------------------------------------------------------------------------- + +.. function:: int fmpcb_calc_integrate_taylor(fmpcb_t res, fmpcb_calc_func_t func, void * param, const fmpcb_t a, const fmpcb_t b, const fmpr_t inner_radius, const fmpr_t outer_radius, long accuracy_goal, long prec) + + Computes the integral + + .. math :: + + I = \int_a^b f(t) dt + + where *f* is specified by (*func*, *param*), following a straight-line + path between the complex numbers *a* and *b* which both must be finite. + + The integral is approximated by piecewise centered Taylor polynomials. + Rigorous truncation error bounds are calculated using the Cauchy integral + formula. More precisely, if the Taylor series of *f* centered at the point + *m* is `f(m+x) = \sum_{n=0}^{\infty} a_n x^n`, then + + .. math :: + + \int f(m+x) = \left( \sum_{n=0}^{N-1} a_n \frac{x^{n+1}}{n+1} \right) + + \left( \sum_{n=N}^{\infty} a_n \frac{x^{n+1}}{n+1} \right). + + For sufficiently small *x*, the second series converges and its + absolute value is bounded by + + .. math :: + + \sum_{n=N}^{\infty} \frac{C(m,R)}{R^n} \frac{|x|^{n+1}}{N+1} + = \frac{C(m,R) R x}{(R-x)(N+1)} \left( \frac{x}{R} \right)^N. + + It is required that any singularities of *f* are + isolated from the path of integration by a distance strictly + greater than the positive value *outer_radius* (which is the integration + radius used for the Cauchy bound). Taylor series step lengths are + chosen so as not to + exceed *inner_radius*, which must be strictly smaller than *outer_radius* + for convergence. A smaller *inner_radius* gives more rapid convergence + of each Taylor series but means that more series might have to be used. + A reasonable choice might be to set *inner_radius* to half the value of + *outer_radius*, giving roughly one accurate bit per term. + + The truncation point of each Taylor series is chosen so that the absolute + truncation error is roughly `2^{-p}` where *p* is given by *accuracy_goal* + (in the future, this might change to a relative accuracy). + Arithmetic operations and function + evaluations are performed at a precision of *prec* bits. Note that due + to accumulation of numerical errors, both values may have to be set + higher (and the endpoints may have to be computed more accurately) + to achieve a desired accuracy. + + This function chooses the evaluation points uniformly rather + than implementing adaptive subdivision. + diff --git a/fmpcb_calc.h b/fmpcb_calc.h index 4c47dd8d..6ca400db 100644 --- a/fmpcb_calc.h +++ b/fmpcb_calc.h @@ -44,6 +44,15 @@ void fmpcb_calc_cauchy_bound(fmprb_t bound, fmpcb_calc_func_t func, void * param, const fmpcb_t x, const fmprb_t radius, long maxdepth, long prec); +/* Integration */ + +int fmpcb_calc_integrate_taylor(fmpcb_t res, + fmpcb_calc_func_t func, void * param, + const fmpcb_t a, const fmpcb_t b, + const fmpr_t inner_radius, + const fmpr_t outer_radius, + long accuracy_goal, long prec); + #ifdef __cplusplus } #endif diff --git a/fmpcb_calc/cauchy_bound.c b/fmpcb_calc/cauchy_bound.c index 78b29857..a555f740 100644 --- a/fmpcb_calc/cauchy_bound.c +++ b/fmpcb_calc/cauchy_bound.c @@ -69,7 +69,6 @@ fmpcb_calc_cauchy_bound(fmprb_t bound, fmpcb_calc_func_t func, void * param, } fmprb_div_ui(b, b, n, prec); - fmprb_div(b, b, radius, prec); if (fmprb_is_exact(b) || fmpr_cmp(fmprb_radref(b), fmprb_midref(b)) < 0) break; diff --git a/fmpcb_calc/integrate_taylor.c b/fmpcb_calc/integrate_taylor.c new file mode 100644 index 00000000..9e80798c --- /dev/null +++ b/fmpcb_calc/integrate_taylor.c @@ -0,0 +1,219 @@ +/*============================================================================= + + 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) 2013 Fredrik Johansson + +******************************************************************************/ + +#include "fmpcb_calc.h" +#include "math.h" + +int +fmpcb_calc_integrate_taylor(fmpcb_t res, + fmpcb_calc_func_t func, void * param, + const fmpcb_t a, const fmpcb_t b, + const fmpr_t inner_radius, + const fmpr_t outer_radius, + long accuracy_goal, long prec) +{ + long num_steps, step, N, bp; + int result; + + fmpcb_t delta, m, x, y1, y2, sum; + fmpcb_ptr taylor_poly; + fmpr_t err; + + fmpcb_init(delta); + fmpcb_init(m); + fmpcb_init(x); + fmpcb_init(y1); + fmpcb_init(y2); + fmpcb_init(sum); + fmpr_init(err); + + fmpcb_sub(delta, b, a, prec); + + /* precision used for bounds calculations */ + bp = FMPRB_RAD_PREC; + + /* compute the number of steps */ + { + fmpr_t t; + fmpr_init(t); + fmpcb_get_abs_ubound_fmpr(t, delta, bp); + fmpr_div(t, t, inner_radius, bp, FMPR_RND_UP); + fmpr_mul_2exp_si(t, t, -1); + num_steps = (long) (fmpr_get_d(t, FMPR_RND_UP) + 1.0); + /* make sure it's not something absurd */ + num_steps = FLINT_MIN(num_steps, 10 * prec); + num_steps = FLINT_MAX(num_steps, 1); + fmpr_clear(t); + } + + result = FMPRB_CALC_SUCCESS; + + fmpcb_zero(sum); + + for (step = 0; step < num_steps; step++) + { + /* midpoint of subinterval */ + fmpcb_mul_ui(m, delta, 2 * step + 1, prec); + fmpcb_div_ui(m, m, 2 * num_steps, prec); + fmpcb_add(m, m, a, prec); + + if (fmprb_calc_verbose) + { + printf("integration point %ld/%ld: ", 2 * step + 1, 2 * num_steps); + fmpcb_printd(m, 15); printf("\n"); + } + + /* evaluate at +/- x */ + /* TODO: exactify m, and include error in x? */ + fmpcb_div_ui(x, delta, 2 * num_steps, prec); + + /* compute bounds and number of terms to use */ + { + fmprb_t cbound, xbound, rbound; + fmpr_t C, D, R, X, T; + double DD, TT, NN; + + fmprb_init(cbound); + fmprb_init(xbound); + fmprb_init(rbound); + fmpr_init(C); + fmpr_init(D); + fmpr_init(R); + fmpr_init(X); + fmpr_init(T); + + /* R is the outer radius */ + fmpr_set(R, outer_radius); + + /* X = upper bound for |x| */ + fmpcb_get_abs_ubound_fmpr(X, x, bp); + fmprb_set_fmpr(xbound, X); + + /* Compute C(m,R). Important subtlety: due to rounding when + computing m, we will in general be farther than R away from + the integration path. But since fmpcb_calc_cauchy_bound + actually integrates over the area traced by a complex + interval, it will catch any extra singularities (giving + an infinite bound). */ + fmprb_set_fmpr(rbound, outer_radius); + fmpcb_calc_cauchy_bound(cbound, func, param, m, rbound, 8, bp); + fmpr_add(C, fmprb_midref(cbound), fmprb_radref(cbound), bp, FMPR_RND_UP); + + /* Sanity check: we need C < inf and R > X */ + if (fmpr_is_finite(C) && fmpr_cmp(R, X) > 0) + { + /* Compute upper bound for D = C * R * X / (R - X) */ + fmpr_mul(D, C, R, bp, FMPR_RND_UP); + fmpr_mul(D, D, X, bp, FMPR_RND_UP); + fmpr_sub(T, R, X, bp, FMPR_RND_DOWN); + fmpr_div(D, D, T, bp, FMPR_RND_UP); + + /* Compute upper bound for T = (X / R) */ + fmpr_div(T, X, R, bp, FMPR_RND_UP); + + /* Choose N */ + /* TODO: use fmpr arithmetic to avoid overflow */ + /* TODO: use relative accuracy (look at |f(m)|?) */ + DD = fmpr_get_d(D, FMPR_RND_UP); + TT = fmpr_get_d(T, FMPR_RND_UP); + NN = -(accuracy_goal * 0.69314718055994530942 + log(DD)) / log(TT); + N = NN + 0.5; + N = FLINT_MIN(N, 100 * prec); + N = FLINT_MAX(N, 1); + + /* Tail bound: D / (N + 1) * T^N */ + fmpr_pow_sloppy_ui(T, T, N, bp, FMPR_RND_UP); + fmpr_mul(D, D, T, bp, FMPR_RND_UP); + fmpr_div_ui(err, D, N + 1, bp, FMPR_RND_UP); + } + else + { + N = 1; + fmpr_pos_inf(err); + result = FMPRB_CALC_NO_CONVERGENCE; + } + + if (fmprb_calc_verbose) + { + printf("N = %ld; bound: ", N); fmpr_printd(err, 15); printf("\n"); + printf("R: "); fmpr_printd(R, 15); printf("\n"); + printf("C: "); fmpr_printd(C, 15); printf("\n"); + printf("X: "); fmpr_printd(X, 15); printf("\n"); + } + + fmprb_clear(cbound); + fmprb_clear(xbound); + fmprb_clear(rbound); + fmpr_clear(C); + fmpr_clear(D); + fmpr_clear(R); + fmpr_clear(X); + fmpr_clear(T); + } + + /* evaluate Taylor polynomial */ + taylor_poly = _fmpcb_vec_init(N + 1); + func(taylor_poly, m, param, N, prec); + _fmpcb_poly_integral(taylor_poly, taylor_poly, N + 1, prec); + _fmpcb_poly_evaluate(y2, taylor_poly, N + 1, x, prec); + fmpcb_neg(x, x); + _fmpcb_poly_evaluate(y1, taylor_poly, N + 1, x, prec); + fmpcb_neg(x, x); + + /* add truncation error */ + fmprb_add_error_fmpr(fmpcb_realref(y1), err); + fmprb_add_error_fmpr(fmpcb_imagref(y1), err); + fmprb_add_error_fmpr(fmpcb_realref(y2), err); + fmprb_add_error_fmpr(fmpcb_imagref(y2), err); + + fmpcb_add(sum, sum, y2, prec); + fmpcb_sub(sum, sum, y1, prec); + + if (fmprb_calc_verbose) + { + printf("values: "); + fmpcb_printd(y1, 15); printf(" "); + fmpcb_printd(y2, 15); printf("\n"); + } + + _fmpcb_vec_clear(taylor_poly, N + 1); + + if (result == FMPRB_CALC_NO_CONVERGENCE) + break; + } + + fmpcb_set(res, sum); + + fmpcb_clear(delta); + fmpcb_clear(m); + fmpcb_clear(x); + fmpcb_clear(y1); + fmpcb_clear(y2); + fmpcb_clear(sum); + fmpr_clear(err); + + return result; +} + diff --git a/fmpcb_calc/test/t-cauchy_bound.c b/fmpcb_calc/test/t-cauchy_bound.c index e313c4f2..9c919c5d 100644 --- a/fmpcb_calc/test/t-cauchy_bound.c +++ b/fmpcb_calc/test/t-cauchy_bound.c @@ -40,10 +40,10 @@ sin_x(fmpcb_ptr out, const fmpcb_t inp, void * params, long order, long prec) } static const double answers[10] = { - 1.04570093561423, 1.01793337548660, 1.60902133468552, - 2.80867584626455, 5.42667597855619, 11.2017176650732, - 24.0806215966243, 53.4370258520436, 121.507732609920, - 281.570144392142 + 1.04570093561423094, 2.0358667496686487, 4.82706400405656566, + 11.2347033850581819, 27.1331828778516522, 67.210305990439042, + 168.564351176369878, 427.496180295369909, 1093.56959348921777, + 2815.70144392142227 }; int main() diff --git a/fmpcb_calc/test/t-integrate_taylor.c b/fmpcb_calc/test/t-integrate_taylor.c new file mode 100644 index 00000000..87b7df2e --- /dev/null +++ b/fmpcb_calc/test/t-integrate_taylor.c @@ -0,0 +1,109 @@ +/*============================================================================= + + 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) 2013 Fredrik Johansson + +******************************************************************************/ + +#include "fmpcb_calc.h" + +/* sin(x) */ +int +sin_x(fmpcb_ptr out, const fmpcb_t inp, void * params, long order, long prec) +{ + int xlen = FLINT_MIN(2, order); + + fmpcb_set(out, inp); + if (xlen > 1) + fmpcb_one(out + 1); + + _fmpcb_poly_sin_series(out, out, xlen, order, prec); + return 0; +} + +int main() +{ + long iter; + flint_rand_t state; + + printf("integrate_taylor...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 150; iter++) + { + fmpcb_t ans, res, a, b; + fmpr_t inr, outr; + double t; + long goal, prec; + + fmpcb_init(ans); + fmpcb_init(res); + fmpcb_init(a); + fmpcb_init(b); + fmpr_init(inr); + fmpr_init(outr); + + goal = 2 + n_randint(state, 300); + prec = 2 + n_randint(state, 300); + + fmpcb_randtest(a, state, 1 + n_randint(state, 200), 2); + fmpcb_randtest(b, state, 1 + n_randint(state, 200), 2); + + fmpcb_cos(ans, a, prec); + fmpcb_cos(res, b, prec); + fmpcb_sub(ans, ans, res, prec); + + t = (1 + n_randint(state, 20)) / 10.0; + fmpr_set_d(inr, t); + fmpr_set_d(outr, t + (1 + n_randint(state, 20)) / 5.0); + + fmpcb_calc_integrate_taylor(res, sin_x, NULL, + a, b, inr, outr, goal, prec); + + if (!fmpcb_overlaps(res, ans)) + { + printf("FAIL! (iter = %ld)\n", iter); + printf("prec = %ld, goal = %ld\n", prec, goal); + printf("inr = "); fmpr_printd(inr, 15); printf("\n"); + printf("outr = "); fmpr_printd(outr, 15); printf("\n"); + printf("a = "); fmpcb_printd(a, 15); printf("\n"); + printf("b = "); fmpcb_printd(b, 15); printf("\n"); + printf("res = "); fmpcb_printd(res, 15); printf("\n\n"); + printf("ans = "); fmpcb_printd(ans, 15); printf("\n\n"); + abort(); + } + + fmpcb_clear(ans); + fmpcb_clear(res); + fmpcb_clear(a); + fmpcb_clear(b); + fmpr_clear(inr); + fmpr_clear(outr); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} +