hiro/arb
1
0
Fork 0
mirror of https://github.com/vale981/arb synced 2025-03-06 01:41:39 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions
arb/examples/quadrature.c

528 lines
14 KiB
C
Raw Blame History

/*
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;
}
/* ------------------------------------------------------------------------- */
/* Main test program */
/* ------------------------------------------------------------------------- */
#define NUM_INTEGRALS 13
const char * descr[NUM_INTEGRALS] =
{
"int_0^100 sin(x) dx",
"4 int_0^1 1/(1+x^2) dx",
"4 int_1^{10^100} 1/(1+x^2) dx",
"4 int_0^1 sqrt(1-x^2) dx",
"4 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)",
"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 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;
}
}
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, 1);
acb_set_ui(b, 10);
acb_pow_ui(b, b, 100, prec);
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 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;
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;
}
Reference in a new issue View git blame Copy permalink