diff --git a/acb_hypgeom.h b/acb_hypgeom.h
index 4c963653..591ff256 100644
--- a/acb_hypgeom.h
+++ b/acb_hypgeom.h
@@ -244,6 +244,14 @@ void acb_hypgeom_chebyshev_t(acb_t res, const acb_t n, const acb_t z, slong prec
void acb_hypgeom_chebyshev_u(acb_t res, const acb_t n, const acb_t z, slong prec);
void acb_hypgeom_spherical_y(acb_t res, slong n, slong m, const acb_t theta, const acb_t phi, slong prec);
+void acb_hypgeom_dilog_bernoulli(acb_t res, const acb_t z, slong prec);
+void acb_hypgeom_dilog_continuation(acb_t res, const acb_t a, const acb_t z, slong prec);
+void acb_hypgeom_dilog_bitburst(acb_t res, acb_t z0, const acb_t z, slong prec);
+void acb_hypgeom_dilog_transform(acb_t res, const acb_t z, int algorithm, slong prec);
+void acb_hypgeom_dilog_zero_taylor(acb_t res, const acb_t z, slong prec);
+void acb_hypgeom_dilog_zero(acb_t res, const acb_t z, slong prec);
+void acb_hypgeom_dilog(acb_t res, const acb_t z, slong prec);
+
#ifdef __cplusplus
}
#endif
diff --git a/acb_hypgeom/dilog.c b/acb_hypgeom/dilog.c
new file mode 100644
index 00000000..25e9a3e7
--- /dev/null
+++ b/acb_hypgeom/dilog.c
@@ -0,0 +1,153 @@
+/*
+ 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 .
+*/
+
+#include "acb_hypgeom.h"
+
+void
+acb_hypgeom_dilog(acb_t res, const acb_t z, slong prec)
+{
+ double a, b, best, mz, mz1, t, u;
+ int algorithm;
+ slong acc;
+
+ if (!acb_is_finite(z))
+ {
+ acb_indeterminate(res);
+ return;
+ }
+
+ if (acb_is_zero(z))
+ {
+ acb_zero(res);
+ return;
+ }
+
+ acc = acb_rel_accuracy_bits(z);
+ acc = FLINT_MAX(acc, 0);
+ acc = FLINT_MIN(acc, prec);
+ prec = FLINT_MIN(prec, acc + 30);
+
+ /* first take care of exponents that may overflow doubles */
+ if (arf_cmpabs_2exp_si(arb_midref(acb_realref(z)), -20) <= 0 &&
+ arf_cmpabs_2exp_si(arb_midref(acb_imagref(z)), -20) <= 0)
+ {
+ acb_hypgeom_dilog_zero(res, z, prec);
+ return;
+ }
+
+ if (arf_cmpabs_2exp_si(arb_midref(acb_realref(z)), 20) >= 0 ||
+ arf_cmpabs_2exp_si(arb_midref(acb_imagref(z)), 20) >= 0)
+ {
+ acb_hypgeom_dilog_transform(res, z, 1, prec);
+ return;
+ }
+
+ a = arf_get_d(arb_midref(acb_realref(z)), ARF_RND_DOWN);
+ b = arf_get_d(arb_midref(acb_imagref(z)), ARF_RND_DOWN);
+
+ best = mz = a * a + b * b;
+ algorithm = 0;
+
+ /* if |z| > 0.25, consider expanding somewhere other than the origin */
+ if (best > 0.25 * 0.25)
+ {
+ if (1.0 / mz < best) /* use 1/z */
+ {
+ best = 1.0 / mz;
+ algorithm = 1;
+ }
+
+ mz1 = (a - 1.0) * (a - 1.0) + b * b;
+
+ if (mz1 < best) /* use 1-z */
+ {
+ best = mz1;
+ algorithm = 2;
+ }
+
+ if (mz1 > 0.001 && mz / mz1 < best) /* use z/(z-1) */
+ {
+ best = mz / mz1;
+ algorithm = 3;
+ }
+
+ if (mz1 > 0.001 && 1.0 / mz1 < best) /* use 1/(1-z) */
+ {
+ best = 1.0 / mz1;
+ algorithm = 4;
+ }
+ }
+
+ /* do we still have |z| > 0.25 after transforming? */
+ if (best > 0.25 * 0.25)
+ {
+ /* use series with bernoulli numbers (if not too many!) */
+ if (prec < 10000)
+ {
+ /* t = |log(a+bi)|^2 / (2 pi)^2 */
+ t = log(a * a + b * b);
+ u = atan2(b, a);
+ t = (t * t + u * u) * 0.02533029591;
+
+ if (prec > 1000)
+ t *= 4.0; /* penalty at high precision */
+ else
+ t *= 1.1; /* small penalty... also helps avoid this
+ method at negative reals where the log branch
+ cut enters (todo: combine with 1-z formula?) */
+ if (t < best)
+ {
+ algorithm = 8;
+ best = t;
+ }
+ }
+ }
+
+ /* fall back on expanding at another special point
+ (this should only happen at high precision, where we will use
+ the bit-burst algorithm) */
+ if (best > 0.75 * 0.75)
+ {
+ b = fabs(b); /* reduce to upper half plane */
+
+ /* expanding at z0 = i: effective radius |z-i|/sqrt(2) */
+ t = ((b - 1.0) * (b - 1.0) + a * a) * 0.5;
+ if (t < best)
+ {
+ best = t;
+ algorithm = 5;
+ }
+
+ /* expanding at z0 = (1+i)/2: effective radius |z-(1+i)/2|*sqrt(2) */
+ t = 1.0 + 2.0 * (a * (a - 1.0) + b * (b - 1.0));
+ if (t < best)
+ {
+ best = t;
+ algorithm = 6;
+ }
+
+ /* expanding at z0 = 1+i: effective radius |z-(1+i)| */
+ t = 2.0 + (a - 2.0) * a + (b - 2.0) * b;
+ if (t < best)
+ {
+ best = t;
+ algorithm = 7;
+ }
+ }
+
+ if (algorithm == 0)
+ acb_hypgeom_dilog_zero(res, z, prec);
+ else if (algorithm >= 1 && algorithm <= 7)
+ acb_hypgeom_dilog_transform(res, z, algorithm, prec);
+ else /* (algorithm == 8) */
+ acb_hypgeom_dilog_bernoulli(res, z, prec);
+}
+
diff --git a/acb_hypgeom/dilog_bernoulli.c b/acb_hypgeom/dilog_bernoulli.c
new file mode 100644
index 00000000..79cb525f
--- /dev/null
+++ b/acb_hypgeom/dilog_bernoulli.c
@@ -0,0 +1,92 @@
+/*
+ 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 .
+*/
+
+#include "acb_hypgeom.h"
+#include "bernoulli.h"
+
+/* todo: use log(1-z) when this is better? would also need to
+ adjust strategy in the main function */
+void
+acb_hypgeom_dilog_bernoulli(acb_t res, const acb_t z, slong prec)
+{
+ acb_t s, w, w2;
+ slong n, k;
+ fmpz_t c, d;
+ mag_t m, err;
+ double lm;
+
+ acb_init(s);
+ acb_init(w);
+ acb_init(w2);
+ fmpz_init(c);
+ fmpz_init(d);
+ mag_init(m);
+ mag_init(err);
+
+ acb_log(w, z, prec);
+ acb_get_mag(m, w);
+
+ /* for k >= 4, the terms are bounded by (|w| / (2 pi))^k */
+ mag_set_ui_2exp_si(err, 2670177, -24); /* upper bound for 1/(2pi) */
+ mag_mul(err, err, m);
+ lm = mag_get_d_log2_approx(err);
+
+ if (lm < -0.25)
+ {
+ n = prec / (-lm) + 1;
+ n = FLINT_MAX(n, 4);
+ mag_geom_series(err, err, n);
+
+ BERNOULLI_ENSURE_CACHED(n)
+
+ acb_mul(w2, w, w, prec);
+
+ for (k = n - (n % 2 == 0); k >= 3; k -= 2)
+ {
+ fmpz_mul_ui(c, fmpq_denref(bernoulli_cache + k - 1), k - 1);
+ fmpz_mul_ui(d, c, (k + 1) * (k + 2));
+ acb_mul(s, s, w2, prec);
+ acb_mul_fmpz(s, s, c, prec);
+ fmpz_mul_ui(c, fmpq_numref(bernoulli_cache + k - 1), (k + 1) * (k + 2));
+ acb_sub_fmpz(s, s, c, prec);
+ acb_div_fmpz(s, s, d, prec);
+ }
+
+ acb_mul(s, s, w, prec);
+ acb_mul_2exp_si(s, s, 1);
+ acb_sub_ui(s, s, 3, prec);
+ acb_mul(s, s, w2, prec);
+ acb_mul_2exp_si(s, s, -1);
+ acb_const_pi(w2, prec);
+ acb_addmul(s, w2, w2, prec);
+ acb_div_ui(s, s, 6, prec);
+
+ acb_neg(w2, w);
+ acb_log(w2, w2, prec);
+ acb_submul(s, w2, w, prec);
+ acb_add(res, s, w, prec);
+
+ acb_add_error_mag(res, err);
+ }
+ else
+ {
+ acb_indeterminate(res);
+ }
+
+ acb_clear(s);
+ acb_clear(w);
+ acb_clear(w2);
+ fmpz_clear(c);
+ fmpz_clear(d);
+ mag_clear(m);
+ mag_clear(err);
+}
+
diff --git a/acb_hypgeom/dilog_bitburst.c b/acb_hypgeom/dilog_bitburst.c
new file mode 100644
index 00000000..d4dc66d3
--- /dev/null
+++ b/acb_hypgeom/dilog_bitburst.c
@@ -0,0 +1,88 @@
+/*
+ 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 .
+*/
+
+#include "acb_hypgeom.h"
+#include "bernoulli.h"
+
+static void
+_arf_trunc(arf_t x)
+{
+ if (arf_sgn(x) < 0)
+ arf_ceil(x, x);
+ else
+ arf_floor(x, x);
+}
+
+static void
+acb_extract_bits(acb_t t, const acb_t z, slong b)
+{
+ acb_mul_2exp_si(t, z, b);
+ _arf_trunc(arb_midref(acb_realref(t)));
+ _arf_trunc(arb_midref(acb_imagref(t)));
+ mag_zero(arb_radref(acb_realref(t)));
+ mag_zero(arb_radref(acb_imagref(t)));
+ acb_mul_2exp_si(t, t, -b);
+}
+
+void
+acb_hypgeom_dilog_bitburst(acb_t res, acb_t z0, const acb_t z, slong prec)
+{
+ acb_t s, t, tprev, u;
+ slong w;
+ slong start = 16;
+
+ acb_init(s);
+ acb_init(t);
+ acb_init(tprev);
+ acb_init(u);
+
+ acb_sub_ui(t, z, 1, 30);
+ arb_abs(acb_imagref(t), acb_imagref(t));
+
+ /* we don't want to end up on the branch cut */
+ if (arb_contains_nonnegative(acb_realref(t))
+ && !arb_gt(acb_imagref(t), acb_realref(t)))
+ {
+ acb_set(z0, z);
+ acb_zero(res);
+ }
+ else
+ {
+ acb_extract_bits(t, z, start);
+ acb_set(z0, t);
+ acb_set(tprev, t);
+
+ for (w = 2 * start; w < prec; w *= 2)
+ {
+ acb_extract_bits(t, z, w);
+ acb_sub(u, t, z, 30);
+
+ if (arf_cmpabs_2exp_si(arb_midref(acb_realref(u)), -prec / 8) < 0 &&
+ arf_cmpabs_2exp_si(arb_midref(acb_realref(u)), -prec / 8) < 0)
+ break;
+
+ acb_hypgeom_dilog_continuation(u, tprev, t, prec);
+ acb_add(s, s, u, prec);
+ acb_set(tprev, t);
+ }
+
+ acb_hypgeom_dilog_continuation(u, tprev, z, prec);
+ acb_add(s, s, u, prec);
+
+ acb_set(res, s);
+ }
+
+ acb_clear(s);
+ acb_clear(t);
+ acb_clear(tprev);
+ acb_clear(u);
+}
+
diff --git a/acb_hypgeom/dilog_continuation.c b/acb_hypgeom/dilog_continuation.c
new file mode 100644
index 00000000..05d11d30
--- /dev/null
+++ b/acb_hypgeom/dilog_continuation.c
@@ -0,0 +1,262 @@
+/*
+ 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 .
+*/
+
+#include "acb_hypgeom.h"
+
+static void
+bsplit(acb_ptr VA, const acb_t z, const acb_t z2,
+ const acb_t a, const acb_t a1a, slong k, slong h, slong prec)
+{
+ if (h - k == 1)
+ {
+ acb_zero(VA + 0);
+ acb_mul_ui(VA + 1, a1a, (k+1)*(k+2), prec);
+ acb_mul_si(VA + 2, z2, -k*k, prec);
+ acb_mul_ui(VA + 3, a, (k+1)*(2*k+1), prec);
+ acb_sub_ui(VA + 3, VA + 3, (k+1)*(k+1), prec);
+ acb_mul(VA + 3, VA + 3, z, prec);
+ acb_neg(VA + 3, VA + 3);
+ acb_set(VA + 4, VA + 1);
+ acb_zero(VA + 5);
+ acb_set(VA + 6, VA + 1);
+ }
+ else
+ {
+ slong m;
+ acb_ptr VB;
+
+ if (h <= k) abort();
+
+ m = k + (h - k) / 2;
+ VB = _acb_vec_init(7);
+
+ bsplit(VA, z, z2, a, a1a, k, m, prec);
+ bsplit(VB, z, z2, a, a1a, m, h, prec);
+
+ acb_mul(VA + 6, VA + 6, VB + 6, prec);
+ acb_mul(VA + 4, VA + 4, VB + 6, prec);
+ acb_addmul(VA + 4, VA + 0, VB + 4, prec);
+ acb_addmul(VA + 4, VA + 2, VB + 5, prec);
+ acb_mul(VA + 5, VA + 5, VB + 6, prec);
+ acb_addmul(VA + 5, VA + 1, VB + 4, prec);
+ acb_addmul(VA + 5, VA + 3, VB + 5, prec);
+ acb_set(VB + 6, VA + 3);
+ acb_mul(VA + 3, VA + 3, VB + 3, prec);
+ acb_addmul(VA + 3, VA + 1, VB + 2, prec);
+ acb_set(VB + 5, VA + 2);
+ acb_mul(VA + 2, VA + 2, VB + 3, prec);
+ acb_addmul(VA + 2, VA + 0, VB + 2, prec);
+ acb_mul(VA + 1, VA + 1, VB + 0, prec);
+ acb_addmul(VA + 1, VB + 6, VB + 1, prec);
+ acb_mul(VA + 0, VA + 0, VB + 0, prec);
+ acb_addmul(VA + 0, VB + 5, VB + 1, prec);
+
+ _acb_vec_clear(VB, 7);
+ }
+}
+
+/*
+Some possible approaches to bounding the Taylor coefficients c_k at
+the expansion point a:
+
+1. Inspection of the symbolic derivatives gives the trivial bound
+
+ |c_k| <= (1+|log(1-a)|) / min(|a|,|a-1|)^k
+
+ which is good enough when not close to 0.
+
+2. Using Cauchy's integral formula, some explicit computation gives
+ |c_k| <= 4/|1-a|^k when |a| <= 1/2 or a = +/- i. The constant
+ could certainly be improved.
+
+3. For k >= 1, c_k = 2F1(k,k,k+1,a) / k^2. Can we use monotonicity to get
+ good estimates here when a is complex? Note that 2F1(k,k,k,a) = (1-a)^-k.
+
+*/
+
+void
+acb_hypgeom_dilog_continuation(acb_t res, const acb_t a, const acb_t z, slong prec)
+{
+ acb_t za, a1, a1a, za2, log1a;
+ acb_ptr V;
+ slong n;
+ double tr;
+ mag_t tm, err, am;
+ int real;
+
+ if (acb_is_zero(a))
+ {
+ acb_hypgeom_dilog_zero_taylor(res, z, prec);
+ return;
+ }
+
+ if (acb_eq(a, z))
+ {
+ acb_zero(res);
+ return;
+ }
+
+ acb_init(za);
+ acb_init(a1);
+ acb_init(a1a);
+ acb_init(za2);
+ acb_init(log1a);
+ mag_init(tm);
+ mag_init(err);
+ mag_init(am);
+
+ acb_sub(za, z, a, prec); /* z-a */
+ acb_sub_ui(a1, a, 1, prec); /* a-1 */
+ acb_mul(a1a, a1, a, prec); /* (a-1)a */
+ acb_mul(za2, za, za, prec); /* (z-a)^2 */
+ acb_neg(log1a, a1);
+ acb_log(log1a, log1a, prec); /* log(1-a) */
+
+ acb_get_mag(am, a);
+ if (mag_cmp_2exp_si(am, -1) <= 0 ||
+ (acb_is_exact(a) && arb_is_zero(acb_realref(a)) &&
+ arf_cmpabs_ui(arb_midref(acb_imagref(a)), 1) == 0))
+ {
+ acb_get_mag_lower(am, a1);
+ acb_get_mag(tm, za);
+ mag_div(tm, tm, am); /* tm = ratio */
+ mag_set_ui(am, 4); /* am = prefactor */
+ }
+ else
+ {
+ acb_get_mag_lower(am, a);
+ acb_get_mag_lower(tm, a1);
+ mag_min(am, am, tm);
+ acb_get_mag(tm, za);
+ mag_div(tm, tm, am); /* tm = ratio */
+ acb_get_mag(am, log1a);
+ mag_add_ui(am, am, 1); /* am = prefactor */
+ }
+
+ tr = mag_get_d_log2_approx(tm);
+ if (tr < -0.1)
+ {
+ arf_srcptr rr, ii;
+
+ rr = arb_midref(acb_realref(z));
+ ii = arb_midref(acb_imagref(z));
+ if (arf_cmpabs(ii, rr) > 0)
+ rr = ii;
+
+ /* target relative accuracy near 0 */
+ if (arf_cmpabs_2exp_si(rr, -2) < 0 && arf_cmpabs_2exp_si(rr, -prec) > 0)
+ n = (prec - arf_abs_bound_lt_2exp_si(rr)) / (-tr) + 1;
+ else
+ n = prec / (-tr) + 1;
+
+ n = FLINT_MAX(n, 2);
+ }
+ else
+ {
+ n = 2;
+ }
+
+ mag_geom_series(err, tm, n);
+ mag_mul(err, err, am);
+
+ real = acb_is_real(a) && acb_is_real(z) &&
+ arb_is_negative(acb_realref(a1)) &&
+ mag_is_finite(err);
+
+ if (n < 10)
+ {
+ /* forward recurrence - faster for small n and/or low precision, but
+ must be avoided for large n since complex intervals blow up */
+ acb_t s, t, u, v;
+ slong k;
+
+ acb_init(s);
+ acb_init(t);
+ acb_init(u);
+ acb_init(v);
+
+ acb_div(u, log1a, a, prec);
+ acb_neg(u, u);
+
+ if (n >= 3)
+ {
+ acb_inv(v, a1, prec);
+ acb_add(v, v, u, prec);
+ acb_div(v, v, a, prec);
+ acb_mul_2exp_si(v, v, -1);
+ acb_neg(v, v);
+ acb_mul(v, v, za2, prec);
+ }
+
+ acb_mul(u, u, za, prec);
+ acb_add(s, u, v, prec);
+
+ for (k = 3; k < n; k++)
+ {
+ acb_mul_ui(u, u, (k - 2) * (k - 2), prec);
+ acb_mul(u, u, za2, prec);
+ acb_mul_ui(t, a, (k - 1) * (2 * k - 3), prec);
+ acb_sub_ui(t, t, (k - 1) * (k - 1), prec);
+ acb_mul(t, t, v, prec);
+ acb_addmul(u, t, za, prec);
+ acb_mul_ui(t, a1a, (k - 1) * k, prec);
+ acb_neg(t, t);
+ acb_div(u, u, t, prec);
+ acb_add(s, s, u, prec);
+ acb_swap(v, u);
+ }
+
+ acb_set(res, s);
+
+ acb_clear(s);
+ acb_clear(t);
+ acb_clear(u);
+ acb_clear(v);
+ }
+ else
+ {
+ /* binary splitting */
+ V = _acb_vec_init(7);
+
+ bsplit(V, za, za2, a, a1a, 1, 1 + n, prec);
+
+ acb_mul(V + 1, V + 4, log1a, prec);
+ acb_neg(V + 1, V + 1);
+ acb_mul(V + 2, V + 5, za2, prec);
+ acb_mul_2exp_si(V + 2, V + 2, -1);
+ acb_mul(V + 1, V + 1, za, prec);
+ acb_div(V + 3, V + 2, a1, prec);
+ acb_sub(V + 1, V + 1, V + 3, prec);
+ acb_div(V + 0, log1a, a, prec);
+ acb_addmul(V + 1, V + 2, V + 0, prec);
+ acb_mul(V + 6, V + 6, a, prec);
+
+ acb_div(V + 0, V + 1, V + 6, prec);
+ acb_set(res, V + 0);
+
+ _acb_vec_clear(V, 7);
+ }
+
+ if (real)
+ arb_add_error_mag(acb_realref(res), err);
+ else
+ acb_add_error_mag(res, err);
+
+ acb_clear(za);
+ acb_clear(a1);
+ acb_clear(a1a);
+ acb_clear(za2);
+ acb_clear(log1a);
+ mag_clear(tm);
+ mag_clear(err);
+ mag_clear(am);
+}
+
diff --git a/acb_hypgeom/dilog_transform.c b/acb_hypgeom/dilog_transform.c
new file mode 100644
index 00000000..8ebb25d3
--- /dev/null
+++ b/acb_hypgeom/dilog_transform.c
@@ -0,0 +1,192 @@
+/*
+ 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 .
+*/
+
+#include "acb_hypgeom.h"
+
+void
+acb_hypgeom_dilog_transform(acb_t res, const acb_t z, int algorithm, slong prec)
+{
+ acb_t t, u;
+
+ acb_init(t);
+ acb_init(u);
+
+ if (algorithm == 1)
+ {
+ /* Li_2(z) = -Li_2(1/z) - log(-z)^2/2 - pi^2/6, z not in (0,1) */
+ arf_set_ui_2exp_si(arb_midref(acb_realref(t)), 1, -1);
+ mag_set_ui_2exp_si(arb_radref(acb_realref(t)), 1, -1);
+
+ if (acb_overlaps(z, t))
+ {
+ acb_indeterminate(res);
+ }
+ else
+ {
+ acb_inv(t, z, prec);
+ acb_hypgeom_dilog_zero(t, t, prec);
+ acb_neg(u, z);
+ acb_log(u, u, prec);
+ acb_mul(u, u, u, prec);
+ acb_mul_2exp_si(u, u, -1);
+ acb_add(t, t, u, prec);
+ acb_const_pi(u, prec);
+ acb_mul(u, u, u, prec);
+ acb_div_ui(u, u, 6, prec);
+ acb_add(t, t, u, prec);
+ acb_neg(res, t);
+ }
+ }
+ else if (algorithm == 2)
+ {
+ /* Li_2(z) = -Li_2(1-z) - log(1-z) log(z) + pi^2/6 */
+ if (acb_is_one(z))
+ {
+ acb_zero(res);
+ }
+ else
+ {
+ acb_sub_ui(t, z, 1, prec);
+ acb_neg(t, t);
+ acb_hypgeom_dilog_zero(u, t, prec);
+ acb_log(t, t, prec);
+ acb_log(res, z, prec);
+ acb_mul(res, res, t, prec);
+ acb_add(res, res, u, prec);
+ }
+
+ acb_const_pi(t, prec);
+ acb_mul(t, t, t, prec);
+ acb_div_ui(t, t, 6, prec);
+ acb_sub(res, t, res, prec);
+ }
+ else if (algorithm == 3)
+ {
+ /* Li_2(z) = -Li_2(z/(z-1)) - log(1-z)^2/2, z not in (1,inf) */
+ acb_sub_ui(t, z, 1, prec);
+
+ if (!arb_is_negative(acb_realref(t)))
+ {
+ acb_indeterminate(res);
+ }
+ else
+ {
+ acb_div(u, z, t, prec);
+ acb_hypgeom_dilog_zero(u, u, prec);
+ acb_neg(t, t);
+ acb_log(t, t, prec);
+ acb_mul(t, t, t, prec);
+ acb_mul_2exp_si(t, t, -1);
+ acb_add(t, t, u, prec);
+ acb_neg(res, t);
+ }
+ }
+ else if (algorithm == 4)
+ {
+ /* Li_2(z) = Li_2(1/(1-z)) + log(1-z) [log(1-z)/2 - log(-z)] - pi^2/6 */
+ acb_sub_ui(t, z, 1, prec);
+ acb_neg(t, t);
+
+ acb_inv(u, t, prec);
+ acb_hypgeom_dilog_zero(u, u, prec);
+
+ acb_log(t, t, prec);
+ acb_neg(res, z);
+ acb_log(res, res, prec);
+ acb_mul_2exp_si(res, res, 1);
+ acb_sub(res, t, res, prec);
+ acb_mul_2exp_si(res, res, -1);
+ acb_addmul(u, res, t, prec);
+
+ acb_const_pi(t, prec);
+ acb_mul(t, t, t, prec);
+ acb_div_ui(t, t, 6, prec);
+ acb_sub(res, u, t, prec);
+ }
+ else if (algorithm >= 5 && algorithm <= 7)
+ {
+ if (arb_contains_zero(acb_imagref(z)))
+ {
+ acb_indeterminate(res);
+ }
+ else
+ {
+ acb_t a;
+ acb_init(a);
+
+ if (algorithm == 5)
+ {
+ acb_onei(a);
+ /* Li_2(i) = -pi^2/48 + C i */
+ arb_const_catalan(acb_imagref(u), prec);
+ arb_const_pi(acb_realref(u), prec);
+ arb_mul(acb_realref(u), acb_realref(u), acb_realref(u), prec);
+ arb_div_si(acb_realref(u), acb_realref(u), -48, prec);
+ }
+ else if (algorithm == 6)
+ {
+ /* Li_2((1+i)/2) = (5 pi^2 / 96 - log(2)^2/8) + (C - pi log(2) / 8) i */
+ arb_t t;
+ arb_init(t);
+ acb_set_d_d(a, 0.5, 0.5);
+ arb_const_pi(t, prec);
+ arb_const_log2(acb_imagref(u), prec);
+ arb_mul(acb_realref(u), acb_imagref(u), acb_imagref(u), prec);
+ arb_mul(acb_imagref(u), acb_imagref(u), t, prec);
+ acb_mul_2exp_si(u, u, -3);
+ arb_mul(t, t, t, prec);
+ arb_mul_ui(t, t, 5, prec);
+ arb_div_ui(t, t, 96, prec);
+ arb_sub(acb_realref(u), t, acb_realref(u), prec);
+ arb_const_catalan(t, prec);
+ arb_sub(acb_imagref(u), t, acb_imagref(u), prec);
+ arb_clear(t);
+ }
+ else
+ {
+ /* Li_2(1+i) = pi^2/16 + (C + pi log(2)/4) i */
+ arb_t t;
+ arb_init(t);
+ acb_set_d_d(a, 1.0, 1.0);
+ arb_const_pi(acb_realref(u), prec);
+ arb_mul_2exp_si(acb_realref(u), acb_realref(u), -2);
+ arb_const_log2(t, prec);
+ arb_mul(acb_imagref(u), acb_realref(u), t, prec);
+ arb_const_catalan(t, prec);
+ arb_add(acb_imagref(u), acb_imagref(u), t, prec);
+ arb_mul(acb_realref(u), acb_realref(u), acb_realref(u), prec);
+ arb_clear(t);
+ }
+
+ if (arf_sgn(arb_midref(acb_imagref(z))) < 0)
+ {
+ acb_conj(a, a);
+ acb_conj(u, u);
+ }
+
+ acb_hypgeom_dilog_bitburst(res, t, z, prec);
+ acb_add(res, res, u, prec);
+ acb_hypgeom_dilog_continuation(t, a, t, prec);
+ acb_add(res, res, t, prec);
+
+ acb_clear(a);
+ }
+ }
+ else
+ {
+ flint_printf("unknown algorithm\n");
+ abort();
+ }
+
+ acb_clear(t);
+ acb_clear(u);
+}
+
diff --git a/acb_hypgeom/dilog_zero.c b/acb_hypgeom/dilog_zero.c
new file mode 100644
index 00000000..eb4315c2
--- /dev/null
+++ b/acb_hypgeom/dilog_zero.c
@@ -0,0 +1,33 @@
+/*
+ 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 .
+*/
+
+#include "acb_hypgeom.h"
+
+void
+acb_hypgeom_dilog_zero(acb_t res, const acb_t z, slong prec)
+{
+ if (prec < 40000 ||
+ (arf_cmpabs_2exp_si(arb_midref(acb_realref(z)), -prec / 1000) < 0
+ && arf_cmpabs_2exp_si(arb_midref(acb_imagref(z)), -prec / 1000) < 0))
+ {
+ acb_hypgeom_dilog_zero_taylor(res, z, prec);
+ }
+ else
+ {
+ acb_t z0;
+ acb_init(z0);
+ acb_hypgeom_dilog_bitburst(res, z0, z, prec);
+ acb_hypgeom_dilog_zero_taylor(z0, z0, prec);
+ acb_add(res, res, z0, prec);
+ acb_clear(z0);
+ }
+}
+
diff --git a/acb_hypgeom/dilog_zero_taylor.c b/acb_hypgeom/dilog_zero_taylor.c
new file mode 100644
index 00000000..cb61af05
--- /dev/null
+++ b/acb_hypgeom/dilog_zero_taylor.c
@@ -0,0 +1,236 @@
+/*
+ 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 .
+*/
+
+#include "acb_hypgeom.h"
+
+static void
+bsplit_zero(acb_t P, acb_t R, acb_t Q, const acb_t z,
+ slong a, slong b, slong prec)
+{
+ if (b - a == 1)
+ {
+ acb_mul_ui(P, z, a * a, prec);
+ acb_set_ui(R, (a + 1) * (a + 1));
+ acb_set(Q, R);
+ }
+ else
+ {
+ acb_t P2, R2, Q2;
+ slong m;
+
+ acb_init(P2);
+ acb_init(R2);
+ acb_init(Q2);
+
+ m = a + (b - a) / 2;
+
+ bsplit_zero(P, R, Q, z, a, m, prec);
+ bsplit_zero(P2, R2, Q2, z, m, b, prec);
+
+ acb_mul(R, R, Q2, prec);
+ acb_addmul(R, P, R2, prec);
+ acb_mul(P, P, P2, prec);
+ acb_mul(Q, Q, Q2, prec);
+
+ acb_clear(P2);
+ acb_clear(R2);
+ acb_clear(Q2);
+ }
+}
+
+#if FLINT_BITS == 64
+#define DIVLIM 1625 /* just a tuning value */
+#define DIVLIM2 1000000000 /* avoid overflow */
+#else
+#define DIVLIM 40
+#define DIVLIM2 30000
+#endif
+
+void
+acb_hypgeom_dilog_taylor_sum(acb_t res, const acb_t z, slong n, slong prec)
+{
+ slong k, qk, m, power;
+ ulong q;
+ acb_t s, t, u;
+ acb_ptr zpow;
+ int real;
+
+ if (n <= 3)
+ {
+ if (n <= 1)
+ acb_zero(res);
+ else if (n == 2)
+ acb_set_round(res, z, prec);
+ else
+ {
+ acb_init(t);
+ acb_mul(t, z, z, prec);
+ acb_mul_2exp_si(t, t, -2);
+ acb_add(res, z, t, prec);
+ acb_clear(t);
+ }
+ return;
+ }
+
+ /* use binary splitting */
+ if (prec > 4000 && acb_bits(z) < prec * 0.02)
+ {
+ acb_init(s);
+ acb_init(t);
+ acb_init(u);
+ bsplit_zero(s, t, u, z, 1, n, prec);
+ acb_add(s, s, t, prec);
+ acb_mul(s, s, z, prec);
+ acb_div(res, s, u, prec);
+ acb_clear(s);
+ acb_clear(t);
+ acb_clear(u);
+ return;
+ }
+
+ real = acb_is_real(z);
+ k = n - 1;
+ m = n_sqrt(n);
+
+ acb_init(s);
+ acb_init(t);
+ zpow = _acb_vec_init(m + 1);
+
+ _acb_vec_set_powers(zpow, z, m + 1, prec);
+
+ power = (n - 1) % m;
+
+ while (k >= DIVLIM)
+ {
+ if (k < DIVLIM2) /* todo: write a div_uiui function? */
+ {
+ acb_div_ui(t, zpow + power, k * k, prec);
+ }
+ else
+ {
+ acb_div_ui(t, zpow + power, k, prec);
+ acb_div_ui(t, t, k, prec);
+ }
+
+ acb_add(s, s, t, prec);
+
+ if (power == 0)
+ {
+ acb_mul(s, s, zpow + m, prec);
+ power = m - 1;
+ }
+ else
+ {
+ power--;
+ }
+
+ k--;
+ }
+
+ qk = k; /* k at which to change denominator */
+ q = 1;
+
+ while (k >= 2)
+ {
+ /* find next qk such that the consecutive denominators can be
+ collected in a single word */
+ if (k == qk)
+ {
+ if (q != 1)
+ acb_div_ui(s, s, q, prec);
+
+ q = qk * qk;
+ qk--;
+
+
+ while (qk > 1)
+ {
+ ulong hi, lo;
+ umul_ppmm(hi, lo, q, qk * qk);
+ if (hi != 0)
+ break;
+ q *= qk * qk;
+ qk--;
+ }
+
+ acb_mul_ui(s, s, q, prec);
+ }
+
+ if (power == 0)
+ {
+ acb_add_ui(s, s, q / (k * k), prec);
+ acb_mul(s, s, zpow + m, prec);
+ power = m - 1;
+ }
+ else
+ {
+ if (real) /* minor optimization */
+ arb_addmul_ui(acb_realref(s), acb_realref(zpow + power), q / (k * k), prec);
+ else
+ acb_addmul_ui(s, zpow + power, q / (k * k), prec);
+ power--;
+ }
+
+ k--;
+ }
+
+ if (q != 1)
+ acb_div_ui(s, s, q, prec);
+ acb_add(s, s, z, prec);
+ acb_swap(res, s);
+
+ _acb_vec_clear(zpow, m + 1);
+ acb_clear(s);
+ acb_clear(t);
+}
+
+void
+acb_hypgeom_dilog_zero_taylor(acb_t res, const acb_t z, slong prec)
+{
+ mag_t m;
+ slong n;
+ double x;
+ int real;
+
+ mag_init(m);
+ acb_get_mag(m, z);
+ real = acb_is_real(z);
+ x = -mag_get_d_log2_approx(m);
+
+ n = 2;
+ if (x > 0.01)
+ {
+ n = prec / x + 1;
+ n += (x > 2.0); /* relative error for very small |z| */
+ }
+
+ n = FLINT_MAX(n, 2);
+
+ mag_geom_series(m, m, n);
+ mag_div_ui(m, m, n);
+ mag_div_ui(m, m, n);
+
+ if (mag_is_finite(m))
+ {
+ acb_hypgeom_dilog_taylor_sum(res, z, n, prec);
+ if (real)
+ arb_add_error_mag(acb_realref(res), m);
+ else
+ acb_add_error_mag(res, m);
+ }
+ else
+ {
+ acb_indeterminate(res);
+ }
+
+ mag_clear(m);
+}
+
diff --git a/acb_hypgeom/test/t-dilog.c b/acb_hypgeom/test/t-dilog.c
new file mode 100644
index 00000000..7443a721
--- /dev/null
+++ b/acb_hypgeom/test/t-dilog.c
@@ -0,0 +1,128 @@
+/*
+ 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 .
+*/
+
+#include "acb_hypgeom.h"
+
+int main()
+{
+ slong iter;
+ flint_rand_t state;
+
+ flint_printf("dilog....");
+ fflush(stdout);
+
+ flint_randinit(state);
+
+ for (iter = 0; iter < 10000 * arb_test_multiplier(); iter++)
+ {
+ acb_t z0, z1, w0, w1, a;
+ slong prec0, prec1;
+ int alg0, alg1;
+
+ acb_init(z0);
+ acb_init(z1);
+ acb_init(w0);
+ acb_init(w1);
+ acb_init(a);
+
+ prec0 = 2 + n_randint(state, 500);
+ prec1 = 2 + n_randint(state, 500);
+
+ acb_randtest(z0, state, 1 + n_randint(state, 500), 1 + n_randint(state, 100));
+ acb_randtest(w0, state, 1 + n_randint(state, 500), 1 + n_randint(state, 100));
+ acb_randtest(w1, state, 1 + n_randint(state, 500), 1 + n_randint(state, 100));
+
+ acb_set(z1, z0);
+ if (n_randint(state, 2))
+ {
+ acb_add(z1, z1, w0, prec0);
+ acb_sub(z1, z1, w0, prec0);
+ }
+
+ switch ((alg0 = n_randint(state, 15)))
+ {
+ case 0:
+ acb_hypgeom_dilog_zero(w0, z0, prec0);
+ break;
+ case 1:
+ case 2:
+ case 3:
+ case 4:
+ case 5:
+ case 6:
+ case 7:
+ acb_hypgeom_dilog_transform(w0, z0, alg0, prec0);
+ break;
+ case 8:
+ acb_hypgeom_dilog_bernoulli(w0, z0, prec0);
+ break;
+ case 9:
+ acb_hypgeom_dilog_bitburst(w0, a, z0, prec0);
+ acb_hypgeom_dilog(a, a, prec0);
+ acb_add(w0, w0, a, prec0);
+ break;
+ default:
+ acb_hypgeom_dilog(w0, z0, prec0);
+ }
+
+ acb_set(w1, z1); /* also test aliasing */
+
+ switch ((alg1 = n_randint(state, 15)))
+ {
+ case 0:
+ acb_hypgeom_dilog_zero(w1, w1, prec1);
+ break;
+ case 1:
+ case 2:
+ case 3:
+ case 4:
+ case 5:
+ case 6:
+ case 7:
+ acb_hypgeom_dilog_transform(w1, w1, alg1, prec1);
+ break;
+ case 8:
+ acb_hypgeom_dilog_bernoulli(w1, z1, prec1);
+ break;
+ case 9:
+ acb_hypgeom_dilog_bitburst(w1, a, z1, prec1);
+ acb_hypgeom_dilog(a, a, prec1);
+ acb_add(w1, w1, a, prec1);
+ break;
+ default:
+ acb_hypgeom_dilog(w1, z1, prec1);
+ }
+
+ if (!acb_overlaps(w0, w1))
+ {
+ flint_printf("FAIL: consistency\n\n");
+ flint_printf("alg0 = %d, alg1 = %d\n\n", alg0, alg1);
+ flint_printf("prec0 = %wd, prec1 = %wd\n\n", prec0, prec1);
+ flint_printf("z0 = "); acb_printd(z0, 30); flint_printf("\n\n");
+ flint_printf("z1 = "); acb_printd(z1, 30); flint_printf("\n\n");
+ flint_printf("w0 = "); acb_printd(w0, 30); flint_printf("\n\n");
+ flint_printf("w1 = "); acb_printd(w1, 30); flint_printf("\n\n");
+ abort();
+ }
+
+ acb_clear(z0);
+ acb_clear(z1);
+ acb_clear(w0);
+ acb_clear(w1);
+ acb_clear(a);
+ }
+
+ flint_randclear(state);
+ flint_cleanup();
+ flint_printf("PASS\n");
+ return EXIT_SUCCESS;
+}
+
diff --git a/acb_poly/polylog_series.c b/acb_poly/polylog_series.c
index 9a4bb086..7cfc79ab 100644
--- a/acb_poly/polylog_series.c
+++ b/acb_poly/polylog_series.c
@@ -10,6 +10,7 @@
*/
#include "acb_poly.h"
+#include "acb_hypgeom.h"
/* note: will not return a wrong value, as arf_get_si aborts on overflow */
slong
@@ -57,30 +58,15 @@ int
polylog_is_real(const acb_t s, const acb_t z)
{
if (!arb_is_zero(acb_imagref(s)))
- {
return 0;
- }
else if (!arb_is_zero(acb_imagref(z)))
- {
return 0;
- }
+ else if (arb_contains_si(acb_realref(z), 1))
+ return 0;
+ else if (acb_is_int(s) && arb_is_nonpositive(acb_realref(s)))
+ return 1;
else
- {
- fmpz_t one;
- int res;
-
- fmpz_init(one);
- fmpz_one(one);
-
- if (arb_contains_fmpz(acb_realref(z), one))
- res = 0;
- else
- res = (arf_cmp_2exp_si(arb_midref(acb_realref(z)), 0) < 0);
-
- fmpz_clear(one);
-
- return res;
- }
+ return (arf_cmp_2exp_si(arb_midref(acb_realref(z)), 0) < 0);
}
void
@@ -89,7 +75,7 @@ _acb_poly_polylog_cpx_zeta(acb_ptr w, const acb_t s, const acb_t z, slong len, s
acb_ptr e1, e2, z1, z2, e1z1, e2z2;
acb_t t, u, v;
slong k, len2;
- int deflate_zeta, deflate_gamma;
+ int deflate_zeta, deflate_gamma, is_real;
if (!acb_is_finite(s) || !acb_is_finite(z))
{
@@ -112,6 +98,8 @@ _acb_poly_polylog_cpx_zeta(acb_ptr w, const acb_t s, const acb_t z, slong len, s
return;
}
+ is_real = polylog_is_real(s, z);
+
acb_init(t);
acb_init(u);
acb_init(v);
@@ -218,6 +206,10 @@ _acb_poly_polylog_cpx_zeta(acb_ptr w, const acb_t s, const acb_t z, slong len, s
for (k = 1; k < len; k += 2)
acb_neg(w + k, w + k);
+ if (is_real)
+ if (acb_is_finite(w))
+ arb_zero(acb_imagref(w));
+
_acb_vec_clear(e1, len + 1);
_acb_vec_clear(e2, len + 1);
_acb_vec_clear(z1, len + 1);
@@ -260,7 +252,7 @@ _acb_poly_polylog_cpx_small(acb_ptr w, const acb_t s, const acb_t z, slong len,
mag_rfac_ui(errf, k);
mag_mul(err, err, errf);
- if (is_real)
+ if (is_real && mag_is_finite(err))
arb_add_error_mag(acb_realref(w + k), err);
else
acb_add_error_mag(w + k, err);
@@ -277,6 +269,12 @@ _acb_poly_polylog_cpx(acb_ptr w, const acb_t s, const acb_t z, slong len, slong
{
mag_t zmag;
+ if (len == 1 && acb_equal_si(s, 2))
+ {
+ acb_hypgeom_dilog(w, z, prec);
+ return;
+ }
+
mag_init(zmag);
acb_get_mag(zmag, z);
diff --git a/doc/source/acb_hypgeom.rst b/doc/source/acb_hypgeom.rst
index 535e5c68..d6aceac3 100644
--- a/doc/source/acb_hypgeom.rst
+++ b/doc/source/acb_hypgeom.rst
@@ -1076,3 +1076,56 @@ Orthogonal polynomials and functions
This function is a polynomial in `\cos(\theta)` and `\sin(\theta)`.
We evaluate it using :func:`acb_hypgeom_legendre_p_uiui_rec`.
+Dilogarithm
+-------------------------------------------------------------------------------
+
+The dilogarithm function
+is given by `\operatorname{Li}_2(z) = -\int_0^z \frac{\log(1-t)}{t} dt = z {}_3F_2(1,1,1,2,2,z)`.
+
+.. function :: void acb_hypgeom_dilog_bernoulli(acb_t res, const acb_t z, slong prec)
+
+ Computes the dilogarithm using a series expansion in `w = \log(z)`,
+ with rate of convergence `|w/(2\pi)|^n`. This provides good convergence
+ near `z = e^{\pm i \pi / 3}`, where hypergeometric series expansions fail.
+ Since the coefficients involve Bernoulli numbers, this method should
+ only be used at moderate precision.
+
+.. function:: void acb_hypgeom_dilog_zero_taylor(acb_t res, const acb_t z, slong prec)
+
+ Computes the dilogarithm for *z* close to 0 using the hypergeometric series
+ (effective only when `|z| \ll 1`).
+
+.. function:: void acb_hypgeom_dilog_zero(acb_t res, const acb_t z, slong prec)
+
+ Computes the dilogarithm for *z* close to 0, using the bit-burst algorithm
+ instead of the hypergeometric series directly at very high precision.
+
+.. function:: void acb_hypgeom_dilog_transform(acb_t res, const acb_t z, int algorithm, slong prec)
+
+ Computes the dilogarithm by applying one of the transformations
+ `1/z`, `1-z`, `z/(z-1)`, `1/(1-z)`, indexed by *algorithm* from 1 to 4,
+ and calling :func:`acb_hypgeom_dilog_zero` with the reduced variable.
+ Alternatively, for *algorithm* between 5 and 7, starts from the
+ respective point `\pm i`, `(1\pm i)/2`, `(1\pm i)/2` (with the sign
+ chosen according to the midpoint of *z*)
+ and computes the dilogarithm by the bit-burst method.
+
+.. function:: void acb_hypgeom_dilog_continuation(acb_t res, const acb_t a, const acb_t z, slong prec)
+
+ Computes `\operatorname{Li}_2(z) - \operatorname{Li}_2(a)` using
+ Taylor expansion at *a*. Binary splitting is used. Both *a* and *z*
+ should be well isolated from the points 0 and 1, except that *a* may
+ be exactly 0. If the straight line path from *a* to *b* crosses the branch
+ cut, this method provides continuous analytic continuation instead of
+ computing the principal branch.
+
+.. function:: void acb_hypgeom_dilog_bitburst(acb_t res, acb_t z0, const acb_t z, slong prec)
+
+ Sets *z0* to a point with short bit expansion close to *z* and sets
+ *res* to `\operatorname{Li}_2(z) - \operatorname{Li}_2(z_0)`, computed
+ using the bit-burst algorithm.
+
+.. function:: void acb_hypgeom_dilog(acb_t res, const acb_t z, slong prec)
+
+ Computes the dilogarithm using a default algorithm choice.
+