diff --git a/acb_calc.h b/acb_calc.h index 7a723a35..3f6ba8bf 100644 --- a/acb_calc.h +++ b/acb_calc.h @@ -44,7 +44,7 @@ int acb_calc_integrate_taylor(acb_t res, const arf_t outer_radius, slong accuracy_goal, slong prec); -void +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, @@ -52,10 +52,11 @@ acb_calc_integrate(acb_t res, acb_calc_func_t f, void * param, int flags, slong prec); -slong -acb_calc_integrate_gl_auto_deg(acb_t res, 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); +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); #ifdef __cplusplus } diff --git a/acb_calc/integrate.c b/acb_calc/integrate.c index 5e58afe1..14237267 100644 --- a/acb_calc/integrate.c +++ b/acb_calc/integrate.c @@ -104,7 +104,7 @@ _acb_overlaps(acb_t tmp, const acb_t a, const acb_t b, slong prec) return acb_contains_zero(tmp); } -void +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, @@ -117,7 +117,10 @@ acb_calc_integrate(acb_t res, acb_calc_func_t f, void * param, acb_t s, t, u; mag_t tmpm, tmpn, new_tol; 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(t); @@ -155,6 +158,7 @@ acb_calc_integrate(acb_t res, acb_calc_func_t f, void * param, depth = depth_max = 1; eval = 1; stopping = 0; + leaf_interval_count = 0; /* Adjust absolute tolerance based on new information. */ 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) flint_printf("stopping at eval_limit %wd\n", eval_limit); + status = ARB_CALC_NO_CONVERGENCE; stopping = 1; 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_add(s, s, vs + top, prec); + leaf_interval_count++; + depth--; 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. */ 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); eval += feval; /* We are done with this subinterval. */ - if (feval > 0) + if (gl_status == ARB_CALC_SUCCESS) { if (real_error) arb_zero(acb_imagref(u)); acb_add(s, s, u, prec); + leaf_interval_count++; + /* Adjust absolute tolerance based on new information. */ acb_get_mag_lower(tmpm, u); 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) flint_printf("stopping at depth_limit %wd\n", depth_limit); + status = ARB_CALC_NO_CONVERGENCE; stopping = 1; continue; } @@ -247,14 +257,6 @@ acb_calc_integrate(acb_t res, acb_calc_func_t f, void * param, /* Interval [top] becomes [a, mid]. */ 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] */ 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))); @@ -294,8 +296,8 @@ acb_calc_integrate(acb_t res, acb_calc_func_t f, void * param, if (flags & ACB_CALC_VERBOSE) { - flint_printf("depth %wd/%wd, eval %wd/%wd\n", - depth_max, depth_limit, eval, eval_limit); + flint_printf("depth %wd/%wd, eval %wd/%wd, %wd leaf intervals\n", + depth_max, depth_limit, eval, eval_limit, leaf_interval_count); } 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(tmpn); mag_clear(new_tol); + + return status; } diff --git a/acb_calc/integrate_gl_auto_deg.c b/acb_calc/integrate_gl_auto_deg.c index cb0859f2..c1811909 100644 --- a/acb_calc/integrate_gl_auto_deg.c +++ b/acb_calc/integrate_gl_auto_deg.c @@ -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); } -slong -acb_calc_integrate_gl_auto_deg(acb_t res, 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) +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) { acb_t mid, delta, wide; mag_t tmpm; - slong success; + slong status; 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 i, n, best_n; - success = 0; + status = ARB_CALC_NO_CONVERGENCE; if (deg_limit <= 0) { acb_indeterminate(res); - return 0; + return status; } 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(t); mag_init(err); + mag_init(best_rho); best_n = -1; + eval_count[0] = 0; 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); f(v, wide, param, 1, prec); + eval_count[0]++; /* no chance */ 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) { - success = 1; + status = ARB_CALC_SUCCESS; /* The best so far. */ if (best_n == -1 || n < best_n) { mag_set(err, t); + mag_set(best_rho, rho); 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. */ - if (success) + if (status == ARB_CALC_SUCCESS) { arb_t x, w; arb_init(x); arb_init(w); - /* Return value. */ - success = best_n; - if ((flags & ACB_CALC_VERY_VERBOSE) == ACB_CALC_VERY_VERBOSE) { + acb_get_mag(tmpm, delta); flint_printf(" {GL deg %ld on [", best_n); - acb_printn(a, 20, ARB_STR_NO_RADIUS); flint_printf(", "); - acb_printn(b, 20, ARB_STR_NO_RADIUS); - flint_printf("], tol "); mag_printd(tol, 10); + acb_printn(a, 10, ARB_STR_NO_RADIUS); flint_printf(", "); + acb_printn(b, 10, ARB_STR_NO_RADIUS); + 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"); } @@ -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); } + eval_count[0] += best_n; + acb_mul(res, s, delta, prec); acb_add_error_mag(res, err); arb_clear(x); arb_clear(w); } + else + { + acb_indeterminate(res); + } acb_clear(s); 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(t); mag_clear(err); + mag_clear(best_rho); acb_clear(mid); acb_clear(delta); acb_clear(wide); mag_clear(tmpm); - return success; + return status; } diff --git a/doc/source/acb_calc.rst b/doc/source/acb_calc.rst index 11c3cb0e..ffcd9838 100644 --- a/doc/source/acb_calc.rst +++ b/doc/source/acb_calc.rst @@ -98,7 +98,7 @@ Types, macros and constants 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 @@ -114,6 +114,9 @@ Integration 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. + By default, the integrand *func* will only be called with *order* = 0 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 *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 of Gauss-Legendre quadrature with automatic determination of the quadrature degree so that the error is smaller than *tol*. - Returns a positive integer *n* indicating that the integral has - been evaluated successfully, or returns 0 if the tolerance could not be met. + Returns *ARB_CALC_SUCCESS* if the integral has been evaluated successfully + 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 rule is bounded by diff --git a/examples/integrals.c b/examples/integrals.c index 48c1db70..152d4262 100644 --- a/examples/integrals.c +++ b/examples/integrals.c @@ -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 #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); - acb_set(out, inp); - if (xlen > 1) - acb_one(out + 1); - _acb_poly_sin_series(out, out, xlen, order, prec); - return 0; + 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); + } + } } -int -elliptic(acb_ptr out, const acb_t inp, void * params, slong order, slong prec) +/* 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) { - 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); + 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 -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 z; - ulong n; + acb_t t; - t = _acb_vec_init(order); - acb_init(z); - - acb_set(t, inp); if (order > 1) - acb_one(t + 1); + flint_abort(); /* Would be needed for Taylor method. */ - n = 10; - arb_set_si(acb_realref(z), 20); - arb_set_si(acb_imagref(z), 10); + acb_init(t); - /* z sin(t) */ - _acb_poly_sin_series(out, t, FLINT_MIN(2, order), order, prec); - _acb_vec_scalar_mul(out, out, order, z, prec); + acb_exp(t, z, prec); + acb_holomorphic_floor(res, t, order != 0, prec); - /* t n */ - _acb_vec_scalar_mul_ui(t, t, FLINT_MIN(2, order), n, 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_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; } +/* 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[]) { - acb_t r, s, a, b; - arf_t inr, outr; - slong digits, prec; + 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; - if (argc < 2) - { - flint_printf("integrals d\n"); - flint_printf("compute integrals using d decimal digits of precision\n"); - return 1; - } + 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(r); - acb_init(s); acb_init(a); acb_init(b); - arf_init(inr); - arf_init(outr); + acb_init(s); + 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]); - prec = digits * 3.32193; - flint_printf("Digits: %wd\n", digits); + goal = prec; - 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 + if (!havetol) + mag_set_ui_2exp_si(tol, 1, -prec); - 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 + for (integral = ifrom; integral <= ito; integral++) + { + flint_printf("I%d = %s ...\n", integral, descr[integral]); - 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 + 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, 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(b); - arf_clear(inr); - arf_clear(outr); + acb_clear(s); + acb_clear(t); + mag_clear(tol); flint_cleanup(); return 0; diff --git a/examples/quadrature.c b/examples/quadrature.c deleted file mode 100644 index 1503517d..00000000 --- a/examples/quadrature.c +++ /dev/null @@ -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 -#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; -} - diff --git a/examples/taylor_integrals.c b/examples/taylor_integrals.c new file mode 100644 index 00000000..48c1db70 --- /dev/null +++ b/examples/taylor_integrals.c @@ -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; +} +