first stab at some integration code; fix cauchy_bound

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

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

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;

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;
+}
+

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,
+  1.04570093561423094, 2.0358667496686487, 4.82706400405656566,
-  2.80867584626455, 5.42667597855619, 11.2017176650732,
+  11.2347033850581819, 27.1331828778516522, 67.210305990439042,
-  24.0806215966243, 53.4370258520436, 121.507732609920,
+  168.564351176369878, 427.496180295369909, 1093.56959348921777,
-  281.570144392142
+  2815.70144392142227
 };
 
 int main()

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;
+}
+