From c9cca291d57e59ae67e2eb5b35a1cedf403d8f4b Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Wed, 22 Nov 2017 00:24:20 +0100 Subject: [PATCH] minor improvements to integration code; change interface; more examples --- acb_calc.h | 30 ++- acb_calc/integrate.c | 71 ++++-- acb_calc/integrate_gl_auto_deg.c | 5 +- acb_calc/integrate_opt_init.c | 23 ++ doc/source/acb_calc.rst | 188 +++++++++------ doc/source/credits.rst | 1 + examples/integrals.c | 396 +++++++++++++++++++++++++------ 7 files changed, 539 insertions(+), 175 deletions(-) create mode 100644 acb_calc/integrate_opt_init.c diff --git a/acb_calc.h b/acb_calc.h index 3f6ba8bf..1f9c0fc8 100644 --- a/acb_calc.h +++ b/acb_calc.h @@ -21,22 +21,15 @@ extern "C" { #endif -#define ACB_CALC_VERBOSE 1 -#define ACB_CALC_VERY_VERBOSE 3 - -#define ACB_CALC_INTEGRATE_HEAP 8 - typedef int (*acb_calc_func_t)(acb_ptr out, const acb_t inp, void * param, slong order, slong prec); -/* Bounds */ +/* Integration (old) */ void acb_calc_cauchy_bound(arb_t bound, acb_calc_func_t func, void * param, const acb_t x, const arb_t radius, slong maxdepth, slong prec); -/* Integration */ - int acb_calc_integrate_taylor(acb_t res, acb_calc_func_t func, void * param, const acb_t a, const acb_t b, @@ -44,19 +37,34 @@ int acb_calc_integrate_taylor(acb_t res, const arf_t outer_radius, slong accuracy_goal, slong prec); +/* Integration */ + +typedef struct +{ + slong deg_limit; + slong eval_limit; + slong depth_limit; + int use_heap; + int verbose; +} +acb_calc_integrate_opt_struct; + +typedef acb_calc_integrate_opt_struct acb_calc_integrate_opt_t[1]; + +void acb_calc_integrate_opt_init(acb_calc_integrate_opt_t options); + int acb_calc_integrate(acb_t res, acb_calc_func_t f, void * param, const acb_t a, const acb_t b, slong goal, const mag_t tol, - slong deg_limit, slong eval_limit, slong depth_limit, - int flags, + const acb_calc_integrate_opt_t options, slong prec); int acb_calc_integrate_gl_auto_deg(acb_t res, slong * eval_count, acb_calc_func_t f, void * param, const acb_t a, const acb_t b, const mag_t tol, - slong deg_limit, int flags, slong prec); + slong deg_limit, int verbose, slong prec); #ifdef __cplusplus } diff --git a/acb_calc/integrate.c b/acb_calc/integrate.c index 14237267..e3f60bed 100644 --- a/acb_calc/integrate.c +++ b/acb_calc/integrate.c @@ -108,17 +108,25 @@ int acb_calc_integrate(acb_t res, acb_calc_func_t f, void * param, const acb_t a, const acb_t b, slong goal, const mag_t tol, - slong deg_limit, slong eval_limit, slong depth_limit, - int flags, + const acb_calc_integrate_opt_t options, slong prec) { acb_ptr as, bs, vs; mag_ptr ms; acb_t s, t, u; mag_t tmpm, tmpn, new_tol; + slong depth_limit, eval_limit, deg_limit; slong depth, depth_max, eval, feval, top; slong leaf_interval_count; - int stopping, real_error, use_heap, status, gl_status; + slong alloc; + int stopping, real_error, use_heap, status, gl_status, verbose; + + if (options == NULL) + { + acb_calc_integrate_opt_t opt; + acb_calc_integrate_opt_init(opt); + return acb_calc_integrate(res, f, param, a, b, goal, tol, opt, prec); + } status = ARB_CALC_SUCCESS; @@ -129,25 +137,29 @@ acb_calc_integrate(acb_t res, acb_calc_func_t f, void * param, mag_init(tmpn); mag_init(new_tol); + depth_limit = options->depth_limit; if (depth_limit <= 0) depth_limit = 2 * prec; depth_limit = FLINT_MAX(depth_limit, 1); + eval_limit = options->eval_limit; if (eval_limit <= 0) eval_limit = 1000 * prec + prec * prec; eval_limit = FLINT_MAX(eval_limit, 1); goal = FLINT_MAX(goal, 0); + deg_limit = options->deg_limit; if (deg_limit <= 0) - deg_limit = 0.5 * goal + 10; + deg_limit = 0.5 * FLINT_MIN(goal, prec) + 10; - use_heap = (flags & ACB_CALC_INTEGRATE_HEAP); + verbose = options->verbose; + use_heap = options->use_heap; - /* todo: allocate dynamically */ - as = _acb_vec_init(depth_limit); - bs = _acb_vec_init(depth_limit); - vs = _acb_vec_init(depth_limit); - ms = _mag_vec_init(depth_limit); + alloc = 4; + as = _acb_vec_init(alloc); + bs = _acb_vec_init(alloc); + vs = _acb_vec_init(alloc); + ms = _mag_vec_init(alloc); /* Compute initial crude estimate for the whole interval. */ acb_set(as, a); @@ -171,7 +183,7 @@ acb_calc_integrate(acb_t res, acb_calc_func_t f, void * param, { if (stopping == 0 && eval >= eval_limit - 1) { - if (flags & ACB_CALC_VERBOSE) + if (verbose > 0) flint_printf("stopping at eval_limit %wd\n", eval_limit); status = ARB_CALC_NO_CONVERGENCE; stopping = 1; @@ -205,16 +217,16 @@ acb_calc_integrate(acb_t res, acb_calc_func_t f, void * param, /* Attempt using Gauss-Legendre rule. */ if (acb_is_finite(vs + top)) { - /* We know that the result is real. */ - real_error = acb_is_finite(vs + top) && acb_is_real(vs + top); - gl_status = acb_calc_integrate_gl_auto_deg(u, &feval, f, param, - as + top, bs + top, new_tol, deg_limit, flags, prec); + as + top, bs + top, new_tol, deg_limit, verbose > 1, prec); eval += feval; /* We are done with this subinterval. */ if (gl_status == ARB_CALC_SUCCESS) { + /* We know that the result is real. */ + real_error = acb_is_finite(vs + top) && acb_is_real(vs + top); + if (real_error) arb_zero(acb_imagref(u)); @@ -241,13 +253,30 @@ acb_calc_integrate(acb_t res, acb_calc_func_t f, void * param, if (depth >= depth_limit - 1) { - if (flags & ACB_CALC_VERBOSE) + if (verbose > 0) flint_printf("stopping at depth_limit %wd\n", depth_limit); status = ARB_CALC_NO_CONVERGENCE; stopping = 1; continue; } + if (depth >= alloc - 1) + { + slong k; + as = flint_realloc(as, 2 * alloc * sizeof(acb_struct)); + bs = flint_realloc(bs, 2 * alloc * sizeof(acb_struct)); + vs = flint_realloc(vs, 2 * alloc * sizeof(acb_struct)); + ms = flint_realloc(ms, 2 * alloc * sizeof(mag_struct)); + for (k = alloc; k < 2 * alloc; k++) + { + acb_init(as + k); + acb_init(bs + k); + acb_init(vs + k); + mag_init(ms + k); + } + alloc *= 2; + } + /* Bisection. */ /* Interval [depth] becomes [mid, b]. */ acb_set(bs + depth, bs + top); @@ -294,7 +323,7 @@ acb_calc_integrate(acb_t res, acb_calc_func_t f, void * param, depth_max = FLINT_MAX(depth, depth_max); } - if (flags & ACB_CALC_VERBOSE) + if (verbose > 0) { flint_printf("depth %wd/%wd, eval %wd/%wd, %wd leaf intervals\n", depth_max, depth_limit, eval, eval_limit, leaf_interval_count); @@ -302,10 +331,10 @@ acb_calc_integrate(acb_t res, acb_calc_func_t f, void * param, acb_set(res, s); - _acb_vec_clear(as, depth_limit); - _acb_vec_clear(bs, depth_limit); - _acb_vec_clear(vs, depth_limit); - _mag_vec_clear(ms, depth_limit); + _acb_vec_clear(as, alloc); + _acb_vec_clear(bs, alloc); + _acb_vec_clear(vs, alloc); + _mag_vec_clear(ms, alloc); acb_clear(s); acb_clear(t); acb_clear(u); diff --git a/acb_calc/integrate_gl_auto_deg.c b/acb_calc/integrate_gl_auto_deg.c index c1811909..9b79e7e0 100644 --- a/acb_calc/integrate_gl_auto_deg.c +++ b/acb_calc/integrate_gl_auto_deg.c @@ -115,7 +115,7 @@ int acb_calc_integrate_gl_auto_deg(acb_t res, slong * eval_count, acb_calc_func_t f, void * param, const acb_t a, const acb_t b, const mag_t tol, - slong deg_limit, int flags, slong prec) + slong deg_limit, int verbose, slong prec) { acb_t mid, delta, wide; mag_t tmpm; @@ -130,6 +130,7 @@ acb_calc_integrate_gl_auto_deg(acb_t res, slong * eval_count, if (deg_limit <= 0) { acb_indeterminate(res); + eval_count[0] = 0; return status; } @@ -239,7 +240,7 @@ acb_calc_integrate_gl_auto_deg(acb_t res, slong * eval_count, arb_init(x); arb_init(w); - if ((flags & ACB_CALC_VERY_VERBOSE) == ACB_CALC_VERY_VERBOSE) + if (verbose) { acb_get_mag(tmpm, delta); flint_printf(" {GL deg %ld on [", best_n); diff --git a/acb_calc/integrate_opt_init.c b/acb_calc/integrate_opt_init.c new file mode 100644 index 00000000..b8793329 --- /dev/null +++ b/acb_calc/integrate_opt_init.c @@ -0,0 +1,23 @@ +/* + Copyright (C) 2017 Fredrik Johansson + + 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 "acb_calc.h" + +void +acb_calc_integrate_opt_init(acb_calc_integrate_opt_t options) +{ + options->deg_limit = 0; + options->eval_limit = 0; + options->depth_limit = 0; + options->use_heap = 0; + options->verbose = 0; +} + diff --git a/doc/source/acb_calc.rst b/doc/source/acb_calc.rst index ffcd9838..06251c30 100644 --- a/doc/source/acb_calc.rst +++ b/doc/source/acb_calc.rst @@ -47,8 +47,9 @@ Types, macros and constants However, manual action is needed for bounded functions with branch cuts. For example, when evaluating `\sqrt{z}`, the output must be set to an non-finite value if `z` overlaps with the branch cut `[-\infty,0]`. - The easiest solution is to use holomorphy-testing versions of atomic - functions. + The easiest way to accomplish this is to use versions of basic + functions (sqrt, log, pow, etc.) that test holomorphicity of their + arguments individually. For example, here we define a piecewise holomorphic extension of the function @@ -98,7 +99,7 @@ Types, macros and constants Integration ------------------------------------------------------------------------------- -.. function:: int acb_calc_integrate(acb_t res, acb_calc_func_t func, void * param, const acb_t a, const acb_t b, slong goal, const mag_t tol, slong deg_limit, slong eval_limit, slong depth_limit, int flags, slong prec) +.. function:: int acb_calc_integrate(acb_t res, acb_calc_func_t func, void * param, const acb_t a, const acb_t b, slong rel_goal, const mag_t abs_tol, const acb_calc_integrate_opt_t options, slong prec) Computes a rigorous enclosure of the integral @@ -113,9 +114,9 @@ Integration To compute improper integrals, the user should therefore truncate the path of integration manually (or make a regularizing change of variables, if possible). - - Returns *ARB_CALC_SUCCESS* (0) if the integration converged on all - subintervals and *ARB_CALC_NO_CONVERGENCE* otherwise. + Returns *ARB_CALC_SUCCESS* if the integration converged to the + target accuracy on all subintervals, and returns + *ARB_CALC_NO_CONVERGENCE* otherwise. By default, the integrand *func* will only be called with *order* = 0 or *order* = 1; that is, derivatives are not required. @@ -137,90 +138,136 @@ Integration discontinuities while providing exponential convergence for typical piecewise holomorphic integrands. - On each subinterval, the algorithm attempts to - achieve a relative error of `2^{-\text{goal}}` or an absolute error of - *tol*, whichever is larger. The parameters *goal* and *tol* are only - guidelines; the cumulative error may be larger than both the prescribed + The following parameters control accuracy: + + - *rel_goal* - relative accuracy goal as a number of bits, i.e. + target a relative error of `\varepsilon_{rel} = 2^{-\text{rel_goal}}` + (note the sign: *rel_goal* should be nonnegative). + + - *abs_tol* - absolute accuracy goal as a :type:`mag_t` describing + the error tolerance, i.e. + target an absolute error of `\varepsilon_{abs} = \text{abs_tol}`. + + - *prec* - working precision. This is the working precision used to + evaluate the integrand and manipulate interval endpoints. + As currently implemented, the algorithm does not attempt to adjust the + working precision by itself, and adaptive + control of the working precision must be handled by the user. + + For typical usage, set *rel_goal* = *prec* and *abs_tol* = `2^{-\text{prec}}`. + It usually only makes sense to have *rel_goal* between 0 and *prec*. + + The algorithm attempts to achieve an error of + `\max(\varepsilon_{abs}, M \varepsilon_{rel})` on each subinterval, + where *M* is the magnitude of the integral. + These parameters are only guidelines; the cumulative error may be larger + than both the prescribed absolute and relative error goals, depending on the number of subdivisions, cancellation between segments of the integral, and numerical errors in the evaluation of the integrand. - The following parameters control the integration. + To compute tiny integrals with high relative accuracy, one should set + `\varepsilon_{abs} \approx M \varepsilon_{rel}` where *M* is a known + estimate of the magnitude. Setting `\varepsilon_{abs}` to 0 is also + allowed, forcing use of a relative instead of an absolute tolerance goal. + This can be handy for exponentially small or + large functions of unknown magnitude. It is recommended to avoid + setting `\varepsilon_{abs}` very small + if possible since the algorithm might need many extra + subdivisions to estimate *M* automatically; if the approximate + magnitude can be estimated by some external means (for example if + a midpoint-width or endpoint-width estimate is known to be accurate), + providing an appropriate `\varepsilon_{abs} \approx M \varepsilon_{rel}` + will be more efficient. - - *prec* - working precision. This is the working precision used to - evaluate the integrand and manipulate interval endpoints. + If the integral has very large magnitude, setting the absolute + tolerance to a corresponding large value is recommended for best + performance, but it is not necessary for convergence since the absolute + tolerance is increased automatically during the execution of the + algorithm if the partial integrals are found to have larger error. - - *goal* - relative accuracy goal as a nonnegative number of bits, i.e. - target a relative error of `2^{-\text{goal}}`. This parameter can - simply be set to *prec* (or a slightly smaller value) in most situations. + Additional options for the integration can be provided via the *options* + parameter (documented below). To use all defaults, *NULL* can be passed + for *options*. - - *tol* - absolute tolerance goal (specified as a :type:`mag_t`). - In general, *tol* should be set to about `2^{-\text{goal}}` times - an estimate for the magnitude of the integral. For typical - integrals which will have magnitude within a few orders of magnitude - of unity, `\text{tol} = 2^{-\text{goal}}` works well enough. - Of course, if the integral has very small magnitude, the absolute - tolerance must be set to a small value to get high relative accuracy. +Options for integration +............................................................................... - If the integral has very large magnitude, setting the absolute - tolerance to a corresponding large value is recommended for best - performance, but it is not necessary for convergence since the absolute - tolerance is increased automatically during the execution of the - algorithm if the partial integrals are found to have larger error. +.. type:: acb_calc_integrate_opt_struct - Setting *tol* to 0 is allowed and forces use of a relative instead of an - absolute tolerance goal, which can be handy for exponentially small or - large functions of unknown magnitude. It is recommended to avoid this - solution if possible since the algorithm might need many extra - subdivisions to determine an initial scale; if the approximate - magnitude can be estimated by some external means (for example if - a midpoint-width or endpoint-width estimate is known to be accurate), - the caller should set *tol* to `2^{-\text{goal}}` - times such an external estimate for best performance. +.. type:: acb_calc_integrate_opt_t - - *deg_limit* - maximum quadrature degree for each subinterval. - If a zero or negative value is provided, the limit is set to a default - value which currently equals `0.5 \cdot \text{goal} + 10` for - Gauss-Legendre quadrature. + This structure contains several fields, explained below. + An *acb_calc_integrate_opt_t* is defined as an array of + *acb_calc_integrate_opt_struct* + of length 1, permitting it to be passed by reference. + An *acb_calc_integrate_opt_t* must be initialized before use, which sets + all fields to 0 or *NULL*. For fields that have not been set to other + values, the integration algorithm will choose defaults automatically + (based on the precision and accuracy goals). + This structure will most likely be extended in the future to + accommodate more options. - A higher quadrature degree can be beneficial for functions that - are holomorphic on a large domain around the integration path - and yet behave irregularly, such as entire functions with a high - amount of oscillation. The drawback of increasing the degree is that - the precomputation time for quadrature nodes increases. + .. member:: slong deg_limit - - *eval_limit* - maximum number of function evaluations. - If a zero or negative value is provided, the limit is set to a default - value which currently equals `1000 \cdot \text{prec} + \text{prec}^2`. + Maximum quadrature degree for each subinterval. + If a zero or negative value is provided, the limit is set to a default + value which currently equals `0.5 \cdot \min(\text{prec}, \text{rel_goal}) + 10` for + Gauss-Legendre quadrature. + A higher quadrature degree can be beneficial for functions that + are holomorphic on a large domain around the integration path + and yet behave irregularly, such as oscillatory entire functions. + The drawback of increasing the degree is that + the precomputation time for quadrature nodes increases. - This is the main parameter used to limit the amount of work before - aborting due to possible slow convergence or non-convergence. - (This limit is only taken as a rough guideline, and the actual number of - function evaluations may be slightly higher depending on the - actual subdivisions.) - A lower limit allows aborting faster. A higher limit may be needed - for integrands with many discontinuities or many singularities - close to the integration path. + .. member:: slong eval_limit - - *depth_limit* - maximum subdivision depth. - If a zero or negative value is provided, the limit is set to the - default value `2 \cdot \text{prec}`. - This limits the number of recursive bisections around any single point. - Warning: the memory usage increases proportionally. + Maximum number of function evaluations. + If a zero or negative value is provided, the limit is set to a default + value which currently equals `1000 \cdot \text{prec} + \text{prec}^2`. + This is the main parameter used to limit the amount of work before + aborting due to possible slow convergence or non-convergence. + A lower limit allows aborting faster. A higher limit may be needed + for integrands with many discontinuities or many singularities + close to the integration path. + This limit is only taken as a rough guideline, and the actual number of + function evaluations may be slightly higher depending on the + actual subdivisions. - - *flags* - additional options + .. member:: slong depth_limit - *ACB_CALC_VERBOSE* - print some information + Maximum search depth for adaptive subdivision. Technically, this is not + the limit on the local bisection depth but the limit on the number + of simultaneously queued subintervals. + If a zero or negative value is provided, the limit is set to the + default value `2 \cdot \text{prec}`. + Warning: memory usage may increase in proportion to this limit. - *ACB_CALC_VERY_VERBOSE* - print even more information + .. member:: int use_heap - *ACB_CALC_INTEGRATE_HEAP* - Use a heap to prioritize the segment - with the largest error globally. In this case *depth_limit* becomes - the maximum size of the heap. This sometimes gives better results - in case of convergence failure around an isolated point, but can + By default (if set to 0), new subintervals generated by adaptive + bisection will be appended to the top of a stack. + If set to 1, a binary heap will be used to maintain a priority queue + where the subintervals with larger error have higher priority. + This sometimes gives better results + in case of convergence failure, but can lead to a much larger array of subintervals (requiring a higher *depth_limit*) when many global bisections are needed. + .. member:: int verbose + + If set to 1, some information about the overall integration process + is printed to standard output. If set to 2, information about each + subinterval is printed. + +.. function:: void acb_calc_integrate_opt_init(acb_calc_integrate_opt_t options) + + Initializes *options* for use, setting all fields to 0 indicating + default values. + +Local integration algorithms +------------------------------------------------------------------------------- + .. function:: int acb_calc_integrate_gl_auto_deg(acb_t res, slong * num_eval, acb_calc_func_t func, void * param, const acb_t a, const acb_t b, const mag_t tol, slong deg_limit, int flags, slong prec) Attempts to compute `I = \int_a^b f(t) dt` using a single application @@ -240,8 +287,7 @@ Integration if `f` is holomorphic with `|f(z)| \le M` inside the ellipse *E* with foci `\pm 1` and semiaxes `X` and `Y = \sqrt{X^2 - 1}` such that `\rho = X + Y` - with `\rho > 1` - (See Trefethen, "Is Gauss Quadrature Better than Clenshaw-Curtis?"). + with `\rho > 1` [Tre2008]_. For an arbitrary interval, we use `\int_a^b f(t) dt = \int_{-1}^1 g(t) dt` where `g(t) = \Delta f(\Delta t + m)`, diff --git a/doc/source/credits.rst b/doc/source/credits.rst index a86be9a8..604d7e2a 100644 --- a/doc/source/credits.rst +++ b/doc/source/credits.rst @@ -226,4 +226,5 @@ Bibliography .. [Tak2000] \D. Takahashi, "A fast algorithm for computing large Fibonacci numbers", Information Processing Letters 75 (2000) 243-246, http://www.ii.uni.wroc.pl/~lorys/IPL/article75-6-1.pdf +.. [Tre2008] \L. N. Trefethen, "Is Gauss Quadrature Better than Clenshaw-Curtis?", SIAM Review, 50:1 (2008), 67-87, https://doi.org/10.1137/060659831 diff --git a/examples/integrals.c b/examples/integrals.c index 152d4262..b7b024f4 100644 --- a/examples/integrals.c +++ b/examples/integrals.c @@ -12,6 +12,7 @@ #include "arb_hypgeom.h" #include "acb_hypgeom.h" #include "acb_dirichlet.h" +#include "acb_modular.h" #include "acb_calc.h" /* ------------------------------------------------------------------------- */ @@ -155,7 +156,7 @@ f_rump(acb_ptr res, const acb_t z, void * param, slong order, slong prec) return 0; } -/* f(z) = |z^4 + 10z^3 + 19z^2 + 6z - 6| exp(z) (for real z) -- +/* f(z) = |z^4 + 10z^3 + 19z^2 - 6z - 6| exp(z) (for real z) -- Helfgott's integral on MathOverflow */ int f_helfgott(acb_ptr res, const acb_t z, void * param, slong order, slong prec) @@ -199,7 +200,7 @@ f_zeta(acb_ptr res, const acb_t z, void * param, slong order, slong prec) /* f(z) = z sin(1/z), assume on real interval */ int -f_essing(acb_ptr res, const acb_t z, void * param, slong order, slong prec) +f_essing2(acb_ptr res, const acb_t z, void * param, slong order, slong prec) { if (order > 1) flint_abort(); /* Would be needed for Taylor method. */ @@ -221,6 +222,28 @@ f_essing(acb_ptr res, const acb_t z, void * param, slong order, slong prec) return 0; } +/* f(z) = sin(1/z), assume on real interval */ +int +f_essing(acb_ptr res, const acb_t z, void * param, slong order, slong prec) +{ + if (order > 1) + flint_abort(); /* Would be needed for Taylor method. */ + + if ((order == 0) && acb_is_real(z) && arb_contains_zero(acb_realref(z))) + { + /* todo: arb_zero_pm_one, arb_unit_interval? */ + acb_zero(res); + mag_one(arb_radref(acb_realref(res))); + } + else + { + acb_inv(res, z, prec); + acb_sin(res, res, prec); + } + + return 0; +} + /* f(z) = exp(-z) z^1000 */ int f_factorial1000(acb_ptr res, const acb_t z, void * param, slong order, slong prec) @@ -388,11 +411,106 @@ f_sech3(acb_ptr res, const acb_t z, void * param, slong order, slong prec) return 0; } +/* f(z) = -log(z) / (1 + z) */ +int +f_log_div1p(acb_ptr res, const acb_t z, void * param, slong order, slong prec) +{ + acb_t t; + + if (order > 1) + flint_abort(); /* Would be needed for Taylor method. */ + + acb_init(t); + acb_add_ui(t, z, 1, prec); + acb_log(res, z, prec); + acb_div(res, res, t, prec); + acb_neg(res, res); + + acb_clear(t); + + return 0; +} + +/* f(z) = z exp(-z) / (1 + exp(-z)) */ +int +f_log_div1p_transformed(acb_ptr res, const acb_t z, void * param, slong order, slong prec) +{ + acb_t t; + + if (order > 1) + flint_abort(); /* Would be needed for Taylor method. */ + + acb_init(t); + acb_neg(t, z); + acb_exp(t, t, prec); + acb_add_ui(res, t, 1, prec); + acb_div(res, t, res, prec); + acb_mul(res, res, z, prec); + + acb_clear(t); + + return 0; +} + +int +f_elliptic_p_laurent_n(acb_ptr res, const acb_t z, void * param, slong order, slong prec) +{ + slong n; + acb_t tau; + + if (order > 1) + flint_abort(); /* Would be needed for Taylor method. */ + + n = ((slong *)(param))[0]; + + acb_init(tau); + acb_onei(tau); + + acb_modular_elliptic_p(res, z, tau, prec); + + if (n + 1 < 0) + { + acb_pow_ui(tau, z, -n - 1, prec); + acb_mul(res, res, tau, prec); + } + else /* note: acb_pow_si is currently bad; do it this way */ + { + acb_inv(tau, z, prec); + acb_pow_ui(tau, tau, n + 1, prec); + acb_mul(res, res, tau, prec); + } + + acb_clear(tau); + + return 0; +} + +/* f(z) = zeta'(z) / zeta(z) */ +int +f_zeta_frac(acb_ptr res, const acb_t z, void * param, slong order, slong prec) +{ + acb_struct t[2]; + + if (order > 1) + flint_abort(); /* Would be needed for Taylor method. */ + + acb_init(t); + acb_init(t + 1); + + acb_dirichlet_zeta_jet(t, z, 0, 2, prec); + acb_div(res, t + 1, t, prec); + + acb_clear(t); + acb_clear(t + 1); + + return 0; +} + /* ------------------------------------------------------------------------- */ /* Main test program */ /* ------------------------------------------------------------------------- */ -#define NUM_INTEGRALS 18 +#define NUM_INTEGRALS 23 const char * descr[NUM_INTEGRALS] = { @@ -404,6 +522,7 @@ const char * descr[NUM_INTEGRALS] = "int_0^100 floor(x) dx", "int_0^1 |x^4+10x^3+19x^2-6x-6| exp(x) dx", "1/(2 pi i) int zeta(s) ds (closed path around s = 1)", + "int_0^1 sin(1/x) dx (slow convergence, use -heap and/or -tol)", "int_0^1 x sin(1/x) dx (slow convergence, use -heap and/or -tol)", "int_0^10000 x^1000 exp(-x) dx", "int_1^{1+1000i} gamma(x) dx", @@ -414,33 +533,70 @@ const char * descr[NUM_INTEGRALS] = "int_0^8 (exp(x)-floor(exp(x))) sin(x+exp(x)) dx (use higher -eval)", "int_0^{inf} sech(x) dx (using domain truncation)", "int_0^{inf} sech^3(x) dx (using domain truncation)", + "int_0^1 -log(x)/(1+x) dx (using domain truncation)", + "int_0^{inf} x exp(-x)/(1+exp(-x)) dx (using domain truncation)", + "int_C wp(x)/x^(11) dx (contour for 10th Laurent coefficient of Weierstrass p-function)", + "N(1000) = count zeros with 0 < t <= 1000 of zeta(s) using argument principle" }; int main(int argc, char *argv[]) { acb_t s, t, a, b; mag_t tol; - slong prec, goal, deg_limit, eval_limit, depth_limit; + slong prec, goal; + slong N; int integral, ifrom, ito; - int i, twice, havetol, flags; + int i, twice, havegoal, havetol; + acb_calc_integrate_opt_t options; - flint_printf("Compute integrals using subdivision and Gauss-Legendre quadrature.\n"); - flint_printf("Usage: quadrature [-i n] [-prec p] [-tol eps] [-twice]\n\n"); - flint_printf("-i n - compute integral In (0 <= n <= %d)\n", NUM_INTEGRALS - 1); - flint_printf("-prec p - precision in bits (default p = 333)\n"); - flint_printf("-tol eps - approximate absolute error goal (default 2^-p)\n"); - flint_printf("-twice - run twice (to see overhead of computing nodes)\n"); - flint_printf("\n\n"); + ifrom = ito = -1; - prec = 333; + for (i = 1; i < argc; i++) + { + if (!strcmp(argv[i], "-i")) + { + if (!strcmp(argv[i+1], "all")) + { + ifrom = 0; + ito = NUM_INTEGRALS - 1; + } + else + { + ifrom = ito = atol(argv[i+1]); + if (ito < 0 || ito >= NUM_INTEGRALS) + flint_abort(); + } + } + } + + if (ifrom == -1) + { + flint_printf("Compute integrals using acb_calc_integrate.\n"); + flint_printf("Usage: integrals -i n [-prec p] [-tol eps] [-twice] [...]\n\n"); + flint_printf("-i n - compute integral n (0 <= n <= %d), or \"-i all\"\n", NUM_INTEGRALS - 1); + flint_printf("-prec p - precision in bits (default p = 64)\n"); + flint_printf("-goal p - approximate relative accuracy goal (default p)\n"); + flint_printf("-tol eps - approximate absolute error goal (default 2^-p)\n"); + flint_printf("-twice - run twice (to see overhead of computing nodes)\n"); + flint_printf("-heap - use heap for subinterval queue\n"); + flint_printf("-verbose - show information\n"); + flint_printf("-verbose2 - show more information\n"); + flint_printf("-deg n - use quadrature degree up to n\n"); + flint_printf("-eval n - limit number of function evaluations to n\n"); + flint_printf("-depth n - limit subinterval queue size to n\n\n"); + flint_printf("Implemented integrals:\n"); + for (integral = 0; integral < NUM_INTEGRALS; integral++) + flint_printf("I%d = %s\n", integral, descr[integral]); + flint_printf("\n"); + return 1; + } + + acb_calc_integrate_opt_init(options); + + prec = 64; twice = 0; - ifrom = 0; - ito = NUM_INTEGRALS - 1; - havetol = 0; - deg_limit = -1; - eval_limit = -1; - depth_limit = -1; - flags = 0; + goal = 0; + havetol = havegoal = 0; acb_init(a); acb_init(b); @@ -458,6 +614,16 @@ int main(int argc, char *argv[]) { twice = 1; } + else if (!strcmp(argv[i], "-goal")) + { + goal = atol(argv[i+1]); + if (goal < 0) + { + flint_printf("expected goal >= 0\n"); + return 1; + } + havegoal = 1; + } else if (!strcmp(argv[i], "-tol")) { arb_t x; @@ -467,39 +633,34 @@ int main(int argc, char *argv[]) arb_clear(x); havetol = 1; } - else if (!strcmp(argv[i], "-i")) - { - ifrom = ito = atol(argv[i+1]); - if (ito < 0 || ito >= NUM_INTEGRALS) - flint_abort(); - } else if (!strcmp(argv[i], "-deg")) { - deg_limit = atol(argv[i+1]); + options->deg_limit = atol(argv[i+1]); } else if (!strcmp(argv[i], "-eval")) { - eval_limit = atol(argv[i+1]); + options->eval_limit = atol(argv[i+1]); } else if (!strcmp(argv[i], "-depth")) { - depth_limit = atol(argv[i+1]); + options->depth_limit = atol(argv[i+1]); } else if (!strcmp(argv[i], "-verbose")) { - flags |= ACB_CALC_VERBOSE; + options->verbose = 1; } else if (!strcmp(argv[i], "-verbose2")) { - flags |= ACB_CALC_VERY_VERBOSE; + options->verbose = 2; } else if (!strcmp(argv[i], "-heap")) { - flags |= ACB_CALC_INTEGRATE_HEAP; + options->use_heap = 1; } } - goal = prec; + if (!havegoal) + goal = prec; if (!havetol) mag_set_ui_2exp_si(tol, 1, -prec); @@ -516,13 +677,13 @@ int main(int argc, char *argv[]) case 0: acb_set_d(a, 0); acb_set_d(b, 100); - acb_calc_integrate(s, f_sin, NULL, a, b, goal, tol, deg_limit, eval_limit, depth_limit, flags, prec); + acb_calc_integrate(s, f_sin, NULL, a, b, goal, tol, options, prec); break; case 1: acb_set_d(a, 0); acb_set_d(b, 1); - acb_calc_integrate(s, f_atanderiv, NULL, a, b, goal, tol, deg_limit, eval_limit, depth_limit, flags, prec); + acb_calc_integrate(s, f_atanderiv, NULL, a, b, goal, tol, options, prec); acb_mul_2exp_si(s, s, 2); break; @@ -530,7 +691,7 @@ int main(int argc, char *argv[]) acb_set_d(a, 0); acb_one(b); acb_mul_2exp_si(b, b, goal); - acb_calc_integrate(s, f_atanderiv, NULL, a, b, goal, tol, deg_limit, eval_limit, depth_limit, flags, prec); + acb_calc_integrate(s, f_atanderiv, NULL, a, b, goal, tol, options, prec); arb_add_error_2exp_si(acb_realref(s), -goal); acb_mul_2exp_si(s, s, 1); break; @@ -538,26 +699,26 @@ int main(int argc, char *argv[]) case 3: acb_set_d(a, 0); acb_set_d(b, 1); - acb_calc_integrate(s, f_circle, NULL, a, b, goal, tol, deg_limit, eval_limit, depth_limit, flags, prec); + acb_calc_integrate(s, f_circle, NULL, a, b, goal, tol, options, prec); acb_mul_2exp_si(s, s, 2); break; case 4: acb_set_d(a, 0); acb_set_d(b, 8); - acb_calc_integrate(s, f_rump, NULL, a, b, goal, tol, deg_limit, eval_limit, depth_limit, flags, prec); + acb_calc_integrate(s, f_rump, NULL, a, b, goal, tol, options, prec); break; case 5: - acb_set_d(a, 0); - acb_set_d(b, 100); - acb_calc_integrate(s, f_floor, NULL, a, b, goal, tol, deg_limit, eval_limit, depth_limit, flags, prec); + acb_set_d(a, 1); + acb_set_d(b, 101); + acb_calc_integrate(s, f_floor, NULL, a, b, goal, tol, options, prec); break; case 6: acb_set_d(a, 0); acb_set_d(b, 1); - acb_calc_integrate(s, f_helfgott, NULL, a, b, goal, tol, deg_limit, eval_limit, depth_limit, flags, prec); + acb_calc_integrate(s, f_helfgott, NULL, a, b, goal, tol, options, prec); break; case 7: @@ -565,22 +726,22 @@ int main(int argc, char *argv[]) acb_set_d_d(a, -1.0, -1.0); acb_set_d_d(b, 2.0, -1.0); - acb_calc_integrate(t, f_zeta, NULL, a, b, goal, tol, deg_limit, eval_limit, depth_limit, flags, prec); + acb_calc_integrate(t, f_zeta, NULL, a, b, goal, tol, options, prec); acb_add(s, s, t, prec); acb_set_d_d(a, 2.0, -1.0); acb_set_d_d(b, 2.0, 1.0); - acb_calc_integrate(t, f_zeta, NULL, a, b, goal, tol, deg_limit, eval_limit, depth_limit, flags, prec); + acb_calc_integrate(t, f_zeta, NULL, a, b, goal, tol, options, prec); acb_add(s, s, t, prec); acb_set_d_d(a, 2.0, 1.0); acb_set_d_d(b, -1.0, 1.0); - acb_calc_integrate(t, f_zeta, NULL, a, b, goal, tol, deg_limit, eval_limit, depth_limit, flags, prec); + acb_calc_integrate(t, f_zeta, NULL, a, b, goal, tol, options, prec); acb_add(s, s, t, prec); acb_set_d_d(a, -1.0, 1.0); acb_set_d_d(b, -1.0, -1.0); - acb_calc_integrate(t, f_zeta, NULL, a, b, goal, tol, deg_limit, eval_limit, depth_limit, flags, prec); + acb_calc_integrate(t, f_zeta, NULL, a, b, goal, tol, options, prec); acb_add(s, s, t, prec); acb_const_pi(t, prec); @@ -592,69 +753,75 @@ int main(int argc, char *argv[]) case 8: acb_set_d(a, 0); acb_set_d(b, 1); - acb_calc_integrate(s, f_essing, NULL, a, b, goal, tol, deg_limit, eval_limit, depth_limit, flags, prec); + acb_calc_integrate(s, f_essing, NULL, a, b, goal, tol, options, prec); break; case 9: acb_set_d(a, 0); - acb_set_d(b, 10000); - acb_calc_integrate(s, f_factorial1000, NULL, a, b, goal, tol, deg_limit, eval_limit, depth_limit, flags, prec); + acb_set_d(b, 1); + acb_calc_integrate(s, f_essing2, NULL, a, b, goal, tol, options, prec); break; case 10: - acb_set_d_d(a, 1.0, 0.0); - acb_set_d_d(b, 1.0, 1000.0); - acb_calc_integrate(s, f_gamma, NULL, a, b, goal, tol, deg_limit, eval_limit, depth_limit, flags, prec); + acb_set_d(a, 0); + acb_set_d(b, 10000); + acb_calc_integrate(s, f_factorial1000, NULL, a, b, goal, tol, options, prec); break; case 11: - acb_set_d(a, -10.0); - acb_set_d(b, 10.0); - acb_calc_integrate(s, f_sin_plus_small, NULL, a, b, goal, tol, deg_limit, eval_limit, depth_limit, flags, prec); + acb_set_d_d(a, 1.0, 0.0); + acb_set_d_d(b, 1.0, 1000.0); + acb_calc_integrate(s, f_gamma, NULL, a, b, goal, tol, options, prec); break; case 12: - acb_set_d(a, -1020.0); - acb_set_d(b, -1010.0); - acb_calc_integrate(s, f_exp, NULL, a, b, goal, tol, deg_limit, eval_limit, depth_limit, flags, prec); + acb_set_d(a, -10.0); + acb_set_d(b, 10.0); + acb_calc_integrate(s, f_sin_plus_small, NULL, a, b, goal, tol, options, prec); break; case 13: + acb_set_d(a, -1020.0); + acb_set_d(b, -1010.0); + acb_calc_integrate(s, f_exp, NULL, a, b, goal, tol, options, prec); + break; + + case 14: acb_set_d(a, 0); acb_set_d(b, ceil(sqrt(goal * 0.693147181) + 1.0)); - acb_calc_integrate(s, f_gaussian, NULL, a, b, goal, tol, deg_limit, eval_limit, depth_limit, flags, prec); + acb_calc_integrate(s, f_gaussian, NULL, a, b, goal, tol, options, prec); acb_mul(b, b, b, prec); acb_neg(b, b); acb_exp(b, b, prec); arb_add_error(acb_realref(s), acb_realref(b)); break; - case 14: - acb_set_d(a, 0.0); - acb_set_d(b, 1.0); - acb_calc_integrate(s, f_wolfram, NULL, a, b, goal, tol, deg_limit, eval_limit, depth_limit, flags, prec); - break; - case 15: acb_set_d(a, 0.0); - acb_set_d(b, 8.0); - acb_calc_integrate(s, f_monster, NULL, a, b, goal, tol, deg_limit, eval_limit, depth_limit, flags, prec); + acb_set_d(b, 1.0); + acb_calc_integrate(s, f_wolfram, NULL, a, b, goal, tol, options, prec); break; case 16: + acb_set_d(a, 0.0); + acb_set_d(b, 8.0); + acb_calc_integrate(s, f_monster, NULL, a, b, goal, tol, options, prec); + break; + + case 17: acb_set_d(a, 0); acb_set_d(b, ceil(goal * 0.693147181 + 1.0)); - acb_calc_integrate(s, f_sech, NULL, a, b, goal, tol, deg_limit, eval_limit, depth_limit, flags, prec); + acb_calc_integrate(s, f_sech, NULL, a, b, goal, tol, options, prec); acb_neg(b, b); acb_exp(b, b, prec); acb_mul_2exp_si(b, b, 1); arb_add_error(acb_realref(s), acb_realref(b)); break; - case 17: + case 18: acb_set_d(a, 0); acb_set_d(b, ceil(goal * 0.693147181 / 3.0 + 2.0)); - acb_calc_integrate(s, f_sech3, NULL, a, b, goal, tol, deg_limit, eval_limit, depth_limit, flags, prec); + acb_calc_integrate(s, f_sech3, NULL, a, b, goal, tol, options, prec); acb_neg(b, b); acb_mul_ui(b, b, 3, prec); acb_exp(b, b, prec); @@ -663,9 +830,98 @@ int main(int argc, char *argv[]) arb_add_error(acb_realref(s), acb_realref(b)); break; + case 19: + if (goal < 0) + abort(); + /* error bound 2^-N (1+N) when truncated at 2^-N */ + N = goal + FLINT_BIT_COUNT(goal); + acb_one(a); + acb_mul_2exp_si(a, a, -N); + acb_one(b); + acb_calc_integrate(s, f_log_div1p, NULL, a, b, goal, tol, options, prec); + acb_set_ui(b, N + 1); + acb_mul_2exp_si(b, b, -N); + arb_add_error(acb_realref(s), acb_realref(b)); + break; + + case 20: + if (goal < 0) + abort(); + /* error bound (N+1) exp(-N) when truncated at N */ + N = goal + FLINT_BIT_COUNT(goal); + acb_zero(a); + acb_set_ui(b, N); + acb_calc_integrate(s, f_log_div1p_transformed, NULL, a, b, goal, tol, options, prec); + acb_neg(b, b); + acb_exp(b, b, prec); + acb_mul_ui(b, b, N + 1, prec); + arb_add_error(acb_realref(s), acb_realref(b)); + break; + + case 21: + + acb_zero(s); + + N = 10; + + acb_set_d_d(a, 0.5, -0.5); + acb_set_d_d(b, 0.5, 0.5); + acb_calc_integrate(t, f_elliptic_p_laurent_n, &N, a, b, goal, tol, options, prec); + acb_add(s, s, t, prec); + + acb_set_d_d(a, 0.5, 0.5); + acb_set_d_d(b, -0.5, 0.5); + acb_calc_integrate(t, f_elliptic_p_laurent_n, &N, a, b, goal, tol, options, prec); + acb_add(s, s, t, prec); + + acb_set_d_d(a, -0.5, 0.5); + acb_set_d_d(b, -0.5, -0.5); + acb_calc_integrate(t, f_elliptic_p_laurent_n, &N, a, b, goal, tol, options, prec); + acb_add(s, s, t, prec); + + acb_set_d_d(a, -0.5, -0.5); + acb_set_d_d(b, 0.5, -0.5); + acb_calc_integrate(t, f_elliptic_p_laurent_n, &N, a, b, goal, tol, options, prec); + acb_add(s, s, t, prec); + + acb_const_pi(t, prec); + acb_div(s, s, t, prec); + acb_mul_2exp_si(s, s, -1); + acb_div_onei(s, s); + break; + + case 22: + + acb_zero(s); + + N = 1000; + + acb_set_d_d(a, 100.0, 0.0); + acb_set_d_d(b, 100.0, N); + acb_calc_integrate(t, f_zeta_frac, NULL, a, b, goal, tol, options, prec); + acb_add(s, s, t, prec); + + acb_set_d_d(a, 100, N); + acb_set_d_d(b, 0.5, N); + acb_calc_integrate(t, f_zeta_frac, NULL, a, b, goal, tol, options, prec); + acb_add(s, s, t, prec); + + acb_div_onei(s, s); + arb_zero(acb_imagref(s)); + + acb_set_ui(t, N); + acb_dirichlet_hardy_theta(t, t, NULL, NULL, 1, prec); + acb_add(s, s, t, prec); + + acb_const_pi(t, prec); + acb_div(s, s, t, prec); + acb_add_ui(s, s, 1, prec); + break; + default: abort(); } + TIMEIT_ONCE_STOP } flint_printf("I%d = ", integral);