diff --git a/examples/integrals.c b/examples/integrals.c index b6458a63ca..d93d3ae580 100644 --- a/examples/integrals.c +++ b/examples/integrals.c @@ -732,11 +732,22 @@ scaled_bessel_select_N(arb_t N, ulong k, slong prec) arb_mul_2exp_si(N, N, e); } +/* f(z) = Si(z) */ +int +f_si(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_hypgeom_si(res, z, prec); + return 0; +} + /* ------------------------------------------------------------------------- */ /* Main test program */ /* ------------------------------------------------------------------------- */ -#define NUM_INTEGRALS 37 +#define NUM_INTEGRALS 38 const char * descr[NUM_INTEGRALS] = { @@ -777,6 +788,7 @@ const char * descr[NUM_INTEGRALS] = "int_0^{inf} exp(-x) I_0(x/15)^{15} dx (using domain truncation)", "int_{-1-i}^{-1+i} 1/sqrt(x) dx", "int_0^{inf} 1/gamma(x) dx (using domain truncation)", + "int_0^{1000} Si(x) dx", }; int main(int argc, char *argv[]) @@ -1273,6 +1285,13 @@ int main(int argc, char *argv[]) arb_add_error_2exp_si(acb_realref(s), -goal); break; + case 37: + acb_zero(a); + acb_set_ui(b, 1000); + acb_calc_integrate(s, f_si, NULL, a, b, goal, tol, options, prec); + break; + + default: abort(); } diff --git a/src/acb_hypgeom/chi.c b/src/acb_hypgeom/chi.c index 32a79e63b4..1a1b86f852 100644 --- a/src/acb_hypgeom/chi.c +++ b/src/acb_hypgeom/chi.c @@ -206,6 +206,7 @@ acb_hypgeom_chi(acb_t res, const acb_t z, slong prec) wp = prec + FLINT_MAX(0, cancellation); wp = wp * 1.001 + 5; acb_hypgeom_chi_2f3(res, m, wp); + acb_set_round(res, res, prec); } if (pure_real) diff --git a/src/acb_hypgeom/ci.c b/src/acb_hypgeom/ci.c index b5556dc1e0..2aeac193e2 100644 --- a/src/acb_hypgeom/ci.c +++ b/src/acb_hypgeom/ci.c @@ -271,6 +271,7 @@ acb_hypgeom_ci(acb_t res, const acb_t z, slong prec) wp = prec + FLINT_MAX(0, cancellation); wp = wp * 1.001 + 5; acb_hypgeom_ci_2f3(res, m, wp); + acb_set_round(res, res, prec); } acb_add_error_mag(res, err); diff --git a/src/acb_hypgeom/ei.c b/src/acb_hypgeom/ei.c index ccfcb8078c..fc1abb8416 100644 --- a/src/acb_hypgeom/ei.c +++ b/src/acb_hypgeom/ei.c @@ -236,6 +236,7 @@ acb_hypgeom_ei(acb_t res, const acb_t z, slong prec) wp = prec + FLINT_MAX(0, cancellation); wp = wp * 1.001 + 5; acb_hypgeom_ei_2f2(res, m, wp); + acb_set_round(res, res, prec); } if (pure_real) diff --git a/src/acb_hypgeom/si.c b/src/acb_hypgeom/si.c index 0a66e56dfd..0f2db1bdfd 100644 --- a/src/acb_hypgeom/si.c +++ b/src/acb_hypgeom/si.c @@ -215,6 +215,7 @@ acb_hypgeom_si(acb_t res, const acb_t z, slong prec) wp = prec + FLINT_MAX(0, cancellation); wp = wp * 1.001 + 5; acb_hypgeom_si_1f2(res, m, wp); + acb_set_round(res, res, prec); } if (pure_imag) diff --git a/src/python/flint_ctypes.py b/src/python/flint_ctypes.py index f12975d88c..2cc0ac1698 100644 --- a/src/python/flint_ctypes.py +++ b/src/python/flint_ctypes.py @@ -1913,7 +1913,7 @@ def exp_integral(ctx, x, y): def exp_integral_ei(ctx, x): """ >>> RR.exp_integral_ei(1) - [1.89511781635594 +/- 5.11e-15] + [1.895117816355937 +/- 6.91e-16] """ return ctx._unary_op(x, libgr.gr_exp_integral_ei, "exp_integral_ei($x)") @@ -1921,6 +1921,10 @@ def sin_integral(ctx, x): """ >>> RR.sin_integral(1) [0.946083070367183 +/- 1.35e-16] + >>> CC.sin_integral(CC("(5 +/- 1e-10)*I")) + [20.09321183 +/- 5.79e-9]*I + >>> CC.sin_integral(CC("10 + (1 +/- 1e-10)*I")) + ([1.7002629761 +/- 6.33e-11] + [-0.0667638998 +/- 2.17e-11]*I) """ return ctx._unary_op(x, libgr.gr_sin_integral, "sin_integral($x)") @@ -1934,21 +1938,21 @@ def cos_integral(ctx, x): def sinh_integral(ctx, x): """ >>> RR.sinh_integral(1) - [1.05725087537573 +/- 2.77e-15] + [1.057250875375728 +/- 6.29e-16] """ return ctx._unary_op(x, libgr.gr_sinh_integral, "sinh_integral($x)") def cosh_integral(ctx, x): """ >>> RR.cosh_integral(1) - [0.837866940980208 +/- 4.78e-16] + [0.837866940980208 +/- 3.15e-16] """ return ctx._unary_op(x, libgr.gr_cosh_integral, "cosh_integral($x)") def log_integral(ctx, x, offset=False): """ >>> RR.log_integral(2) - [1.04516378011749 +/- 4.01e-15] + [1.04516378011749 +/- 3.31e-15] >>> RR.log_integral(2, offset=True) 0 """