diff --git a/acb_hypgeom/ci.c b/acb_hypgeom/ci.c
index a1239ba7..bfa0ddbb 100644
--- a/acb_hypgeom/ci.c
+++ b/acb_hypgeom/ci.c
@@ -9,6 +9,7 @@
(at your option) any later version. See .
*/
+#include "arb_hypgeom.h"
#include "acb_hypgeom.h"
void
@@ -151,9 +152,28 @@ acb_hypgeom_ci_2f3(acb_t res, const acb_t z, slong prec)
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 (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)))
+ {
+ 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
+ {
+ acb_indeterminate(res);
+ }
+ return;
+ }
+
if (acb_hypgeom_u_use_asymp(z, prec))
acb_hypgeom_ci_asymp(res, z, prec);
else
acb_hypgeom_ci_2f3(res, z, prec);
}
-
diff --git a/arb_hypgeom.h b/arb_hypgeom.h
index 8b202035..377626b9 100644
--- a/arb_hypgeom.h
+++ b/arb_hypgeom.h
@@ -100,7 +100,6 @@ void arb_hypgeom_ei(arb_t res, const arb_t z, slong prec);
void _arb_hypgeom_ei_series(arb_ptr g, arb_srcptr h, slong hlen, slong len, slong prec);
void arb_hypgeom_ei_series(arb_poly_t g, const arb_poly_t h, slong len, slong prec);
-
void _arb_hypgeom_si_asymp(arb_t res, const arb_t z, slong N, slong prec);
void _arb_hypgeom_si_1f2(arb_t res, const arb_t z, slong N, slong wp, slong prec);
void arb_hypgeom_si(arb_t res, const arb_t z, slong prec);
@@ -108,7 +107,10 @@ void arb_hypgeom_si(arb_t res, const arb_t z, slong prec);
void _arb_hypgeom_si_series(arb_ptr g, arb_srcptr h, slong hlen, slong len, slong prec);
void arb_hypgeom_si_series(arb_poly_t g, const arb_poly_t h, slong len, slong prec);
+void _arb_hypgeom_ci_asymp(arb_t res, const arb_t z, slong N, slong prec);
+void _arb_hypgeom_ci_2f3(arb_t res, const arb_t z, slong N, slong wp, slong prec);
void arb_hypgeom_ci(arb_t res, const arb_t z, slong prec);
+
void _arb_hypgeom_ci_series(arb_ptr g, arb_srcptr h, slong hlen, slong len, slong prec);
void arb_hypgeom_ci_series(arb_poly_t g, const arb_poly_t h, slong len, slong prec);
diff --git a/arb_hypgeom/ci.c b/arb_hypgeom/ci.c
new file mode 100644
index 00000000..522b6fb2
--- /dev/null
+++ b/arb_hypgeom/ci.c
@@ -0,0 +1,230 @@
+/*
+ Copyright (C) 2021 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 .
+*/
+
+#include "arb_hypgeom.h"
+
+#define LOG2 0.69314718055994531
+#define INV_LOG2 1.4426950408889634074
+#define EXP1 2.7182818284590452
+
+double arf_get_d_log2_abs_approx_clamped(const arf_t x);
+
+void
+_arb_hypgeom_ci_asymp(arb_t res, const arb_t z, slong N, slong prec)
+{
+ arb_t s, c, sz, cz, u;
+ fmpq a[1];
+ slong wp;
+ mag_t err, t;
+
+ N = FLINT_MAX(N, 1);
+
+ arb_init(s);
+ arb_init(c);
+ arb_init(sz);
+ arb_init(cz);
+ arb_init(u);
+ mag_init(err);
+ mag_init(t);
+
+ /* Error is bounded by first omitted term, N! / z^N */
+ arb_get_mag_lower(err, z);
+
+ mag_inv(err, err);
+ mag_pow_ui(err, err, N);
+ mag_fac_ui(t, N);
+ mag_mul(err, err, t);
+
+ wp = prec * 1.001 + 5;
+
+ arb_set(u, z);
+ *fmpq_numref(&a[0]) = 1;
+ *fmpq_denref(&a[0]) = 1;
+ arb_hypgeom_sum_fmpq_imag_arb(c, s, a, 1, NULL, 0, u, 1, N, wp);
+ arb_add_error_mag(c, err);
+ arb_add_error_mag(s, err);
+
+ arb_sin_cos(sz, cz, z, wp);
+ arb_mul(c, c, sz, wp);
+ arb_submul(c, s, cz, wp);
+ arb_div(res, c, z, prec);
+
+ arb_clear(s);
+ arb_clear(c);
+ arb_clear(sz);
+ arb_clear(cz);
+ arb_clear(u);
+ mag_clear(err);
+ mag_clear(t);
+}
+
+void
+_arb_hypgeom_ci_2f3(arb_t res, const arb_t z, slong N, slong wp, slong prec)
+{
+ mag_t err, t;
+ arb_t s, u;
+ fmpq a[1];
+ fmpq b[3];
+
+ N = FLINT_MAX(N, 1);
+
+ mag_init(err);
+ mag_init(t);
+ arb_init(s);
+ arb_init(u);
+
+ arb_sqr(u, z, wp);
+ arb_mul_2exp_si(u, u, -2);
+ arb_neg(u, u);
+
+ *fmpq_numref(&a[0]) = 1;
+ *fmpq_denref(&a[0]) = 1;
+ *fmpq_numref(&b[0]) = 2;
+ *fmpq_denref(&b[0]) = 1;
+ *fmpq_numref(&b[1]) = 2;
+ *fmpq_denref(&b[1]) = 1;
+ *fmpq_numref(&b[2]) = 3;
+ *fmpq_denref(&b[2]) = 2;
+
+ /* Terms are bounded by u^N / (4 (N!)^2) */
+ arb_get_mag(err, u);
+
+ /* u^N */
+ mag_set(t, err);
+ mag_pow_ui(t, t, N);
+
+ /* geometric factor for u/N^2 */
+ mag_div_ui(err, err, N);
+ mag_div_ui(err, err, N);
+ mag_geom_series(err, err, 0);
+ mag_mul(t, t, err);
+
+ /* 1/(N!)^2 */
+ mag_rfac_ui(err, N);
+ mag_mul(err, err, err);
+ mag_mul(err, err, t);
+
+ /* 1/4 */
+ mag_mul_2exp_si(err, err, -2);
+
+ arb_hypgeom_sum_fmpq_arb(s, a, 1, b, 3, u, 0, N, wp);
+
+ arb_add_error_mag(s, err);
+
+ arb_mul(s, s, u, wp);
+
+ arb_log(u, z, wp);
+ arb_add(s, s, u, wp);
+
+ arb_const_euler(u, wp);
+ arb_add(res, s, u, prec);
+
+ mag_clear(err);
+ mag_clear(t);
+ arb_clear(u);
+ arb_clear(s);
+}
+
+void
+arb_hypgeom_ci(arb_t res, const arb_t z, slong prec)
+{
+ slong wp, N, acc;
+ double dz, du;
+
+ if (!arb_is_positive(z) || !arb_is_finite(z))
+ {
+ arb_indeterminate(res);
+ return;
+ }
+
+ if (ARF_IS_LAGOM(arb_midref(z)))
+ {
+ acc = arb_rel_accuracy_bits(z);
+ acc = FLINT_MAX(acc, 0);
+ acc = FLINT_MIN(acc, prec);
+ acc += FLINT_MAX(0, -ARF_EXP(arb_midref(z)));
+ prec = FLINT_MIN(prec, acc + 32);
+ }
+
+ dz = fabs(arf_get_d(arb_midref(z), ARF_RND_DOWN));
+ dz = FLINT_MIN(dz, 1e300);
+
+ if (dz > 2.0)
+ {
+ double log2_err, err_prev, log2dz;
+
+ log2dz = arf_get_d_log2_abs_approx_clamped(arb_midref(z));
+
+ err_prev = 0.0;
+ for (N = 1; N < 2 * prec; N++)
+ {
+ log2_err = ((N + 1.0) * (log(N + 1.0) - 1.0)) * INV_LOG2 - N * log2dz;
+
+ if (log2_err > err_prev)
+ break;
+
+ if (log2_err < -prec - 2)
+ {
+ _arb_hypgeom_ci_asymp(res, z, N, prec);
+ return;
+ }
+
+ err_prev = log2_err;
+ }
+ }
+
+ if (arf_cmpabs_2exp_si(arb_midref(z), -30) < 0)
+ {
+ N = -arf_abs_bound_lt_2exp_si(arb_midref(z));
+ wp = prec * 1.001 + 10;
+ N = (prec + N - 1) / N;
+ }
+ else
+ {
+ du = 0.25 * dz * dz;
+ wp = prec * 1.001 + 10;
+ if (du > 1.0)
+ wp += dz * 1.4426950408889634;
+
+ N = (prec + 5) * LOG2 / (2 * d_lambertw((prec + 5) * LOG2 / (2 * EXP1 * sqrt(du)))) + 1;
+ }
+
+ if (arb_is_exact(z))
+ {
+ _arb_hypgeom_ci_2f3(res, z, N, wp, prec);
+ }
+ else
+ {
+ mag_t err;
+ mag_init(err);
+
+ /* |ci'(z)| = |cos(z)/z| <= 1/z */
+ arb_get_mag_lower(err, z);
+
+ if (mag_cmp_2exp_si(err, 0) >= 0 || 1)
+ {
+ arb_t zmid;
+ arb_init(zmid);
+ arb_get_mid_arb(zmid, z);
+ mag_inv(err, err);
+ mag_mul(err, err, arb_radref(z));
+ _arb_hypgeom_ci_2f3(res, zmid, N, wp, prec);
+ arb_add_error_mag(res, err);
+ arb_clear(zmid);
+ }
+ else
+ {
+ _arb_hypgeom_ci_2f3(res, z, N, wp, prec);
+ }
+
+ mag_clear(err);
+ }
+}
diff --git a/arb_hypgeom/sum_fmpq_arb_rs.c b/arb_hypgeom/sum_fmpq_arb_rs.c
index d14947e0..4f649d10 100644
--- a/arb_hypgeom/sum_fmpq_arb_rs.c
+++ b/arb_hypgeom/sum_fmpq_arb_rs.c
@@ -34,12 +34,12 @@ tail_precision(slong k, double min_k, slong alen, slong blen, double log2z, doub
new_prec = FLINT_MIN(new_prec, prec);
new_prec = FLINT_MAX(new_prec, 32);
-/* printf("term %ld, max %f, log2x %f, magn %f new_prec %ld\n", k, log2z, log2max, term_magnitude, new_prec); */
+ /* printf("term %ld, max %f, log2x %f, magn %f new_prec %ld\n", k, log2z, log2max, term_magnitude, new_prec); */
return new_prec;
}
-/* Return approximation of log2(|x|), clambed between COEFF_MIN and COEFF_MAX. */
+/* Return approximation of log2(|x|), approximately clamped between COEFF_MIN and COEFF_MAX. */
double
arf_get_d_log2_abs_approx_clamped(const arf_t x)
{
@@ -60,12 +60,21 @@ arf_get_d_log2_abs_approx_clamped(const arf_t x)
}
else
{
+ mp_srcptr tp;
+ mp_size_t tn;
+ double v;
slong e = ARF_EXP(x);
- if (e < -50 || e > 50)
- return e;
+ ARF_GET_MPN_READONLY(tp, tn, x);
+
+ if (tn == 1)
+ v = (double)(tp[0]);
else
- return 1.4426950408889634074 * mag_d_log_upper_bound(fabs(arf_get_d(x, ARF_RND_UP)));
+ v = (double)(tp[tn - 1]) + (double)(tp[tn - 2]) * ldexp(1.0, -FLINT_BITS);
+
+ v *= ldexp(1.0, -FLINT_BITS);
+
+ return 1.4426950408889634074 * mag_d_log_upper_bound(v) + (double) e;
}
}
diff --git a/arb_hypgeom/test/t-ci.c b/arb_hypgeom/test/t-ci.c
new file mode 100644
index 00000000..ec9cf767
--- /dev/null
+++ b/arb_hypgeom/test/t-ci.c
@@ -0,0 +1,92 @@
+/*
+ Copyright (C) 2018 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 .
+*/
+
+#include "arb_hypgeom.h"
+
+int main()
+{
+ slong iter;
+ flint_rand_t state;
+
+ flint_printf("ci....");
+ fflush(stdout);
+
+ flint_randinit(state);
+
+ for (iter = 0; iter < 10000 * arb_test_multiplier(); iter++)
+ {
+ arb_t x, s, t;
+ slong prec1, prec2;
+
+ arb_init(x);
+ arb_init(s);
+ arb_init(t);
+
+ if (n_randint(state, 10) == 0)
+ {
+ prec1 = 2 + n_randint(state, 2000);
+ prec2 = 2 + n_randint(state, 2000);
+ }
+ else
+ {
+ prec1 = 2 + n_randint(state, 200);
+ prec2 = 2 + n_randint(state, 200);
+ }
+
+ arb_randtest(x, state, 2 + n_randint(state, prec1), 2 + n_randint(state, 100));
+ arb_randtest(s, state, 2 + n_randint(state, prec1), 2 + n_randint(state, 100));
+ arb_randtest(t, state, 2 + n_randint(state, prec1), 2 + n_randint(state, 100));
+
+ switch (n_randint(state, 3))
+ {
+ case 0:
+ arb_hypgeom_ci(s, x, prec1);
+ break;
+ case 1:
+ _arb_hypgeom_ci_2f3(s, x, n_randint(state, prec1), prec1, prec1);
+ break;
+ default:
+ _arb_hypgeom_c_asymp(s, x, n_randint(state, prec1 / 2), prec1);
+ break;
+ }
+
+ switch (n_randint(state, 3))
+ {
+ case 0:
+ arb_hypgeom_ci(t, x, prec2);
+ break;
+ case 1:
+ _arb_hypgeom_ci_2f3(t, x, n_randint(state, prec2), prec2, prec2);
+ break;
+ default:
+ _arb_hypgeom_ci_asymp(t, x, n_randint(state, prec2 / 2), prec2);
+ break;
+ }
+
+ if (!arb_overlaps(s, t))
+ {
+ flint_printf("FAIL: overlap\n\n");
+ flint_printf("x = "); arb_printn(x, 100, 0); flint_printf("\n\n");
+ flint_printf("s = "); arb_printn(s, 100, 0); flint_printf("\n\n");
+ flint_printf("t = "); arb_printn(t, 100, 0); flint_printf("\n\n");
+ flint_abort();
+ }
+
+ arb_clear(x);
+ arb_clear(s);
+ arb_clear(t);
+ }
+
+ flint_randclear(state);
+ flint_cleanup();
+ flint_printf("PASS\n");
+ return EXIT_SUCCESS;
+}
diff --git a/arb_hypgeom/test/t-si.c b/arb_hypgeom/test/t-si.c
index 517c2b89..1a032c7a 100644
--- a/arb_hypgeom/test/t-si.c
+++ b/arb_hypgeom/test/t-si.c
@@ -71,9 +71,6 @@ int main()
break;
}
- arb_hypgeom_si(s, x, prec1);
- _arb_hypgeom_si_asymp(t, x, n_randint(state, prec2 / 2), prec2);
-
if (!arb_overlaps(s, t))
{
flint_printf("FAIL: overlap\n\n");
diff --git a/arb_hypgeom/wrappers.c b/arb_hypgeom/wrappers.c
index 057e65f2..9940c0a2 100644
--- a/arb_hypgeom/wrappers.c
+++ b/arb_hypgeom/wrappers.c
@@ -70,24 +70,6 @@ arb_hypgeom_ei(arb_t res, const arb_t z, slong prec)
}
}
-void
-arb_hypgeom_ci(arb_t res, const arb_t z, slong prec)
-{
- if (!arb_is_finite(z) || !arb_is_positive(z))
- {
- arb_indeterminate(res);
- }
- else
- {
- acb_t t;
- acb_init(t);
- arb_set(acb_realref(t), z);
- acb_hypgeom_ci(t, t, prec);
- arb_swap(res, acb_realref(t));
- acb_clear(t);
- }
-}
-
void
arb_hypgeom_shi(arb_t res, const arb_t z, slong prec)
{
diff --git a/doc/source/arb_hypgeom.rst b/doc/source/arb_hypgeom.rst
index 99a6047d..b6e0d368 100644
--- a/doc/source/arb_hypgeom.rst
+++ b/doc/source/arb_hypgeom.rst
@@ -422,7 +422,9 @@ Exponential and trigonometric integrals
Computes the sine integral of the power series *z*,
truncated to length *len*.
-.. function:: void arb_hypgeom_ci(arb_t res, const arb_t z, slong prec)
+.. function:: void _arb_hypgeom_ci_asymp(arb_t res, const arb_t z, slong N, slong prec)
+ void _arb_hypgeom_ci_2f3(arb_t res, const arb_t z, slong N, slong wp, slong prec)
+ void arb_hypgeom_ci(arb_t res, const arb_t z, slong prec)
Computes the cosine integral `\operatorname{Ci}(z)`.
The result is indeterminate if `z < 0` since the value of the function would be complex.