mirror of
https://github.com/vale981/arb
synced 2025-03-06 01:41:39 -05:00
719 lines
21 KiB
C
719 lines
21 KiB
C
/*
|
|
Copyright (C) 2017 Fredrik Johansson
|
|
|
|
This file is part of Arb.
|
|
|
|
Arb is free software: you can redistribute it and/or modify it under
|
|
the terms of the GNU Lesser General Public License (LGPL) as published
|
|
by the Free Software Foundation; either version 2.1 of the License, or
|
|
(at your option) any later version. See <http://www.gnu.org/licenses/>.
|
|
*/
|
|
|
|
#include "acb_calc.h"
|
|
|
|
int
|
|
f_zero(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_zero(res);
|
|
|
|
return 0;
|
|
}
|
|
|
|
int
|
|
f_indet(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_indeterminate(res);
|
|
|
|
return 0;
|
|
}
|
|
|
|
int
|
|
f_cube(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_cube(res, z, prec);
|
|
|
|
return 0;
|
|
}
|
|
|
|
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;
|
|
}
|
|
|
|
int
|
|
f_sqrt(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_sqrt_analytic(res, z, order != 0, prec);
|
|
|
|
return 0;
|
|
}
|
|
|
|
int
|
|
f_rsqrt(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_rsqrt_analytic(res, z, order != 0, prec);
|
|
|
|
return 0;
|
|
}
|
|
|
|
int
|
|
f_log(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_log_analytic(res, z, order != 0, prec);
|
|
|
|
return 0;
|
|
}
|
|
|
|
int
|
|
f_pow_pi(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_const_pi(res, prec);
|
|
acb_pow_analytic(res, z, res, order != 0, prec);
|
|
|
|
return 0;
|
|
}
|
|
|
|
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_real_sqrtpos(res, res, order != 0, 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_real_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) = z sin(1/z), assume on real interval */
|
|
int
|
|
f_essing2(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
|
|
{
|
|
if (order > 1)
|
|
flint_abort(); /* Would be needed for Taylor method. */
|
|
|
|
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) = sech(10(x-0.2))^2 + sech(100(x-0.4))^4 + sech(1000(x-0.6))^6 */
|
|
int
|
|
f_spike(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;
|
|
}
|
|
|
|
/* more tests for the individual acb_real_* functions */
|
|
|
|
/* abs(sin(x))*cos(1+x) */
|
|
int
|
|
f_abs(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
|
|
{
|
|
acb_t s, c;
|
|
acb_init(s);
|
|
acb_init(c);
|
|
acb_sin(s, z, prec);
|
|
acb_add_ui(c, z, 1, prec);
|
|
acb_cos(c, c, prec);
|
|
acb_real_abs(s, s, order != 0, prec);
|
|
acb_mul(res, s, c, prec);
|
|
acb_clear(s);
|
|
acb_clear(c);
|
|
return 0;
|
|
}
|
|
|
|
/* sign(sin(x))*cos(1+x) */
|
|
int
|
|
f_sgn(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
|
|
{
|
|
acb_t s, c;
|
|
acb_init(s);
|
|
acb_init(c);
|
|
acb_sin(s, z, prec);
|
|
acb_add_ui(c, z, 1, prec);
|
|
acb_cos(c, c, prec);
|
|
acb_real_sgn(s, s, order != 0, prec);
|
|
acb_mul(res, s, c, prec);
|
|
acb_clear(s);
|
|
acb_clear(c);
|
|
return 0;
|
|
}
|
|
|
|
/* heaviside(sin(x))*cos(1+x) */
|
|
int
|
|
f_heaviside(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
|
|
{
|
|
acb_t s, c;
|
|
acb_init(s);
|
|
acb_init(c);
|
|
acb_sin(s, z, prec);
|
|
acb_add_ui(c, z, 1, prec);
|
|
acb_cos(c, c, prec);
|
|
acb_real_heaviside(s, s, order != 0, prec);
|
|
acb_mul(res, s, c, prec);
|
|
acb_clear(s);
|
|
acb_clear(c);
|
|
return 0;
|
|
}
|
|
|
|
/* floor(x-5) cos(1+x) */
|
|
int
|
|
f_floor(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
|
|
{
|
|
acb_t c;
|
|
acb_init(c);
|
|
acb_add_ui(c, z, 1, prec);
|
|
acb_cos(c, c, prec);
|
|
|
|
acb_sub_ui(res, z, 5, prec);
|
|
acb_real_floor(res, res, order != 0, prec);
|
|
|
|
acb_mul(res, res, c, prec);
|
|
acb_clear(c);
|
|
return 0;
|
|
}
|
|
|
|
/* ceil(x-5) cos(1+x) */
|
|
int
|
|
f_ceil(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
|
|
{
|
|
acb_t c;
|
|
acb_init(c);
|
|
acb_add_ui(c, z, 1, prec);
|
|
acb_cos(c, c, prec);
|
|
acb_sub_ui(res, z, 5, prec);
|
|
acb_real_ceil(res, res, order != 0, prec);
|
|
acb_mul(res, res, c, prec);
|
|
acb_clear(c);
|
|
return 0;
|
|
}
|
|
|
|
/* max(sin(x),cos(x)) */
|
|
int
|
|
f_max(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
|
|
{
|
|
acb_t s, c;
|
|
acb_init(s);
|
|
acb_init(c);
|
|
acb_sin_cos(s, c, z, prec);
|
|
acb_real_max(res, s, c, order != 0, prec);
|
|
acb_clear(s);
|
|
acb_clear(c);
|
|
return 0;
|
|
}
|
|
|
|
/* min(sin(x),cos(x)) */
|
|
int
|
|
f_min(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
|
|
{
|
|
acb_t s, c;
|
|
acb_init(s);
|
|
acb_init(c);
|
|
acb_sin_cos(s, c, z, prec);
|
|
acb_real_min(res, s, c, order != 0, prec);
|
|
acb_clear(s);
|
|
acb_clear(c);
|
|
return 0;
|
|
}
|
|
|
|
|
|
int main()
|
|
{
|
|
slong iter;
|
|
flint_rand_t state;
|
|
|
|
flint_printf("integrate....");
|
|
fflush(stdout);
|
|
|
|
flint_randinit(state);
|
|
|
|
for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++)
|
|
{
|
|
acb_t ans, res, a, b, t;
|
|
slong goal, prec, abs_goal;
|
|
mag_t tol;
|
|
acb_calc_integrate_opt_t opt;
|
|
int integral;
|
|
|
|
acb_init(ans);
|
|
acb_init(res);
|
|
acb_init(a);
|
|
acb_init(b);
|
|
acb_init(t);
|
|
mag_init(tol);
|
|
acb_calc_integrate_opt_init(opt);
|
|
|
|
goal = 2 + n_randint(state, 300);
|
|
prec = 2 + n_randint(state, 300);
|
|
abs_goal = n_randint(state, 300);
|
|
mag_set_ui_2exp_si(tol, n_randint(state, 2), -abs_goal);
|
|
|
|
opt->eval_limit = n_randint(state, 10000);
|
|
|
|
if (n_randint(state, 2))
|
|
opt->depth_limit = n_randint(state, 100);
|
|
|
|
if (n_randint(state, 2))
|
|
opt->deg_limit = n_randint(state, 100);
|
|
|
|
opt->use_heap = n_randint(state, 2);
|
|
|
|
integral = n_randint(state, 9);
|
|
|
|
if (integral == 0)
|
|
{
|
|
acb_randtest(a, state, 1 + n_randint(state, 200), 2);
|
|
acb_randtest(b, state, 1 + n_randint(state, 200), 2);
|
|
acb_cos(ans, a, prec);
|
|
acb_cos(res, b, prec);
|
|
acb_sub(ans, ans, res, prec);
|
|
acb_calc_integrate(res, f_sin, NULL, a, b, goal, tol, opt, prec);
|
|
}
|
|
else if (integral == 1)
|
|
{
|
|
acb_randtest(a, state, 1 + n_randint(state, 200), 2);
|
|
acb_randtest(b, state, 1 + n_randint(state, 200), 2);
|
|
acb_pow_ui(ans, a, 4, prec);
|
|
acb_pow_ui(res, b, 4, prec);
|
|
acb_sub(ans, res, ans, prec);
|
|
acb_mul_2exp_si(ans, ans, -2);
|
|
acb_calc_integrate(res, f_cube, NULL, a, b, goal, tol, opt, prec);
|
|
}
|
|
else if (integral == 2)
|
|
{
|
|
acb_zero(a);
|
|
acb_one(b);
|
|
acb_const_pi(ans, prec);
|
|
acb_mul_2exp_si(ans, ans, -2);
|
|
if (n_randint(state, 2))
|
|
{
|
|
acb_swap(a, b);
|
|
acb_neg(ans, ans);
|
|
}
|
|
acb_calc_integrate(res, f_circle, NULL, a, b, goal, tol, opt, prec);
|
|
}
|
|
else if (integral == 3)
|
|
{
|
|
acb_zero(a);
|
|
acb_one(b);
|
|
acb_zero(ans);
|
|
arb_set_str(acb_realref(ans), "0.3785300171241613098817352756283519095343133 +/- 1e-40", prec);
|
|
if (n_randint(state, 2))
|
|
{
|
|
acb_swap(a, b);
|
|
acb_neg(ans, ans);
|
|
}
|
|
acb_calc_integrate(res, f_essing2, NULL, a, b, goal, tol, opt, prec);
|
|
}
|
|
else if (integral == 4)
|
|
{
|
|
acb_zero(a);
|
|
acb_one(b);
|
|
acb_zero(ans);
|
|
arb_set_str(acb_realref(ans), "11.147310550057139733915902084255301415775813549800589418261584268232061665808482234384871404010464 +/- 2.28e-97", prec);
|
|
if (n_randint(state, 2))
|
|
{
|
|
acb_swap(a, b);
|
|
acb_neg(ans, ans);
|
|
}
|
|
acb_calc_integrate(res, f_helfgott, NULL, a, b, goal, tol, opt, prec);
|
|
}
|
|
else if (integral == 5)
|
|
{
|
|
acb_zero(a);
|
|
acb_set_ui(b, 8);
|
|
acb_zero(ans);
|
|
arb_set_str(acb_realref(ans), "0.34740017265724780787951215911989312465745625486618018388549271361674821398878532052968510434660 +/- 5.97e-96", prec);
|
|
if (n_randint(state, 2))
|
|
{
|
|
acb_swap(a, b);
|
|
acb_neg(ans, ans);
|
|
}
|
|
acb_calc_integrate(res, f_rump, NULL, a, b, goal, tol, opt, prec);
|
|
}
|
|
else if (integral == 6)
|
|
{
|
|
acb_zero(a);
|
|
acb_one(b);
|
|
acb_zero(ans);
|
|
arb_set_str(acb_realref(ans), "0.21080273550054927737564325570572915436090918643678119034785050587872061312814550020505868926155764 +/- 3.72e-99", prec);
|
|
if (n_randint(state, 2))
|
|
{
|
|
acb_swap(a, b);
|
|
acb_neg(ans, ans);
|
|
}
|
|
acb_calc_integrate(res, f_spike, NULL, a, b, goal, tol, opt, prec);
|
|
}
|
|
else if (integral == 7)
|
|
{
|
|
acb_randtest_special(a, state, 1 + n_randint(state, 200), 2);
|
|
acb_randtest_special(b, state, 1 + n_randint(state, 200), 2);
|
|
acb_zero(ans);
|
|
acb_calc_integrate(res, f_zero, NULL, a, b, goal, tol, opt, prec);
|
|
}
|
|
else if (integral == 8)
|
|
{
|
|
acb_randtest_special(a, state, 1 + n_randint(state, 200), 2);
|
|
acb_randtest_special(b, state, 1 + n_randint(state, 200), 2);
|
|
acb_indeterminate(ans);
|
|
acb_calc_integrate(res, f_indet, NULL, a, b, goal, tol, opt, prec);
|
|
}
|
|
|
|
if (!acb_overlaps(res, ans))
|
|
{
|
|
flint_printf("FAIL! (iter = %wd)\n", iter);
|
|
flint_printf("integral = %d, prec = %wd, goal = %wd\n", integral, prec, goal);
|
|
flint_printf("a = "); acb_printn(a, 150, 0); flint_printf("\n");
|
|
flint_printf("b = "); acb_printn(b, 150, 0); flint_printf("\n");
|
|
flint_printf("res = "); acb_printn(res, 150, 0); flint_printf("\n\n");
|
|
flint_printf("ans = "); acb_printn(ans, 150, 0); flint_printf("\n\n");
|
|
flint_abort();
|
|
}
|
|
|
|
acb_clear(ans);
|
|
acb_clear(res);
|
|
acb_clear(a);
|
|
acb_clear(b);
|
|
acb_clear(t);
|
|
mag_clear(tol);
|
|
}
|
|
|
|
/* more tests for the individual real extensions and branched functions */
|
|
{
|
|
acb_t a, b, z, w;
|
|
slong prec;
|
|
mag_t tol;
|
|
|
|
acb_init(a);
|
|
acb_init(b);
|
|
acb_init(z);
|
|
acb_init(w);
|
|
mag_init(tol);
|
|
|
|
acb_zero(a);
|
|
acb_set_ui(b, 10);
|
|
|
|
prec = 53;
|
|
mag_set_ui_2exp_si(tol, 1, -prec);
|
|
|
|
acb_calc_integrate(z, f_abs, NULL, a, b, prec, tol, NULL, prec);
|
|
arb_set_str(acb_realref(w), "-1.3517710956465318592 +/- 1e-17", prec);
|
|
|
|
if (!acb_overlaps(z, w) || acb_rel_accuracy_bits(z) < prec - 10)
|
|
{
|
|
flint_printf("FAIL (abs)\n");
|
|
flint_printf("z = "); acb_printn(z, 20, 0); flint_printf("\n");
|
|
flint_printf("w = "); acb_printn(w, 20, 0); flint_printf("\n");
|
|
flint_abort();
|
|
}
|
|
|
|
acb_calc_integrate(z, f_sgn, NULL, a, b, prec, tol, NULL, prec);
|
|
arb_set_str(acb_realref(w), "-4.8903066871045720895 +/- 1e-17", prec);
|
|
|
|
if (!acb_overlaps(z, w) || acb_rel_accuracy_bits(z) < prec - 10)
|
|
{
|
|
flint_printf("FAIL (sgn)\n");
|
|
flint_printf("z = "); acb_printn(z, 20, 0); flint_printf("\n");
|
|
flint_printf("w = "); acb_printn(w, 20, 0); flint_printf("\n");
|
|
flint_abort();
|
|
}
|
|
|
|
acb_calc_integrate(z, f_heaviside, NULL, a, b, prec, tol, NULL, prec);
|
|
arb_set_str(acb_realref(w), "-3.3658839392315860266 +/- 1e-17", prec);
|
|
|
|
if (!acb_overlaps(z, w) || acb_rel_accuracy_bits(z) < prec - 10)
|
|
{
|
|
flint_printf("FAIL (heaviside)\n");
|
|
flint_printf("z = "); acb_printn(z, 20, 0); flint_printf("\n");
|
|
flint_printf("w = "); acb_printn(w, 20, 0); flint_printf("\n");
|
|
flint_abort();
|
|
}
|
|
|
|
acb_calc_integrate(z, f_floor, NULL, a, b, prec, tol, NULL, prec);
|
|
arb_set_str(acb_realref(w), "-0.36232328857344524392 +/- 1e-17", prec);
|
|
|
|
if (!acb_overlaps(z, w) || acb_rel_accuracy_bits(z) < prec - 16)
|
|
{
|
|
flint_printf("FAIL (floor)\n");
|
|
flint_printf("z = "); acb_printn(z, 20, 0); flint_printf("\n");
|
|
flint_printf("w = "); acb_printn(w, 20, 0); flint_printf("\n");
|
|
flint_abort();
|
|
}
|
|
|
|
acb_calc_integrate(z, f_ceil, NULL, a, b, prec, tol, NULL, prec);
|
|
arb_set_str(acb_realref(w), "-2.2037844799320452076 +/- 1e-17", prec);
|
|
|
|
if (!acb_overlaps(z, w) || acb_rel_accuracy_bits(z) < prec - 16)
|
|
{
|
|
flint_printf("FAIL (ceil)\n");
|
|
flint_printf("z = "); acb_printn(z, 20, 0); flint_printf("\n");
|
|
flint_printf("w = "); acb_printn(w, 20, 0); flint_printf("\n");
|
|
flint_abort();
|
|
}
|
|
|
|
acb_calc_integrate(z, f_max, NULL, a, b, prec, tol, NULL, prec);
|
|
arb_set_str(acb_realref(w), "5.0817122161957375987 +/- 1e-17", prec);
|
|
|
|
if (!acb_overlaps(z, w) || acb_rel_accuracy_bits(z) < prec - 10)
|
|
{
|
|
flint_printf("FAIL (max)\n");
|
|
flint_printf("z = "); acb_printn(z, 20, 0); flint_printf("\n");
|
|
flint_printf("w = "); acb_printn(w, 20, 0); flint_printf("\n");
|
|
flint_abort();
|
|
}
|
|
|
|
acb_calc_integrate(z, f_min, NULL, a, b, prec, tol, NULL, prec);
|
|
arb_set_str(acb_realref(w), "-3.7866617980086549598 +/- 1e-17", prec);
|
|
|
|
if (!acb_overlaps(z, w) || acb_rel_accuracy_bits(z) < prec - 10)
|
|
{
|
|
flint_printf("FAIL (min)\n");
|
|
flint_printf("z = "); acb_printn(z, 20, 0); flint_printf("\n");
|
|
flint_printf("w = "); acb_printn(w, 20, 0); flint_printf("\n");
|
|
flint_abort();
|
|
}
|
|
|
|
prec = 64;
|
|
mag_set_ui_2exp_si(tol, 1, -prec);
|
|
|
|
acb_set_ui(a, 10);
|
|
acb_set_ui(b, 11);
|
|
acb_calc_integrate(z, f_sqrt, NULL, a, b, prec, tol, NULL, prec);
|
|
arb_set_str(acb_realref(w), "3.2400640614837366802 +/- 1.65e-20", prec);
|
|
|
|
if (!acb_overlaps(z, w) || acb_rel_accuracy_bits(z) < prec - 10)
|
|
{
|
|
flint_printf("FAIL (sqrt)\n");
|
|
flint_printf("z = "); acb_printn(z, 20, 0); flint_printf("\n");
|
|
flint_printf("w = "); acb_printn(w, 20, 0); flint_printf("\n");
|
|
flint_abort();
|
|
}
|
|
|
|
acb_calc_integrate(z, f_rsqrt, NULL, a, b, prec, tol, NULL, prec);
|
|
arb_set_str(acb_realref(w), "0.30869426037404103423 +/- 2.08e-21", prec);
|
|
|
|
if (!acb_overlaps(z, w) || acb_rel_accuracy_bits(z) < prec - 10)
|
|
{
|
|
flint_printf("FAIL (rsqrt)\n");
|
|
flint_printf("z = "); acb_printn(z, 20, 0); flint_printf("\n");
|
|
flint_printf("w = "); acb_printn(w, 20, 0); flint_printf("\n");
|
|
flint_abort();
|
|
}
|
|
|
|
acb_calc_integrate(z, f_log, NULL, a, b, prec, tol, NULL, prec);
|
|
arb_set_str(acb_realref(w), "2.3509970708416191445 +/- 1.47e-21", prec);
|
|
|
|
if (!acb_overlaps(z, w) || acb_rel_accuracy_bits(z) < prec - 10)
|
|
{
|
|
flint_printf("FAIL (log)\n");
|
|
flint_printf("z = "); acb_printn(z, 20, 0); flint_printf("\n");
|
|
flint_printf("w = "); acb_printn(w, 20, 0); flint_printf("\n");
|
|
flint_abort();
|
|
}
|
|
|
|
acb_calc_integrate(z, f_pow_pi, NULL, a, b, prec, tol, NULL, prec);
|
|
arb_set_str(acb_realref(w), "1619.0628341842463436 +/- 2.79e-17", prec);
|
|
|
|
if (!acb_overlaps(z, w) || acb_rel_accuracy_bits(z) < prec - 10)
|
|
{
|
|
flint_printf("FAIL (pow_pi)\n");
|
|
flint_printf("z = "); acb_printn(z, 20, 0); flint_printf("\n");
|
|
flint_printf("w = "); acb_printn(w, 20, 0); flint_printf("\n");
|
|
flint_abort();
|
|
}
|
|
|
|
acb_set_d_d(a, -10, -1);
|
|
acb_set_d_d(b, -10, +1);
|
|
acb_zero(w);
|
|
|
|
acb_calc_integrate(z, f_sqrt, NULL, a, b, prec, tol, NULL, prec);
|
|
arb_set_str(acb_imagref(w), "0.15801534879296271761 +/- 4.91e-21", prec);
|
|
|
|
if (!acb_overlaps(z, w) || acb_rel_accuracy_bits(z) < prec - 25)
|
|
{
|
|
flint_printf("FAIL (sqrt 2)\n");
|
|
flint_printf("z = "); acb_printn(z, 20, 0); flint_printf("\n");
|
|
flint_printf("w = "); acb_printn(w, 20, 0); flint_printf("\n");
|
|
flint_abort();
|
|
}
|
|
|
|
acb_calc_integrate(z, f_rsqrt, NULL, a, b, prec, tol, NULL, prec);
|
|
arb_set_str(acb_imagref(w), "0.015762235473606564294 +/- 2.80e-22", prec);
|
|
|
|
if (!acb_overlaps(z, w) || acb_rel_accuracy_bits(z) < prec - 25)
|
|
{
|
|
flint_printf("FAIL (rsqrt 2)\n");
|
|
flint_printf("z = "); acb_printn(z, 20, 0); flint_printf("\n");
|
|
flint_printf("w = "); acb_printn(w, 20, 0); flint_printf("\n");
|
|
flint_abort();
|
|
}
|
|
|
|
acb_calc_integrate(z, f_log, NULL, a, b, prec, tol, NULL, prec);
|
|
arb_set_str(acb_imagref(w), "4.6084935666644999985 +/- 4.69e-20", prec);
|
|
|
|
if (!acb_overlaps(z, w) || acb_rel_accuracy_bits(z) < prec - 25)
|
|
{
|
|
flint_printf("FAIL (log 2)\n");
|
|
flint_printf("z = "); acb_printn(z, 20, 0); flint_printf("\n");
|
|
flint_printf("w = "); acb_printn(w, 20, 0); flint_printf("\n");
|
|
flint_abort();
|
|
}
|
|
|
|
acb_calc_integrate(z, f_pow_pi, NULL, a, b, prec, tol, NULL, prec);
|
|
arb_set_str(acb_imagref(w), "-2660.1245860252382407 +/- 8.64e-18", prec);
|
|
|
|
if (!acb_overlaps(z, w) || acb_rel_accuracy_bits(z) < prec - 25)
|
|
{
|
|
flint_printf("FAIL (pow_pi 2)\n");
|
|
flint_printf("z = "); acb_printn(z, 20, 0); flint_printf("\n");
|
|
flint_printf("w = "); acb_printn(w, 20, 0); flint_printf("\n");
|
|
flint_abort();
|
|
}
|
|
|
|
prec = 100;
|
|
mag_set_ui_2exp_si(tol, 1, -prec);
|
|
|
|
acb_set_d_d(a, -100, -1);
|
|
acb_set_d_d(b, -100, +2);
|
|
acb_zero(w);
|
|
acb_calc_integrate(z, f_log, NULL, a, b, prec, tol, NULL, prec);
|
|
arb_set_str(acb_realref(w), "-3.12659390337983876281336380373 +/- 3.97e-30", prec);
|
|
arb_set_str(acb_imagref(w), "13.8156605414673448203655975232 +/- 2.39e-29", prec);
|
|
|
|
if (!acb_overlaps(z, w) || acb_rel_accuracy_bits(z) < prec - 25)
|
|
{
|
|
flint_printf("FAIL (log 3)\n");
|
|
flint_printf("z = "); acb_printn(z, 20, 0); flint_printf("\n");
|
|
flint_printf("w = "); acb_printn(w, 20, 0); flint_printf("\n");
|
|
flint_abort();
|
|
}
|
|
|
|
acb_clear(a);
|
|
acb_clear(b);
|
|
acb_clear(z);
|
|
acb_clear(w);
|
|
mag_clear(tol);
|
|
}
|
|
|
|
flint_randclear(state);
|
|
flint_cleanup();
|
|
flint_printf("PASS\n");
|
|
return EXIT_SUCCESS;
|
|
}
|
|
|