hiro/arb
1
0
Fork 0
mirror of https://github.com/vale981/arb synced 2025-03-04 17:01:40 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions

faster dilog implementation

Browse source
This commit is contained in:
Fredrik Johansson 2017-02-24 21:58:23 +01:00
parent 6825a4b912
commit 536ec0faaf
11 changed files with 1265 additions and 22 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

8
acb_hypgeom.h
View file

@ -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

153
acb_hypgeom/dilog.c Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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);
}

92
acb_hypgeom/dilog_bernoulli.c Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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);
}

88
acb_hypgeom/dilog_bitburst.c Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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);
}

262
acb_hypgeom/dilog_continuation.c Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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);
}

192
acb_hypgeom/dilog_transform.c Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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);
}

33
acb_hypgeom/dilog_zero.c Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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);
}
}

236
acb_hypgeom/dilog_zero_taylor.c Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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);
}

128
acb_hypgeom/test/t-dilog.c Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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;
}

42
acb_poly/polylog_series.c
View file

@ -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);

53
doc/source/acb_hypgeom.rst
View file

@ -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.