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 cdec409c6a..1a1b86f852 100644
--- a/src/acb_hypgeom/chi.c
+++ b/src/acb_hypgeom/chi.c
@@ -9,6 +9,7 @@
     (at your option) any later version.  See <https://www.gnu.org/licenses/>.
 */
 
+#include <math.h>
 #include "acb.h"
 #include "acb_hypgeom.h"
 
@@ -101,51 +102,119 @@ acb_hypgeom_chi_asymp(acb_t res, const acb_t z, slong prec)
     acb_clear(one);
 }
 
+/* defined in ci.c */
 void
-acb_hypgeom_chi_2f3(acb_t res, const acb_t z, slong prec)
-{
-    acb_t a, t, u;
-    acb_struct b[3];
-
-    acb_init(a);
-    acb_init(b);
-    acb_init(b + 1);
-    acb_init(b + 2);
-    acb_init(t);
-    acb_init(u);
-
-    acb_one(a);
-    acb_set_ui(b, 2);
-    acb_set(b + 1, b);
-    acb_set_ui(b + 2, 3);
-    acb_mul_2exp_si(b + 2, b + 2, -1);
-
-    acb_mul(t, z, z, prec);
-    acb_mul_2exp_si(t, t, -2);
-    acb_hypgeom_pfq_direct(u, a, 1, b, 3, t, -1, prec);
-    acb_mul(u, u, t, prec);
+acb_hypgeom_chi_2f3(acb_t res, const acb_t z, slong prec);
 
-    acb_log(t, z, prec);
-    acb_add(u, u, t, prec);
+/* Bound propagated error |chi(z)-chi(mid(z))| assuming that we
+   are off the branch cut (not checked here).
+   Uses |chi'(z)| = |cosh(z)/z| <= cosh(|re(z)|)/|z|.
+*/
+static void
+acb_hypgeom_chi_prop_error(mag_t err, const acb_t z)
+{
+    if (acb_is_exact(z))
+    {
+        mag_zero(err);
+    }
+    else
+    {
+        mag_t t, u;
 
-    arb_const_euler(acb_realref(t), prec);
-    arb_add(acb_realref(u), acb_realref(u), acb_realref(t), prec);
+        mag_init(t);
+        mag_init(u);
 
-    acb_swap(res, u);
+        arb_get_mag(t, acb_realref(z));
+        mag_cosh(t, t);
+        acb_get_mag_lower(u, z);
+        mag_div(t, t, u);
+        mag_hypot(u, arb_radref(acb_realref(z)), arb_radref(acb_imagref(z)));
+        mag_mul(err, t, u);
 
-    acb_clear(a);
-    acb_clear(b);
-    acb_clear(b + 1);
-    acb_clear(b + 2);
-    acb_clear(t);
-    acb_clear(u);
+        mag_clear(t);
+        mag_clear(u);
+    }
 }
 
 void
 acb_hypgeom_chi(acb_t res, const acb_t z, slong prec)
 {
-    if (acb_hypgeom_u_use_asymp(z, prec))
-        acb_hypgeom_chi_asymp(res, z, prec);
+    if (!acb_is_finite(z) || acb_contains_zero(z))
+    {
+        acb_indeterminate(res);
+    }
+/* todo -- good real code
+    else if (acb_is_real(z))
+    {
+    }
+*/
+    else if (arb_contains_zero(acb_imagref(z)) &&
+                !arb_is_nonnegative(acb_imagref(z)) &&
+                !arb_is_positive(acb_realref(z)))
+    {
+        /* We straddle the branch cut; do something generic */
+        if (acb_hypgeom_u_use_asymp(z, prec))
+            acb_hypgeom_chi_asymp(res, z, prec);
+        else
+            acb_hypgeom_chi_2f3(res, z, prec);
+    }
     else
-        acb_hypgeom_chi_2f3(res, z, prec);
+    {
+        double asymp_accuracy, a, b, absz, cancellation, rlog2 = 1.4426950408889634;
+        slong wp;
+        int use_asymp;
+        acb_t m;
+        mag_t err;
+        /* since we don't have a special case for real z */
+        int pure_real = arb_is_zero(acb_imagref(z));
+
+        a = arf_get_d(arb_midref(acb_realref(z)), ARF_RND_DOWN);
+        b = arf_get_d(arb_midref(acb_imagref(z)), ARF_RND_DOWN);
+        a = fabs(a);
+        b = fabs(b);
+        absz = sqrt(a * a + b * b);
+
+        if (a <= 1.0 && b <= 1.0)
+        {
+            use_asymp = 0;
+        }
+        else if (a > prec || b > prec)
+        {
+            use_asymp = 1;
+        }
+        else
+        {
+            asymp_accuracy = absz * rlog2 * 0.999 - 5;
+            use_asymp = asymp_accuracy > prec;
+        }
+
+        acb_init(m);
+        mag_init(err);
+
+        /* assumes that we already handled the branch cut */
+        acb_hypgeom_chi_prop_error(err, z);
+        acb_get_mid(m, z);
+
+        if (use_asymp)
+        {
+            acb_hypgeom_chi_asymp(res, m, prec);
+        }
+        else
+        {
+            /* terms grow to ~ exp(|z|), sum is ~ exp(|re(z)|) */
+            cancellation = (absz - a) * rlog2;
+            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)
+            arb_add_error_mag(acb_realref(res), err);
+        else
+            acb_add_error_mag(res, err);
+
+        acb_clear(m);
+        mag_clear(err);
+    }
 }
diff --git a/src/acb_hypgeom/ci.c b/src/acb_hypgeom/ci.c
index 7c9ce2351b..2aeac193e2 100644
--- a/src/acb_hypgeom/ci.c
+++ b/src/acb_hypgeom/ci.c
@@ -1,5 +1,5 @@
 /*
-    Copyright (C) 2015 Fredrik Johansson
+    Copyright (C) 2015, 2024 Fredrik Johansson
 
     This file is part of FLINT.
 
@@ -9,6 +9,7 @@
     (at your option) any later version.  See <https://www.gnu.org/licenses/>.
 */
 
+#include <math.h>
 #include "acb.h"
 #include "arb_hypgeom.h"
 #include "acb_hypgeom.h"
@@ -109,8 +110,8 @@ acb_hypgeom_ci_asymp(acb_t res, const acb_t z, slong prec)
     acb_clear(one);
 }
 
-void
-acb_hypgeom_ci_2f3(acb_t res, const acb_t z, slong prec)
+static void
+_acb_hypgeom_ci_2f3(acb_t res, const acb_t z, int hyperbolic, slong prec)
 {
     acb_t a, t, u;
     acb_struct b[3];
@@ -130,7 +131,8 @@ acb_hypgeom_ci_2f3(acb_t res, const acb_t z, slong prec)
 
     acb_mul(t, z, z, prec);
     acb_mul_2exp_si(t, t, -2);
-    acb_neg(t, t);
+    if (!hyperbolic)
+        acb_neg(t, t);
     acb_hypgeom_pfq_direct(u, a, 1, b, 3, t, -1, prec);
     acb_mul(u, u, t, prec);
 
@@ -150,31 +152,131 @@ acb_hypgeom_ci_2f3(acb_t res, const acb_t z, slong prec)
     acb_clear(u);
 }
 
+void
+acb_hypgeom_ci_2f3(acb_t res, const acb_t z, slong prec)
+{
+    _acb_hypgeom_ci_2f3(res, z, 0, prec);
+}
+
+void
+acb_hypgeom_chi_2f3(acb_t res, const acb_t z, slong prec)
+{
+    _acb_hypgeom_ci_2f3(res, z, 1, prec);
+}
+
+/* Bound propagated error |ci(z)-ci(mid(z))| assuming that we
+   are off the branch cut (not checked here).
+   Uses |ci'(z)| = |cos(z)/z| <= cosh(|im(z)|)/|z|.
+*/
+static void
+acb_hypgeom_ci_prop_error(mag_t err, const acb_t z)
+{
+    if (acb_is_exact(z))
+    {
+        mag_zero(err);
+    }
+    else
+    {
+        mag_t t, u;
+
+        mag_init(t);
+        mag_init(u);
+
+        arb_get_mag(t, acb_imagref(z));
+        mag_cosh(t, t);
+        acb_get_mag_lower(u, z);
+        mag_div(t, t, u);
+        mag_hypot(u, arb_radref(acb_realref(z)), arb_radref(acb_imagref(z)));
+        mag_mul(err, t, u);
+
+        mag_clear(t);
+        mag_clear(u);
+    }
+}
+
 void
 acb_hypgeom_ci(acb_t res, const acb_t z, slong prec)
 {
-    if (acb_is_real(z) && arb_is_finite(acb_realref(z)))
+    if (!acb_is_finite(z) || acb_contains_zero(z))
+    {
+        acb_indeterminate(res);
+    }
+    else if (acb_is_real(z))
     {
         if (arb_is_positive(acb_realref(z)))
         {
             arb_hypgeom_ci(acb_realref(res), acb_realref(z), prec);
             arb_zero(acb_imagref(res));
         }
-        else if (arb_is_negative(acb_realref(z)))
+        else /* arb_is_negative(acb_realref(z)), since we already tested for 0 */
         {
             arb_neg(acb_realref(res), acb_realref(z));
             arb_hypgeom_ci(acb_realref(res), acb_realref(res), prec);
             arb_const_pi(acb_imagref(res), prec);
         }
+    }
+    else if (arb_contains_zero(acb_imagref(z)) &&
+                !arb_is_nonnegative(acb_imagref(z)) &&
+                !arb_is_positive(acb_realref(z)))
+    {
+        /* We straddle the branch cut; do something generic */
+        if (acb_hypgeom_u_use_asymp(z, prec))
+            acb_hypgeom_ci_asymp(res, z, prec);
         else
+            acb_hypgeom_ci_2f3(res, z, prec);
+    }
+    else
+    {
+        double asymp_accuracy, a, b, absz, cancellation, rlog2 = 1.4426950408889634;
+        slong wp;
+        int use_asymp;
+        acb_t m;
+        mag_t err;
+
+        a = arf_get_d(arb_midref(acb_realref(z)), ARF_RND_DOWN);
+        b = arf_get_d(arb_midref(acb_imagref(z)), ARF_RND_DOWN);
+        a = fabs(a);
+        b = fabs(b);
+        absz = sqrt(a * a + b * b);
+
+        if (a <= 1.0 && b <= 1.0)
         {
-            acb_indeterminate(res);
+            use_asymp = 0;
+        }
+        else if (a > prec || b > prec)
+        {
+            use_asymp = 1;
+        }
+        else
+        {
+            asymp_accuracy = absz * rlog2 * 0.999 - 5;
+            use_asymp = asymp_accuracy > prec;
         }
-        return;
-    }
 
-    if (acb_hypgeom_u_use_asymp(z, prec))
-        acb_hypgeom_ci_asymp(res, z, prec);
-    else
-        acb_hypgeom_ci_2f3(res, z, prec);
+        acb_init(m);
+        mag_init(err);
+
+        /* assumes that we already handled the branch cut */
+        acb_hypgeom_ci_prop_error(err, z);
+        acb_get_mid(m, z);
+
+        if (use_asymp)
+        {
+            acb_hypgeom_ci_asymp(res, m, prec);
+        }
+        else
+        {
+            /* terms grow to ~ exp(|z|), sum is ~ exp(|im(z)|) */
+            cancellation = (absz - b) * rlog2;
+            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);
+
+        acb_clear(m);
+        mag_clear(err);
+    }
 }
diff --git a/src/acb_hypgeom/ei.c b/src/acb_hypgeom/ei.c
index e81d7f4ee7..fc1abb8416 100644
--- a/src/acb_hypgeom/ei.c
+++ b/src/acb_hypgeom/ei.c
@@ -9,6 +9,7 @@
     (at your option) any later version.  See <https://www.gnu.org/licenses/>.
 */
 
+#include <math.h>
 #include "acb.h"
 #include "acb_hypgeom.h"
 
@@ -120,11 +121,130 @@ acb_hypgeom_ei_2f2(acb_t res, const acb_t z, slong prec)
     acb_clear(t);
 }
 
+/* Bound propagated error |ei(z)-ei(mid(z))| assuming that we
+   are off the branch cut (not checked here).
+   Uses |ei'(z)| = |exp(z)/z| <= exp(re(z))/|z|.
+*/
+static void
+acb_hypgeom_ei_prop_error(mag_t err, const acb_t z)
+{
+    if (acb_is_exact(z))
+    {
+        mag_zero(err);
+    }
+    else
+    {
+        mag_t t, u;
+        arf_t x;
+
+        mag_init(t);
+        mag_init(u);
+        arf_init(x);
+
+        arb_get_ubound_arf(x, acb_realref(z), MAG_BITS);
+
+        if (arf_sgn(x) >= 0)
+        {
+            arf_get_mag(t, x);
+            mag_exp(t, t);
+        }
+        else
+        {
+            arf_get_mag_lower(t, x);
+            mag_expinv(t, t);
+        }
+
+        acb_get_mag_lower(u, z);
+        mag_div(t, t, u);
+        mag_hypot(u, arb_radref(acb_realref(z)), arb_radref(acb_imagref(z)));
+        mag_mul(err, t, u);
+
+        mag_clear(t);
+        mag_clear(u);
+        arf_clear(x);
+    }
+}
+
 void
 acb_hypgeom_ei(acb_t res, const acb_t z, slong prec)
 {
-    if (acb_hypgeom_u_use_asymp(z, prec))
-        acb_hypgeom_ei_asymp(res, z, prec);
+    if (!acb_is_finite(z) || acb_contains_zero(z))
+    {
+        acb_indeterminate(res);
+    }
+/* todo -- good real code
+    else if (acb_is_real(z))
+    {
+    }
+*/
+    else if (arb_contains_zero(acb_imagref(z)) &&
+                !arb_is_zero(acb_imagref(z)) &&
+                !arb_is_positive(acb_realref(z)))
+    {
+        /* We straddle the branch cut; do something generic */
+        if (acb_hypgeom_u_use_asymp(z, prec))
+            acb_hypgeom_ei_asymp(res, z, prec);
+        else
+            acb_hypgeom_ei_2f2(res, z, prec);
+    }
     else
-        acb_hypgeom_ei_2f2(res, z, prec);
+    {
+        double asymp_accuracy, a, asigned, b, absz, cancellation, rlog2 = 1.4426950408889634;
+        slong wp;
+        int use_asymp;
+        acb_t m;
+        mag_t err;
+        /* since we don't have a special case for real z */
+        int pure_real = arb_is_zero(acb_imagref(z));
+
+        a = arf_get_d(arb_midref(acb_realref(z)), ARF_RND_DOWN);
+        b = arf_get_d(arb_midref(acb_imagref(z)), ARF_RND_DOWN);
+        asigned = a;
+        a = fabs(a);
+        b = fabs(b);
+        absz = sqrt(a * a + b * b);
+
+        if (a <= 1.0 && b <= 1.0)
+        {
+            use_asymp = 0;
+        }
+        else if (a > prec || b > prec)
+        {
+            use_asymp = 1;
+        }
+        else
+        {
+            asymp_accuracy = absz * rlog2 * 0.999 - 5;
+            use_asymp = asymp_accuracy > prec;
+        }
+
+        acb_init(m);
+        mag_init(err);
+
+        /* assumes that we already handled the branch cut */
+        acb_hypgeom_ei_prop_error(err, z);
+        acb_get_mid(m, z);
+
+        if (use_asymp)
+        {
+            acb_hypgeom_ei_asymp(res, m, prec);
+        }
+        else
+        {
+            /* terms grow to ~ exp(|z|), sum is ~ exp(re(z)) */
+            cancellation = (absz - asigned) * rlog2;
+            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)
+            arb_add_error_mag(acb_realref(res), err);
+        else
+            acb_add_error_mag(res, err);
+
+        acb_clear(m);
+        mag_clear(err);
+    }
 }
diff --git a/src/acb_hypgeom/si.c b/src/acb_hypgeom/si.c
index 704f9096b3..0f2db1bdfd 100644
--- a/src/acb_hypgeom/si.c
+++ b/src/acb_hypgeom/si.c
@@ -1,5 +1,5 @@
 /*
-    Copyright (C) 2015 Fredrik Johansson
+    Copyright (C) 2015, 2024 Fredrik Johansson
 
     This file is part of FLINT.
 
@@ -9,6 +9,7 @@
     (at your option) any later version.  See <https://www.gnu.org/licenses/>.
 */
 
+#include <math.h>
 #include "acb.h"
 #include "arb_hypgeom.h"
 #include "acb_hypgeom.h"
@@ -80,6 +81,9 @@ acb_hypgeom_si_asymp(acb_t res, const acb_t z, slong prec)
 
     acb_swap(res, t);
 
+    if (!acb_is_finite(res))
+        acb_indeterminate(res);
+
     acb_clear(t);
     acb_clear(u);
     acb_clear(w);
@@ -121,18 +125,105 @@ acb_hypgeom_si_1f2(acb_t res, const acb_t z, slong prec)
     acb_clear(t);
 }
 
+/* Bound propagated error |si(z)-si(mid(z))| using
+   |si'(z)| = |sin(z)/z| <= cosh(|im(z)|)/max(1, |z|).
+*/
+static void
+acb_hypgeom_si_prop_error(mag_t err, const acb_t z)
+{
+    if (acb_is_exact(z))
+    {
+        mag_zero(err);
+    }
+    else
+    {
+        mag_t t, u;
+
+        mag_init(t);
+        mag_init(u);
+
+        arb_get_mag(t, acb_imagref(z));
+        mag_cosh(t, t);
+        acb_get_mag_lower(u, z);
+        if (mag_cmp_2exp_si(u, 0) < 0)
+            mag_one(u);
+        mag_div(t, t, u);
+        mag_hypot(u, arb_radref(acb_realref(z)), arb_radref(acb_imagref(z)));
+        mag_mul(err, t, u);
+
+        mag_clear(t);
+        mag_clear(u);
+    }
+}
+
 void
 acb_hypgeom_si(acb_t res, const acb_t z, slong prec)
 {
-    if (acb_is_real(z) && arb_is_finite(acb_realref(z)))
+    if (!acb_is_finite(z))
+    {
+        acb_indeterminate(res);
+    }
+    else if (acb_is_real(z))
     {
         arb_hypgeom_si(acb_realref(res), acb_realref(z), prec);
         arb_zero(acb_imagref(res));
-        return;
     }
-
-    if (acb_hypgeom_u_use_asymp(z, prec))
-        acb_hypgeom_si_asymp(res, z, prec);
     else
-        acb_hypgeom_si_1f2(res, z, prec);
+    {
+        double asymp_accuracy, a, b, absz, cancellation, rlog2 = 1.4426950408889634;
+        slong wp;
+        int use_asymp;
+        acb_t m;
+        mag_t err;
+        int pure_imag = arb_is_zero(acb_realref(z));
+
+        a = arf_get_d(arb_midref(acb_realref(z)), ARF_RND_DOWN);
+        b = arf_get_d(arb_midref(acb_imagref(z)), ARF_RND_DOWN);
+        a = fabs(a);
+        b = fabs(b);
+        absz = sqrt(a * a + b * b);
+
+        if (a <= 1.0 && b <= 1.0)
+        {
+            use_asymp = 0;
+        }
+        else if (a > prec || b > prec)
+        {
+            use_asymp = 1;
+        }
+        else
+        {
+            asymp_accuracy = absz * rlog2 * 0.999 - 5;
+            use_asymp = asymp_accuracy > prec;
+        }
+
+        acb_init(m);
+        mag_init(err);
+
+        /* assumes that we already handled the branch cut */
+        acb_hypgeom_si_prop_error(err, z);
+        acb_get_mid(m, z);
+
+        if (use_asymp)
+        {
+            acb_hypgeom_si_asymp(res, m, prec);
+        }
+        else
+        {
+            /* terms grow to ~ exp(|z|), sum is ~ exp(|im(z)|) */
+            cancellation = (absz - b) * rlog2;
+            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)
+            arb_add_error_mag(acb_imagref(res), err);
+        else
+            acb_add_error_mag(res, err);
+
+        acb_clear(m);
+        mag_clear(err);
+    }
 }
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
         """