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

some changes to integration

Browse source
This commit is contained in:
Fredrik Johansson 2017-11-16 13:27:29 +01:00
parent a0ce21c37c
commit 0fbf524a22
7 changed files with 854 additions and 739 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

11
acb_calc.h
View file

@ -44,7 +44,7 @@ int acb_calc_integrate_taylor(acb_t res,
const arf_t outer_radius, const arf_t outer_radius,
slong accuracy_goal, slong prec); slong accuracy_goal, slong prec);
void int
acb_calc_integrate(acb_t res, acb_calc_func_t f, void * param, acb_calc_integrate(acb_t res, acb_calc_func_t f, void * param,
const acb_t a, const acb_t b, const acb_t a, const acb_t b,
slong goal, const mag_t tol, slong goal, const mag_t tol,
@ -52,10 +52,11 @@ acb_calc_integrate(acb_t res, acb_calc_func_t f, void * param,
int flags, int flags,
slong prec); slong prec);
slong int
acb_calc_integrate_gl_auto_deg(acb_t res, acb_calc_func_t f, void * param, acb_calc_integrate_gl_auto_deg(acb_t res, slong * eval_count,
const acb_t a, const acb_t b, const mag_t tol, acb_calc_func_t f, void * param,
slong deg_limit, int flags, slong prec); const acb_t a, const acb_t b, const mag_t tol,
slong deg_limit, int flags, slong prec);
#ifdef __cplusplus #ifdef __cplusplus
} }

32
acb_calc/integrate.c
View file

@ -104,7 +104,7 @@ _acb_overlaps(acb_t tmp, const acb_t a, const acb_t b, slong prec)
return acb_contains_zero(tmp); return acb_contains_zero(tmp);
} }
void int
acb_calc_integrate(acb_t res, acb_calc_func_t f, void * param, acb_calc_integrate(acb_t res, acb_calc_func_t f, void * param,
const acb_t a, const acb_t b, const acb_t a, const acb_t b,
slong goal, const mag_t tol, slong goal, const mag_t tol,
@ -117,7 +117,10 @@ acb_calc_integrate(acb_t res, acb_calc_func_t f, void * param,
acb_t s, t, u; acb_t s, t, u;
mag_t tmpm, tmpn, new_tol; mag_t tmpm, tmpn, new_tol;
slong depth, depth_max, eval, feval, top; slong depth, depth_max, eval, feval, top;
int stopping, real_error, use_heap; slong leaf_interval_count;
int stopping, real_error, use_heap, status, gl_status;
status = ARB_CALC_SUCCESS;
acb_init(s); acb_init(s);
acb_init(t); acb_init(t);
@ -155,6 +158,7 @@ acb_calc_integrate(acb_t res, acb_calc_func_t f, void * param,
depth = depth_max = 1; depth = depth_max = 1;
eval = 1; eval = 1;
stopping = 0; stopping = 0;
leaf_interval_count = 0;
/* Adjust absolute tolerance based on new information. */ /* Adjust absolute tolerance based on new information. */
acb_get_mag_lower(tmpm, vs); acb_get_mag_lower(tmpm, vs);
@ -169,6 +173,7 @@ acb_calc_integrate(acb_t res, acb_calc_func_t f, void * param,
{ {
if (flags & ACB_CALC_VERBOSE) if (flags & ACB_CALC_VERBOSE)
flint_printf("stopping at eval_limit %wd\n", eval_limit); flint_printf("stopping at eval_limit %wd\n", eval_limit);
status = ARB_CALC_NO_CONVERGENCE;
stopping = 1; stopping = 1;
continue; continue;
} }
@ -183,6 +188,8 @@ acb_calc_integrate(acb_t res, acb_calc_func_t f, void * param,
_acb_overlaps(u, as + top, bs + top, prec) || stopping) _acb_overlaps(u, as + top, bs + top, prec) || stopping)
{ {
acb_add(s, s, vs + top, prec); acb_add(s, s, vs + top, prec);
leaf_interval_count++;
depth--; depth--;
if (use_heap && depth > 0) if (use_heap && depth > 0)
{ {
@ -201,17 +208,19 @@ acb_calc_integrate(acb_t res, acb_calc_func_t f, void * param,
/* We know that the result is real. */ /* We know that the result is real. */
real_error = acb_is_finite(vs + top) && acb_is_real(vs + top); real_error = acb_is_finite(vs + top) && acb_is_real(vs + top);
feval = acb_calc_integrate_gl_auto_deg(u, f, param, 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, flags, prec);
eval += feval; eval += feval;
/* We are done with this subinterval. */ /* We are done with this subinterval. */
if (feval > 0) if (gl_status == ARB_CALC_SUCCESS)
{ {
if (real_error) if (real_error)
arb_zero(acb_imagref(u)); arb_zero(acb_imagref(u));
acb_add(s, s, u, prec); acb_add(s, s, u, prec);
leaf_interval_count++;
/* Adjust absolute tolerance based on new information. */ /* Adjust absolute tolerance based on new information. */
acb_get_mag_lower(tmpm, u); acb_get_mag_lower(tmpm, u);
mag_mul_2exp_si(tmpm, tmpm, -goal); mag_mul_2exp_si(tmpm, tmpm, -goal);
@ -234,6 +243,7 @@ acb_calc_integrate(acb_t res, acb_calc_func_t f, void * param,
{ {
if (flags & ACB_CALC_VERBOSE) if (flags & ACB_CALC_VERBOSE)
flint_printf("stopping at depth_limit %wd\n", depth_limit); flint_printf("stopping at depth_limit %wd\n", depth_limit);
status = ARB_CALC_NO_CONVERGENCE;
stopping = 1; stopping = 1;
continue; continue;
} }
@ -247,14 +257,6 @@ acb_calc_integrate(acb_t res, acb_calc_func_t f, void * param,
/* Interval [top] becomes [a, mid]. */ /* Interval [top] becomes [a, mid]. */
acb_set(bs + top, as + depth); acb_set(bs + top, as + depth);
#if 0
printf("new intervals:\n");
acb_printn(as, 10, ARB_STR_NO_RADIUS); printf(" ");
acb_printn(bs, 10, ARB_STR_NO_RADIUS); printf(", ");
acb_printn(as + depth, 10, ARB_STR_NO_RADIUS); printf(" ");
acb_printn(bs + depth, 10, ARB_STR_NO_RADIUS); printf("\n");
#endif
/* Evaluate on [a, mid] */ /* Evaluate on [a, mid] */
quad_simple(vs + top, f, param, as + top, bs + top, prec); quad_simple(vs + top, f, param, as + top, bs + top, prec);
mag_hypot(ms + top, arb_radref(acb_realref(vs + top)), arb_radref(acb_imagref(vs + top))); mag_hypot(ms + top, arb_radref(acb_realref(vs + top)), arb_radref(acb_imagref(vs + top)));
@ -294,8 +296,8 @@ acb_calc_integrate(acb_t res, acb_calc_func_t f, void * param,
if (flags & ACB_CALC_VERBOSE) if (flags & ACB_CALC_VERBOSE)
{ {
flint_printf("depth %wd/%wd, eval %wd/%wd\n", flint_printf("depth %wd/%wd, eval %wd/%wd, %wd leaf intervals\n",
depth_max, depth_limit, eval, eval_limit); depth_max, depth_limit, eval, eval_limit, leaf_interval_count);
} }
acb_set(res, s); acb_set(res, s);
@ -310,5 +312,7 @@ acb_calc_integrate(acb_t res, acb_calc_func_t f, void * param,
mag_clear(tmpm); mag_clear(tmpm);
mag_clear(tmpn); mag_clear(tmpn);
mag_clear(new_tol); mag_clear(new_tol);
return status;
} }

46
acb_calc/integrate_gl_auto_deg.c
View file

@ -111,25 +111,26 @@ acb_calc_gl_node(arb_t x, arb_t w, slong i, slong k, slong prec)
arb_set_round(w, gl_cache->gl_weights[i] + kk, prec); arb_set_round(w, gl_cache->gl_weights[i] + kk, prec);
} }
slong int
acb_calc_integrate_gl_auto_deg(acb_t res, acb_calc_func_t f, void * param, acb_calc_integrate_gl_auto_deg(acb_t res, slong * eval_count,
const acb_t a, const acb_t b, const mag_t tol, acb_calc_func_t f, void * param,
slong deg_limit, int flags, slong prec) const acb_t a, const acb_t b, const mag_t tol,
slong deg_limit, int flags, slong prec)
{ {
acb_t mid, delta, wide; acb_t mid, delta, wide;
mag_t tmpm; mag_t tmpm;
slong success; slong status;
acb_t s, v; acb_t s, v;
mag_t M, X, Y, rho, err, t; mag_t M, X, Y, rho, err, t, best_rho;
slong k, Xexp; slong k, Xexp;
slong i, n, best_n; slong i, n, best_n;
success = 0; status = ARB_CALC_NO_CONVERGENCE;
if (deg_limit <= 0) if (deg_limit <= 0)
{ {
acb_indeterminate(res); acb_indeterminate(res);
return 0; return status;
} }
acb_init(mid); acb_init(mid);
@ -153,8 +154,10 @@ acb_calc_integrate_gl_auto_deg(acb_t res, acb_calc_func_t f, void * param,
mag_init(rho); mag_init(rho);
mag_init(t); mag_init(t);
mag_init(err); mag_init(err);
mag_init(best_rho);
best_n = -1; best_n = -1;
eval_count[0] = 0;
mag_inf(err); mag_inf(err);
@ -185,6 +188,7 @@ acb_calc_integrate_gl_auto_deg(acb_t res, acb_calc_func_t f, void * param,
acb_add(wide, wide, mid, prec); acb_add(wide, wide, mid, prec);
f(v, wide, param, 1, prec); f(v, wide, param, 1, prec);
eval_count[0]++;
/* no chance */ /* no chance */
if (!acb_is_finite(v)) if (!acb_is_finite(v))
@ -211,12 +215,13 @@ acb_calc_integrate_gl_auto_deg(acb_t res, acb_calc_func_t f, void * param,
if (mag_cmp(t, tol) < 0) if (mag_cmp(t, tol) < 0)
{ {
success = 1; status = ARB_CALC_SUCCESS;
/* The best so far. */ /* The best so far. */
if (best_n == -1 || n < best_n) if (best_n == -1 || n < best_n)
{ {
mag_set(err, t); mag_set(err, t);
mag_set(best_rho, rho);
best_n = n; best_n = n;
} }
@ -228,21 +233,21 @@ acb_calc_integrate_gl_auto_deg(acb_t res, acb_calc_func_t f, void * param,
} }
/* Evaluate best found Gauss-Legendre quadrature rule. */ /* Evaluate best found Gauss-Legendre quadrature rule. */
if (success) if (status == ARB_CALC_SUCCESS)
{ {
arb_t x, w; arb_t x, w;
arb_init(x); arb_init(x);
arb_init(w); arb_init(w);
/* Return value. */
success = best_n;
if ((flags & ACB_CALC_VERY_VERBOSE) == ACB_CALC_VERY_VERBOSE) if ((flags & ACB_CALC_VERY_VERBOSE) == ACB_CALC_VERY_VERBOSE)
{ {
acb_get_mag(tmpm, delta);
flint_printf(" {GL deg %ld on [", best_n); flint_printf(" {GL deg %ld on [", best_n);
acb_printn(a, 20, ARB_STR_NO_RADIUS); flint_printf(", "); acb_printn(a, 10, ARB_STR_NO_RADIUS); flint_printf(", ");
acb_printn(b, 20, ARB_STR_NO_RADIUS); acb_printn(b, 10, ARB_STR_NO_RADIUS);
flint_printf("], tol "); mag_printd(tol, 10); flint_printf("], delta "); mag_printd(tmpm, 5);
flint_printf(", rho "); mag_printd(best_rho, 5);
flint_printf(", tol "); mag_printd(tol, 3);
flint_printf("}\n"); flint_printf("}\n");
} }
@ -264,12 +269,18 @@ acb_calc_integrate_gl_auto_deg(acb_t res, acb_calc_func_t f, void * param,
acb_addmul_arb(s, v, w, prec); acb_addmul_arb(s, v, w, prec);
} }
eval_count[0] += best_n;
acb_mul(res, s, delta, prec); acb_mul(res, s, delta, prec);
acb_add_error_mag(res, err); acb_add_error_mag(res, err);
arb_clear(x); arb_clear(x);
arb_clear(w); arb_clear(w);
} }
else
{
acb_indeterminate(res);
}
acb_clear(s); acb_clear(s);
acb_clear(v); acb_clear(v);
@ -279,12 +290,13 @@ acb_calc_integrate_gl_auto_deg(acb_t res, acb_calc_func_t f, void * param,
mag_clear(rho); mag_clear(rho);
mag_clear(t); mag_clear(t);
mag_clear(err); mag_clear(err);
mag_clear(best_rho);
acb_clear(mid); acb_clear(mid);
acb_clear(delta); acb_clear(delta);
acb_clear(wide); acb_clear(wide);
mag_clear(tmpm); mag_clear(tmpm);
return success; return status;
} }

12
doc/source/acb_calc.rst
View file

@ -98,7 +98,7 @@ Types, macros and constants
Integration Integration
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
.. function:: void 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 goal, const mag_t tol, slong deg_limit, slong eval_limit, slong depth_limit, int flags, slong prec)
Computes a rigorous enclosure of the integral Computes a rigorous enclosure of the integral
@ -114,6 +114,9 @@ Integration
of integration manually (or make a regularizing change of variables, of integration manually (or make a regularizing change of variables,
if possible). if possible).
Returns *ARB_CALC_SUCCESS* (0) if the integration converged on all
subintervals and *ARB_CALC_NO_CONVERGENCE* otherwise.
By default, the integrand *func* will only be called with *order* = 0 By default, the integrand *func* will only be called with *order* = 0
or *order* = 1; that is, derivatives are not required. or *order* = 1; that is, derivatives are not required.
@ -218,13 +221,14 @@ Integration
lead to a much larger array of subintervals (requiring a higher lead to a much larger array of subintervals (requiring a higher
*depth_limit*) when many global bisections are needed. *depth_limit*) when many global bisections are needed.
.. function:: slong acb_calc_integrate_gl_auto_deg(acb_t res, 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) .. 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 Attempts to compute `I = \int_a^b f(t) dt` using a single application
of Gauss-Legendre quadrature with automatic determination of the of Gauss-Legendre quadrature with automatic determination of the
quadrature degree so that the error is smaller than *tol*. quadrature degree so that the error is smaller than *tol*.
Returns a positive integer *n* indicating that the integral has Returns *ARB_CALC_SUCCESS* if the integral has been evaluated successfully
been evaluated successfully, or returns 0 if the tolerance could not be met. or *ARB_CALC_NO_CONVERGENCE* if the tolerance could not be met.
The total number of function evaluations is written to *num_eval*.
For the interval `[-1,1]`, the error of the *n*-point Gauss-Legendre For the interval `[-1,1]`, the error of the *n*-point Gauss-Legendre
rule is bounded by rule is bounded by

751
examples/integrals.c
View file

@ -1,148 +1,683 @@
/* This file is public domain. Author: Fredrik Johansson. */ /*
Rigorous numerical integration (with fast convergence for piecewise
holomorphic functions) using Gauss-Legendre quadrature and adaptive
subdivision.
#include "acb_calc.h" Author: Fredrik Johansson.
This file is in the public domain.
*/
#include <string.h>
#include "flint/profiler.h" #include "flint/profiler.h"
#include "arb_hypgeom.h"
#include "acb_hypgeom.h"
#include "acb_dirichlet.h"
#include "acb_calc.h"
int /* ------------------------------------------------------------------------- */
sinx(acb_ptr out, const acb_t inp, void * params, slong order, slong prec) /* Useful helper functions */
/* ------------------------------------------------------------------------- */
/* Absolute value function on R extended to a holomorphic function in the left
and right half planes. */
void
acb_holomorphic_abs(acb_ptr res, const acb_t z, int holomorphic, slong prec)
{ {
int xlen = FLINT_MIN(2, order); if (!acb_is_finite(z) || (holomorphic && arb_contains_zero(acb_realref(z))))
acb_set(out, inp); {
if (xlen > 1) acb_indeterminate(res);
acb_one(out + 1); }
_acb_poly_sin_series(out, out, xlen, order, prec); else
return 0; {
if (arb_is_nonnegative(acb_realref(z)))
{
acb_set_round(res, z, prec);
}
else if (arb_is_negative(acb_realref(z)))
{
acb_neg_round(res, z, prec);
}
else
{
acb_t t;
acb_init(t);
acb_neg(t, res);
acb_union(res, z, t, prec);
acb_clear(t);
}
}
} }
int /* Floor function on R extended to a piecewise holomorphic function in
elliptic(acb_ptr out, const acb_t inp, void * params, slong order, slong prec) vertical strips. */
void
acb_holomorphic_floor(acb_ptr res, const acb_t z, int holomorphic, slong prec)
{
if (!acb_is_finite(z) || (holomorphic && arb_contains_int(acb_realref(z))))
{
acb_indeterminate(res);
}
else
{
arb_floor(acb_realref(res), acb_realref(z), prec);
arb_set_round(acb_imagref(res), acb_imagref(z), prec);
}
}
/* Square root function on C with detection of the branch cut. */
void
acb_holomorphic_sqrt(acb_ptr res, const acb_t z, int holomorphic, slong prec)
{
if (!acb_is_finite(z) || (holomorphic &&
arb_contains_zero(acb_imagref(z)) &&
arb_contains_nonpositive(acb_realref(z))))
{
acb_indeterminate(res);
}
else
{
acb_sqrt(res, z, prec);
}
}
/* ------------------------------------------------------------------------- */
/* Example integrands */
/* ------------------------------------------------------------------------- */
/* f(z) = sin(z) */
int
f_sin(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
{ {
acb_ptr t;
t = _acb_vec_init(order);
acb_set(t, inp);
if (order > 1) if (order > 1)
acb_one(t + 1); flint_abort(); /* Would be needed for Taylor method. */
_acb_poly_sin_series(t, t, FLINT_MIN(2, order), order, prec);
_acb_poly_mullow(out, t, order, t, order, order, prec); acb_sin(res, z, prec);
_acb_vec_scalar_mul_2exp_si(t, out, order, -1);
acb_sub_ui(t, t, 1, prec); return 0;
_acb_vec_neg(t, t, order); }
_acb_poly_rsqrt_series(out, t, order, order, prec);
_acb_vec_clear(t, order); /* f(z) = floor(z) */
int
f_floor(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
{
if (order > 1)
flint_abort(); /* Would be needed for Taylor method. */
acb_holomorphic_floor(res, z, order != 0, prec);
return 0;
}
/* f(z) = sqrt(1-z^2) */
int
f_circle(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
{
if (order > 1)
flint_abort(); /* Would be needed for Taylor method. */
acb_one(res);
acb_submul(res, z, z, prec);
acb_holomorphic_sqrt(res, res, order != 0, prec);
/* Rounding could give |z| = 1 + eps near the endpoints, but we assume
that the interval is [-1,1] which really makes f real. */
if (order == 0)
arb_zero(acb_imagref(res));
return 0;
}
/* f(z) = 1/(1+z^2) */
int
f_atanderiv(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
{
if (order > 1)
flint_abort(); /* Would be needed for Taylor method. */
acb_mul(res, z, z, prec);
acb_add_ui(res, res, 1, prec);
acb_inv(res, res, prec);
return 0;
}
/* f(z) = sin(z + exp(z)) -- Rump's oscillatory example */
int
f_rump(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
{
if (order > 1)
flint_abort(); /* Would be needed for Taylor method. */
acb_exp(res, z, prec);
acb_add(res, res, z, prec);
acb_sin(res, res, prec);
return 0;
}
/* 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)
{
if (order > 1)
flint_abort(); /* Would be needed for Taylor method. */
acb_add_si(res, z, 10, prec);
acb_mul(res, res, z, prec);
acb_add_si(res, res, 19, prec);
acb_mul(res, res, z, prec);
acb_add_si(res, res, -6, prec);
acb_mul(res, res, z, prec);
acb_add_si(res, res, -6, prec);
acb_holomorphic_abs(res, res, order != 0, prec);
if (acb_is_finite(res))
{
acb_t t;
acb_init(t);
acb_exp(t, z, prec);
acb_mul(res, res, t, prec);
acb_clear(t);
}
return 0;
}
/* f(z) = zeta(z) */
int
f_zeta(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
{
if (order > 1)
flint_abort(); /* Would be needed for Taylor method. */
acb_zeta(res, z, prec);
return 0;
}
/* 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)
{
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);
}
acb_mul(res, res, z, 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)
{
acb_t t;
if (order > 1)
flint_abort(); /* Would be needed for Taylor method. */
acb_init(t);
acb_pow_ui(t, z, 1000, prec);
acb_neg(res, z);
acb_exp(res, res, prec);
acb_mul(res, res, t, prec);
acb_clear(t);
return 0;
}
/* f(z) = gamma(z) */
int
f_gamma(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
{
if (order > 1)
flint_abort(); /* Would be needed for Taylor method. */
acb_gamma(res, z, prec);
return 0;
}
/* f(z) = sin(z) + exp(-200-z^2) */
int
f_sin_plus_small(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_mul(t, z, z, prec);
acb_add_ui(t, t, 200, prec);
acb_neg(t, t);
acb_exp(t, t, prec);
acb_sin(res, z, prec);
acb_add(res, res, t, prec);
acb_clear(t);
return 0;
}
/* f(z) = exp(z) */
int
f_exp(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
{
if (order > 1)
flint_abort(); /* Would be needed for Taylor method. */
acb_exp(res, z, prec);
return 0;
}
/* f(z) = exp(-z^2) */
int
f_gaussian(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
{
if (order > 1)
flint_abort(); /* Would be needed for Taylor method. */
acb_mul(res, z, z, prec);
acb_neg(res, res);
acb_exp(res, res, prec);
return 0; return 0;
} }
int int
bessel(acb_ptr out, const acb_t inp, void * params, slong order, slong prec) f_monster(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
{ {
acb_ptr t; acb_t t;
acb_t z;
ulong n;
t = _acb_vec_init(order);
acb_init(z);
acb_set(t, inp);
if (order > 1) if (order > 1)
acb_one(t + 1); flint_abort(); /* Would be needed for Taylor method. */
n = 10; acb_init(t);
arb_set_si(acb_realref(z), 20);
arb_set_si(acb_imagref(z), 10);
/* z sin(t) */ acb_exp(t, z, prec);
_acb_poly_sin_series(out, t, FLINT_MIN(2, order), order, prec); acb_holomorphic_floor(res, t, order != 0, prec);
_acb_vec_scalar_mul(out, out, order, z, prec);
/* t n */ if (acb_is_finite(res))
_acb_vec_scalar_mul_ui(t, t, FLINT_MIN(2, order), n, prec); {
acb_sub(res, t, res, prec);
acb_add(t, t, z, prec);
acb_sin(t, t, prec);
acb_mul(res, res, t, prec);
}
_acb_poly_sub(out, t, FLINT_MIN(2, order), out, order, prec); acb_clear(t);
_acb_poly_cos_series(out, out, order, order, prec);
_acb_vec_clear(t, order);
acb_clear(z);
return 0; return 0;
} }
/* f(z) = sech(10(x-0.2))^2 + sech(100(x-0.4))^4 + sech(1000(x-0.6))^6 */
int
f_wolfram(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
{
acb_t a, b, c;
if (order > 1)
flint_abort(); /* Would be needed for Taylor method. */
acb_init(a);
acb_init(b);
acb_init(c);
acb_mul_ui(a, z, 10, prec);
acb_sub_ui(a, a, 2, prec);
acb_sech(a, a, prec);
acb_pow_ui(a, a, 2, prec);
acb_mul_ui(b, z, 100, prec);
acb_sub_ui(b, b, 40, prec);
acb_sech(b, b, prec);
acb_pow_ui(b, b, 4, prec);
acb_mul_ui(c, z, 1000, prec);
acb_sub_ui(c, c, 600, prec);
acb_sech(c, c, prec);
acb_pow_ui(c, c, 6, prec);
acb_add(res, a, b, prec);
acb_add(res, res, c, prec);
acb_clear(a);
acb_clear(b);
acb_clear(c);
return 0;
}
/* f(z) = sech(z) */
int
f_sech(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
{
if (order > 1)
flint_abort(); /* Would be needed for Taylor method. */
acb_sech(res, z, prec);
return 0;
}
/* f(z) = sech^3(z) */
int
f_sech3(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
{
if (order > 1)
flint_abort(); /* Would be needed for Taylor method. */
acb_sech(res, z, prec);
acb_cube(res, res, prec);
return 0;
}
/* ------------------------------------------------------------------------- */
/* Main test program */
/* ------------------------------------------------------------------------- */
#define NUM_INTEGRALS 18
const char * descr[NUM_INTEGRALS] =
{
"int_0^100 sin(x) dx",
"4 int_0^1 1/(1+x^2) dx",
"2 int_0^{inf} 1/(1+x^2) dx (using domain truncation)",
"4 int_0^1 sqrt(1-x^2) dx",
"int_0^8 sin(x+exp(x)) dx",
"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 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",
"int_{-10}^{10} sin(x) + exp(-200-x^2) dx",
"int_{-1020}^{-1000} exp(x) dx (use -tol 0 for relative error)",
"int_0^{inf} exp(-x^2) dx (using domain truncation)",
"int_0^1 sech(10(x-0.2))^2 + sech(100(x-0.4))^4 + sech(1000(x-0.6))^6 dx",
"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 main(int argc, char *argv[]) int main(int argc, char *argv[])
{ {
acb_t r, s, a, b; acb_t s, t, a, b;
arf_t inr, outr; mag_t tol;
slong digits, prec; slong prec, goal, deg_limit, eval_limit, depth_limit;
int integral, ifrom, ito;
int i, twice, havetol, flags;
if (argc < 2) 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("integrals d\n"); flint_printf("-i n - compute integral In (0 <= n <= %d)\n", NUM_INTEGRALS - 1);
flint_printf("compute integrals using d decimal digits of precision\n"); flint_printf("-prec p - precision in bits (default p = 333)\n");
return 1; 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");
prec = 333;
twice = 0;
ifrom = 0;
ito = NUM_INTEGRALS - 1;
havetol = 0;
deg_limit = -1;
eval_limit = -1;
depth_limit = -1;
flags = 0;
acb_init(r);
acb_init(s);
acb_init(a); acb_init(a);
acb_init(b); acb_init(b);
arf_init(inr); acb_init(s);
arf_init(outr); acb_init(t);
mag_init(tol);
arb_calc_verbose = 0; for (i = 1; i < argc; i++)
{
if (!strcmp(argv[i], "-prec"))
{
prec = atol(argv[i+1]);
}
else if (!strcmp(argv[i], "-twice"))
{
twice = 1;
}
else if (!strcmp(argv[i], "-tol"))
{
arb_t x;
arb_init(x);
arb_set_str(x, argv[i+1], 10);
arb_get_mag(tol, x);
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]);
}
else if (!strcmp(argv[i], "-eval"))
{
eval_limit = atol(argv[i+1]);
}
else if (!strcmp(argv[i], "-depth"))
{
depth_limit = atol(argv[i+1]);
}
else if (!strcmp(argv[i], "-verbose"))
{
flags |= ACB_CALC_VERBOSE;
}
else if (!strcmp(argv[i], "-verbose2"))
{
flags |= ACB_CALC_VERY_VERBOSE;
}
else if (!strcmp(argv[i], "-heap"))
{
flags |= ACB_CALC_INTEGRATE_HEAP;
}
}
digits = atol(argv[1]); goal = prec;
prec = digits * 3.32193;
flint_printf("Digits: %wd\n", digits);
flint_printf("----------------------------------------------------------------\n"); if (!havetol)
flint_printf("Integral of sin(t) from 0 to 100.\n"); mag_set_ui_2exp_si(tol, 1, -prec);
arf_set_d(inr, 0.125);
arf_set_d(outr, 1.0);
TIMEIT_ONCE_START
acb_set_si(a, 0);
acb_set_si(b, 100);
acb_calc_integrate_taylor(r, sinx, NULL, a, b, inr, outr, prec, 1.1 * prec);
flint_printf("RESULT:\n");
acb_printn(r, digits, 0); flint_printf("\n");
TIMEIT_ONCE_STOP
flint_printf("----------------------------------------------------------------\n"); for (integral = ifrom; integral <= ito; integral++)
flint_printf("Elliptic integral F(phi, m) = integral of 1/sqrt(1 - m*sin(t)^2)\n"); {
flint_printf("from 0 to phi, with phi = 6+6i, m = 1/2. Integration path\n"); flint_printf("I%d = %s ...\n", integral, descr[integral]);
flint_printf("0 -> 6 -> 6+6i.\n");
arf_set_d(inr, 0.2);
arf_set_d(outr, 0.5);
TIMEIT_ONCE_START
acb_set_si(a, 0);
acb_set_si(b, 6);
acb_calc_integrate_taylor(r, elliptic, NULL, a, b, inr, outr, prec, 1.1 * prec);
acb_set_si(a, 6);
arb_set_si(acb_realref(b), 6);
arb_set_si(acb_imagref(b), 6);
acb_calc_integrate_taylor(s, elliptic, NULL, a, b, inr, outr, prec, 1.1 * prec);
acb_add(r, r, s, prec);
flint_printf("RESULT:\n");
acb_printn(r, digits, 0); flint_printf("\n");
TIMEIT_ONCE_STOP
flint_printf("----------------------------------------------------------------\n"); for (i = 0; i < 1 + twice; i++)
flint_printf("Bessel function J_n(z) = (1/pi) * integral of cos(t*n - z*sin(t))\n"); {
flint_printf("from 0 to pi. With n = 10, z = 20 + 10i.\n"); TIMEIT_ONCE_START
arf_set_d(inr, 0.1); switch (integral)
arf_set_d(outr, 0.5); {
TIMEIT_ONCE_START case 0:
acb_set_si(a, 0); acb_set_d(a, 0);
acb_const_pi(b, 3 * prec); acb_set_d(b, 100);
acb_calc_integrate_taylor(r, bessel, NULL, a, b, inr, outr, prec, 3 * prec); acb_calc_integrate(s, f_sin, NULL, a, b, goal, tol, deg_limit, eval_limit, depth_limit, flags, prec);
acb_div(r, r, b, prec); break;
flint_printf("RESULT:\n");
acb_printn(r, digits, 0); flint_printf("\n"); case 1:
TIMEIT_ONCE_STOP 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_mul_2exp_si(s, s, 2);
break;
case 2:
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);
arb_add_error_2exp_si(acb_realref(s), -goal);
acb_mul_2exp_si(s, s, 1);
break;
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_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);
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);
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);
break;
case 7:
acb_zero(s);
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_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_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_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_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 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);
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);
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);
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);
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);
break;
case 13:
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_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);
break;
case 16:
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_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:
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_neg(b, b);
acb_mul_ui(b, b, 3, prec);
acb_exp(b, b, prec);
acb_mul_2exp_si(b, b, 3);
acb_div_ui(b, b, 3, prec);
arb_add_error(acb_realref(s), acb_realref(b));
break;
default:
abort();
}
TIMEIT_ONCE_STOP
}
flint_printf("I%d = ", integral);
acb_printn(s, 3.333 * prec, 0);
flint_printf("\n\n");
}
acb_clear(r);
acb_clear(s);
acb_clear(a); acb_clear(a);
acb_clear(b); acb_clear(b);
arf_clear(inr); acb_clear(s);
arf_clear(outr); acb_clear(t);
mag_clear(tol);
flint_cleanup(); flint_cleanup();
return 0; return 0;

591
examples/quadrature.c
View file

@ -1,591 +0,0 @@
/*
Rigorous numerical integration (with fast convergence for piecewise
holomorphic functions) using Gauss-Legendre quadrature and adaptive
subdivision.
Author: Fredrik Johansson.
This file is in the public domain.
*/
#include <string.h>
#include "flint/profiler.h"
#include "arb_hypgeom.h"
#include "acb_hypgeom.h"
#include "acb_dirichlet.h"
#include "acb_calc.h"
/* ------------------------------------------------------------------------- */
/* Useful helper functions */
/* ------------------------------------------------------------------------- */
/* Absolute value function on R extended to a holomorphic function in the left
and right half planes. */
void
acb_holomorphic_abs(acb_ptr res, const acb_t z, int holomorphic, slong prec)
{
if (!acb_is_finite(z) || (holomorphic && arb_contains_zero(acb_realref(z))))
{
acb_indeterminate(res);
}
else
{
if (arb_is_nonnegative(acb_realref(z)))
{
acb_set_round(res, z, prec);
}
else if (arb_is_negative(acb_realref(z)))
{
acb_neg_round(res, z, prec);
}
else
{
acb_t t;
acb_init(t);
acb_neg(t, res);
acb_union(res, z, t, prec);
acb_clear(t);
}
}
}
/* Floor function on R extended to a piecewise holomorphic function in
vertical strips. */
void
acb_holomorphic_floor(acb_ptr res, const acb_t z, int holomorphic, slong prec)
{
if (!acb_is_finite(z) || (holomorphic && arb_contains_int(acb_realref(z))))
{
acb_indeterminate(res);
}
else
{
arb_floor(acb_realref(res), acb_realref(z), prec);
arb_set_round(acb_imagref(res), acb_imagref(z), prec);
}
}
/* Square root function on C with detection of the branch cut. */
void
acb_holomorphic_sqrt(acb_ptr res, const acb_t z, int holomorphic, slong prec)
{
if (!acb_is_finite(z) || (holomorphic &&
arb_contains_zero(acb_imagref(z)) &&
arb_contains_nonpositive(acb_realref(z))))
{
acb_indeterminate(res);
}
else
{
acb_sqrt(res, z, prec);
}
}
/* ------------------------------------------------------------------------- */
/* Example integrands */
/* ------------------------------------------------------------------------- */
/* f(z) = sin(z) */
int
f_sin(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
{
if (order > 1)
flint_abort(); /* Would be needed for Taylor method. */
acb_sin(res, z, prec);
return 0;
}
/* f(z) = floor(z) */
int
f_floor(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
{
if (order > 1)
flint_abort(); /* Would be needed for Taylor method. */
acb_holomorphic_floor(res, z, order != 0, prec);
return 0;
}
/* f(z) = sqrt(1-z^2) */
int
f_circle(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
{
if (order > 1)
flint_abort(); /* Would be needed for Taylor method. */
acb_one(res);
acb_submul(res, z, z, prec);
acb_holomorphic_sqrt(res, res, order != 0, prec);
/* Rounding could give |z| = 1 + eps near the endpoints, but we assume
that the interval is [-1,1] which really makes f real. */
if (order == 0)
arb_zero(acb_imagref(res));
return 0;
}
/* f(z) = 1/(1+z^2) */
int
f_atanderiv(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
{
if (order > 1)
flint_abort(); /* Would be needed for Taylor method. */
acb_mul(res, z, z, prec);
acb_add_ui(res, res, 1, prec);
acb_inv(res, res, prec);
return 0;
}
/* f(z) = sin(z + exp(z)) -- Rump's oscillatory example */
int
f_rump(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
{
if (order > 1)
flint_abort(); /* Would be needed for Taylor method. */
acb_exp(res, z, prec);
acb_add(res, res, z, prec);
acb_sin(res, res, prec);
return 0;
}
/* 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)
{
if (order > 1)
flint_abort(); /* Would be needed for Taylor method. */
acb_add_si(res, z, 10, prec);
acb_mul(res, res, z, prec);
acb_add_si(res, res, 19, prec);
acb_mul(res, res, z, prec);
acb_add_si(res, res, -6, prec);
acb_mul(res, res, z, prec);
acb_add_si(res, res, -6, prec);
acb_holomorphic_abs(res, res, order != 0, prec);
if (acb_is_finite(res))
{
acb_t t;
acb_init(t);
acb_exp(t, z, prec);
acb_mul(res, res, t, prec);
acb_clear(t);
}
return 0;
}
/* f(z) = zeta(z) */
int
f_zeta(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
{
if (order > 1)
flint_abort(); /* Would be needed for Taylor method. */
acb_zeta(res, z, prec);
return 0;
}
/* 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)
{
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);
}
acb_mul(res, res, z, 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)
{
acb_t t;
if (order > 1)
flint_abort(); /* Would be needed for Taylor method. */
acb_init(t);
acb_pow_ui(t, z, 1000, prec);
acb_neg(res, z);
acb_exp(res, res, prec);
acb_mul(res, res, t, prec);
acb_clear(t);
return 0;
}
/* f(z) = gamma(z) */
int
f_gamma(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
{
if (order > 1)
flint_abort(); /* Would be needed for Taylor method. */
acb_gamma(res, z, prec);
return 0;
}
/* f(z) = sin(z) + exp(-200-z^2) */
int
f_sin_plus_small(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_mul(t, z, z, prec);
acb_add_ui(t, t, 200, prec);
acb_neg(t, t);
acb_exp(t, t, prec);
acb_sin(res, z, prec);
acb_add(res, res, t, prec);
acb_clear(t);
return 0;
}
/* f(z) = exp(z) */
int
f_exp(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
{
if (order > 1)
flint_abort(); /* Would be needed for Taylor method. */
acb_exp(res, z, prec);
return 0;
}
/* f(z) = exp(-z^2) */
int
f_gaussian(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
{
if (order > 1)
flint_abort(); /* Would be needed for Taylor method. */
acb_mul(res, z, z, prec);
acb_neg(res, res);
acb_exp(res, res, prec);
return 0;
}
int
f_monster(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_exp(t, z, prec);
acb_holomorphic_floor(res, t, order != 0, prec);
if (acb_is_finite(res))
{
acb_sub(res, t, res, prec);
acb_add(t, t, z, prec);
acb_sin(t, t, prec);
acb_mul(res, res, t, prec);
}
acb_clear(t);
return 0;
}
/* ------------------------------------------------------------------------- */
/* Main test program */
/* ------------------------------------------------------------------------- */
#define NUM_INTEGRALS 15
const char * descr[NUM_INTEGRALS] =
{
"int_0^100 sin(x) dx",
"4 int_0^1 1/(1+x^2) dx",
"2 int_0^{inf} 1/(1+x^2) dx (using domain truncation)",
"4 int_0^1 sqrt(1-x^2) dx",
"int_0^8 sin(x+exp(x)) dx",
"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 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",
"int_{-10}^{10} sin(x) + exp(-200-x^2) dx",
"int_{-1020}^{-1000} exp(x) dx (use -tol 0 for relative error)",
"int_0^{inf} exp(-x^2) dx (using domain truncation)",
"int_0^8 (exp(x)-floor(exp(x))) sin(x+exp(x)) dx (use higher -eval)",
};
int main(int argc, char *argv[])
{
acb_t s, t, a, b;
mag_t tol;
slong prec, goal, deg_limit, eval_limit, depth_limit;
int integral, ifrom, ito;
int i, twice, havetol, flags;
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");
prec = 333;
twice = 0;
ifrom = 0;
ito = NUM_INTEGRALS - 1;
havetol = 0;
deg_limit = -1;
eval_limit = -1;
depth_limit = -1;
flags = 0;
acb_init(a);
acb_init(b);
acb_init(s);
acb_init(t);
mag_init(tol);
for (i = 1; i < argc; i++)
{
if (!strcmp(argv[i], "-prec"))
{
prec = atol(argv[i+1]);
}
else if (!strcmp(argv[i], "-twice"))
{
twice = 1;
}
else if (!strcmp(argv[i], "-tol"))
{
arb_t x;
arb_init(x);
arb_set_str(x, argv[i+1], 10);
arb_get_mag(tol, x);
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]);
}
else if (!strcmp(argv[i], "-eval"))
{
eval_limit = atol(argv[i+1]);
}
else if (!strcmp(argv[i], "-depth"))
{
depth_limit = atol(argv[i+1]);
}
else if (!strcmp(argv[i], "-verbose"))
{
flags |= ACB_CALC_VERBOSE;
}
else if (!strcmp(argv[i], "-verbose2"))
{
flags |= ACB_CALC_VERY_VERBOSE;
}
else if (!strcmp(argv[i], "-heap"))
{
flags |= ACB_CALC_INTEGRATE_HEAP;
}
}
goal = prec;
if (!havetol)
mag_set_ui_2exp_si(tol, 1, -prec);
for (integral = ifrom; integral <= ito; integral++)
{
flint_printf("I%d = %s ...\n", integral, descr[integral]);
for (i = 0; i < 1 + twice; i++)
{
TIMEIT_ONCE_START
switch (integral)
{
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);
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_mul_2exp_si(s, s, 2);
break;
case 2:
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);
arb_add_error_2exp_si(acb_realref(s), -goal);
acb_mul_2exp_si(s, s, 1);
break;
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_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);
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);
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);
break;
case 7:
acb_zero(s);
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_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_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_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_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 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);
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);
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);
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);
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);
break;
case 13:
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_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, 8.0);
acb_calc_integrate(s, f_monster, NULL, a, b, goal, tol, deg_limit, eval_limit, depth_limit, flags, prec);
break;
default:
abort();
}
TIMEIT_ONCE_STOP
}
flint_printf("I%d = ", integral);
acb_printn(s, 3.333 * prec, 0);
flint_printf("\n\n");
}
acb_clear(a);
acb_clear(b);
acb_clear(s);
acb_clear(t);
mag_clear(tol);
flint_cleanup();
return 0;
}

150
examples/taylor_integrals.c Normal file
View file

@ -0,0 +1,150 @@
/* This file is public domain. Author: Fredrik Johansson. */
#include "acb_calc.h"
#include "flint/profiler.h"
int
sinx(acb_ptr out, const acb_t inp, void * params, slong order, slong prec)
{
int xlen = FLINT_MIN(2, order);
acb_set(out, inp);
if (xlen > 1)
acb_one(out + 1);
_acb_poly_sin_series(out, out, xlen, order, prec);
return 0;
}
int
elliptic(acb_ptr out, const acb_t inp, void * params, slong order, slong prec)
{
acb_ptr t;
t = _acb_vec_init(order);
acb_set(t, inp);
if (order > 1)
acb_one(t + 1);
_acb_poly_sin_series(t, t, FLINT_MIN(2, order), order, prec);
_acb_poly_mullow(out, t, order, t, order, order, prec);
_acb_vec_scalar_mul_2exp_si(t, out, order, -1);
acb_sub_ui(t, t, 1, prec);
_acb_vec_neg(t, t, order);
_acb_poly_rsqrt_series(out, t, order, order, prec);
_acb_vec_clear(t, order);
return 0;
}
int
bessel(acb_ptr out, const acb_t inp, void * params, slong order, slong prec)
{
acb_ptr t;
acb_t z;
ulong n;
t = _acb_vec_init(order);
acb_init(z);
acb_set(t, inp);
if (order > 1)
acb_one(t + 1);
n = 10;
arb_set_si(acb_realref(z), 20);
arb_set_si(acb_imagref(z), 10);
/* z sin(t) */
_acb_poly_sin_series(out, t, FLINT_MIN(2, order), order, prec);
_acb_vec_scalar_mul(out, out, order, z, prec);
/* t n */
_acb_vec_scalar_mul_ui(t, t, FLINT_MIN(2, order), n, prec);
_acb_poly_sub(out, t, FLINT_MIN(2, order), out, order, prec);
_acb_poly_cos_series(out, out, order, order, prec);
_acb_vec_clear(t, order);
acb_clear(z);
return 0;
}
int main(int argc, char *argv[])
{
acb_t r, s, a, b;
arf_t inr, outr;
slong digits, prec;
if (argc < 2)
{
flint_printf("integrals d\n");
flint_printf("compute integrals using d decimal digits of precision\n");
return 1;
}
acb_init(r);
acb_init(s);
acb_init(a);
acb_init(b);
arf_init(inr);
arf_init(outr);
arb_calc_verbose = 0;
digits = atol(argv[1]);
prec = digits * 3.32193;
flint_printf("Digits: %wd\n", digits);
flint_printf("----------------------------------------------------------------\n");
flint_printf("Integral of sin(t) from 0 to 100.\n");
arf_set_d(inr, 0.125);
arf_set_d(outr, 1.0);
TIMEIT_ONCE_START
acb_set_si(a, 0);
acb_set_si(b, 100);
acb_calc_integrate_taylor(r, sinx, NULL, a, b, inr, outr, prec, 1.1 * prec);
flint_printf("RESULT:\n");
acb_printn(r, digits, 0); flint_printf("\n");
TIMEIT_ONCE_STOP
flint_printf("----------------------------------------------------------------\n");
flint_printf("Elliptic integral F(phi, m) = integral of 1/sqrt(1 - m*sin(t)^2)\n");
flint_printf("from 0 to phi, with phi = 6+6i, m = 1/2. Integration path\n");
flint_printf("0 -> 6 -> 6+6i.\n");
arf_set_d(inr, 0.2);
arf_set_d(outr, 0.5);
TIMEIT_ONCE_START
acb_set_si(a, 0);
acb_set_si(b, 6);
acb_calc_integrate_taylor(r, elliptic, NULL, a, b, inr, outr, prec, 1.1 * prec);
acb_set_si(a, 6);
arb_set_si(acb_realref(b), 6);
arb_set_si(acb_imagref(b), 6);
acb_calc_integrate_taylor(s, elliptic, NULL, a, b, inr, outr, prec, 1.1 * prec);
acb_add(r, r, s, prec);
flint_printf("RESULT:\n");
acb_printn(r, digits, 0); flint_printf("\n");
TIMEIT_ONCE_STOP
flint_printf("----------------------------------------------------------------\n");
flint_printf("Bessel function J_n(z) = (1/pi) * integral of cos(t*n - z*sin(t))\n");
flint_printf("from 0 to pi. With n = 10, z = 20 + 10i.\n");
arf_set_d(inr, 0.1);
arf_set_d(outr, 0.5);
TIMEIT_ONCE_START
acb_set_si(a, 0);
acb_const_pi(b, 3 * prec);
acb_calc_integrate_taylor(r, bessel, NULL, a, b, inr, outr, prec, 3 * prec);
acb_div(r, r, b, prec);
flint_printf("RESULT:\n");
acb_printn(r, digits, 0); flint_printf("\n");
TIMEIT_ONCE_STOP
acb_clear(r);
acb_clear(s);
acb_clear(a);
acb_clear(b);
arf_clear(inr);
arf_clear(outr);
flint_cleanup();
return 0;
}