hiro/arb
1
0
Fork 0
mirror of https://github.com/vale981/arb synced 2025-03-04 08:51:40 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions

minor improvements to integration code; change interface; more examples

Browse source
This commit is contained in:
Fredrik Johansson 2017-11-22 00:24:20 +01:00
parent 75d8f23521
commit c9cca291d5
7 changed files with 539 additions and 175 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

30
acb_calc.h
View file

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

71
acb_calc/integrate.c
View file

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

5
acb_calc/integrate_gl_auto_deg.c
View file

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

23
acb_calc/integrate_opt_init.c Normal file
View file

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

188
doc/source/acb_calc.rst
View file

@ -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)`,

1
doc/source/credits.rst
View file

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

396
examples/integrals.c
View file

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