mirror of
https://github.com/vale981/arb
synced 2025-03-05 09:21:38 -05:00
initial code for polylogarithms
This commit is contained in:
Fredrik Johansson
parent
4cc48bd0a1
commit
8ad408f251
7 changed files with 947 additions and 24 deletions
23
acb.h
23
acb.h
|
|
@ -256,6 +256,14 @@ acb_add_error_arf(acb_t x, const arf_t err)
|
|||
arb_add_error_arf(acb_imagref(x), err);
|
||||
}
|
||||
|
||||
/* TODO: document */
|
||||
static __inline__ void
|
||||
acb_add_error_mag(acb_t x, const mag_t err)
|
||||
{
|
||||
arb_add_error_mag(acb_realref(x), err);
|
||||
arb_add_error_mag(acb_imagref(x), err);
|
||||
}
|
||||
|
||||
static __inline__ void
|
||||
acb_get_abs_ubound_arf(arf_t u, const acb_t z, long prec)
|
||||
{
|
||||
|
|
@ -660,6 +668,21 @@ void acb_root_newton(acb_t r, const acb_t a, long m, long index, long prec);
|
|||
void acb_root(acb_t r, const acb_t a, long m, long index, long prec);
|
||||
*/
|
||||
|
||||
/* TODO: document */
|
||||
static __inline__ int
|
||||
acb_is_finite(const acb_t x)
|
||||
{
|
||||
return arb_is_finite(acb_realref(x)) && arb_is_finite(acb_imagref(x));
|
||||
}
|
||||
|
||||
/* TODO: document */
|
||||
static __inline__ void
|
||||
acb_indeterminate(acb_t x)
|
||||
{
|
||||
arb_indeterminate(acb_realref(x));
|
||||
arb_indeterminate(acb_imagref(x));
|
||||
}
|
||||
|
||||
static __inline__ void
|
||||
_acb_vec_zero(acb_ptr A, long n)
|
||||
{
|
||||
|
|
|
|||
41
acb_poly.h
41
acb_poly.h
|
|
@ -507,25 +507,32 @@ void _acb_poly_rising_ui_series(acb_ptr res, acb_srcptr f, long flen, ulong r, l
|
|||
void acb_poly_rising_ui_series(acb_poly_t res, const acb_poly_t f, ulong r, long trunc, long prec);
|
||||
|
||||
/* TODO: document */
|
||||
/* series of c^(d+x) */
|
||||
static __inline__ void
|
||||
_acb_poly_pow_cpx(acb_ptr res, const acb_t c, const acb_t d, long trunc, long prec)
|
||||
void
|
||||
_acb_poly_acb_pow_cpx(acb_ptr w, const acb_t a, const acb_t b, long len, long prec)
|
||||
{
|
||||
long i;
|
||||
acb_t logc;
|
||||
|
||||
acb_init(logc);
|
||||
acb_log(logc, c, prec);
|
||||
acb_mul(res + 0, logc, d, prec);
|
||||
acb_exp(res + 0, res + 0, prec);
|
||||
|
||||
for (i = 1; i < trunc; i++)
|
||||
if (len == 1)
|
||||
{
|
||||
acb_mul(res + i, res + i - 1, logc, prec);
|
||||
acb_div_ui(res + i, res + i, i, prec);
|
||||
acb_pow(w, a, b, prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
acb_t log_a;
|
||||
long k;
|
||||
|
||||
acb_init(log_a);
|
||||
|
||||
acb_log(log_a, a, prec);
|
||||
acb_mul(w, log_a, b, prec);
|
||||
acb_exp(w, w, prec);
|
||||
|
||||
for (k = 1; k < len; k++)
|
||||
{
|
||||
acb_mul(w + k, w + k - 1, log_a, prec);
|
||||
acb_div_ui(w + k, w + k, k, prec);
|
||||
}
|
||||
|
||||
acb_clear(logc);
|
||||
acb_clear(log_a);
|
||||
}
|
||||
}
|
||||
|
||||
/* TODO: document */
|
||||
|
|
@ -551,6 +558,10 @@ void _acb_poly_zeta_series(acb_ptr res, acb_srcptr h, long hlen, const acb_t a,
|
|||
|
||||
void acb_poly_zeta_series(acb_poly_t res, const acb_poly_t f, const acb_t a, int deflate, long n, long prec);
|
||||
|
||||
void _acb_poly_polylog_cpx_zeta(acb_ptr w, const acb_t s, const acb_t z, long len, long prec);
|
||||
void _acb_poly_polylog_cpx_small(acb_ptr w, const acb_t s, const acb_t z, long len, long prec);
|
||||
void _acb_poly_polylog_cpx(acb_ptr w, const acb_t s, const acb_t z, long len, long prec);
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
|
|
|||
472
acb_poly/polylog_series.c
Normal file
472
acb_poly/polylog_series.c
Normal file
|
|
@ -0,0 +1,472 @@
|
|||
/*=============================================================================
|
||||
|
||||
This file is part of ARB.
|
||||
|
||||
ARB is free software; you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
ARB is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with ARB; if not, write to the Free Software
|
||||
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
|
||||
|
||||
=============================================================================*/
|
||||
/******************************************************************************
|
||||
|
||||
Copyright (C) 2014 Fredrik Johansson
|
||||
|
||||
******************************************************************************/
|
||||
|
||||
#include "acb_poly.h"
|
||||
|
||||
/* TODO: move and document helper functions */
|
||||
|
||||
/* TODO: implement efficiently */
|
||||
void
|
||||
mag_sqrt(mag_t y, const mag_t x)
|
||||
{
|
||||
if (mag_is_special(x))
|
||||
{
|
||||
mag_set(y, x);
|
||||
}
|
||||
else
|
||||
{
|
||||
arf_t t;
|
||||
arf_init(t);
|
||||
arf_set_mag(t, x);
|
||||
arf_sqrt(t, t, MAG_BITS, ARF_RND_UP);
|
||||
arf_get_mag(y, t);
|
||||
arf_clear(t);
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
acb_get_mag(mag_t u, const acb_t z)
|
||||
{
|
||||
if (arb_is_zero(acb_imagref(z)))
|
||||
{
|
||||
arb_get_mag(u, acb_realref(z));
|
||||
}
|
||||
else if (arb_is_zero(acb_realref(z)))
|
||||
{
|
||||
arb_get_mag(u, acb_imagref(z));
|
||||
}
|
||||
else
|
||||
{
|
||||
mag_t v;
|
||||
mag_init(v);
|
||||
|
||||
arb_get_mag(u, acb_realref(z));
|
||||
arb_get_mag(v, acb_imagref(z));
|
||||
|
||||
mag_mul(u, u, u);
|
||||
mag_addmul(u, v, v);
|
||||
mag_sqrt(u, u);
|
||||
|
||||
mag_clear(v);
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
mag_log_ui(mag_t t, ulong n)
|
||||
{
|
||||
if (n <= 1)
|
||||
{
|
||||
if (n == 1)
|
||||
mag_zero(t);
|
||||
else
|
||||
mag_inf(t);
|
||||
}
|
||||
else
|
||||
{
|
||||
mag_set_ui(t, n-1);
|
||||
mag_log1p(t, t);
|
||||
}
|
||||
}
|
||||
|
||||
/* TODO: needs reimplementing */
|
||||
void mag_sub_lower(mag_t z, const mag_t x, const mag_t y);
|
||||
|
||||
long
|
||||
arb_get_si_lower(const arb_t x)
|
||||
{
|
||||
arf_t t;
|
||||
long v;
|
||||
|
||||
arf_init(t);
|
||||
arf_set_mag(t, arb_radref(x));
|
||||
arf_sub(t, arb_midref(x), t, 2 * FLINT_BITS, ARF_RND_FLOOR);
|
||||
|
||||
v = arf_get_si(t, ARF_RND_FLOOR);
|
||||
|
||||
arf_clear(t);
|
||||
|
||||
return v;
|
||||
}
|
||||
|
||||
/* bound (1 + 1/m)^n, m > 0, n >= 0 */
|
||||
void
|
||||
mag_binpow_uiui(mag_t b, ulong m, ulong n)
|
||||
{
|
||||
mag_t t;
|
||||
mag_init(t);
|
||||
|
||||
/* bound by exp(n/m) <= 1 + (n/m) + (n/m)^2 */
|
||||
if (m > n)
|
||||
{
|
||||
mag_set_ui(t, n); /* x = n/m */
|
||||
mag_div_ui(t, t, m);
|
||||
|
||||
mag_mul(b, t, t); /* x^2 */
|
||||
mag_add(b, b, t); /* x */
|
||||
mag_one(t);
|
||||
mag_add(b, b, t); /* 1 */
|
||||
}
|
||||
else
|
||||
{
|
||||
mag_one(b);
|
||||
mag_div_ui(b, b, m);
|
||||
mag_one(t);
|
||||
mag_add(t, t, b);
|
||||
mag_pow_ui(b, t, n);
|
||||
}
|
||||
|
||||
mag_clear(t);
|
||||
}
|
||||
|
||||
void
|
||||
polylog_remainder_bound(mag_t u, const mag_t z, long sigma, ulong d, ulong N)
|
||||
{
|
||||
mag_t TN, UN, t;
|
||||
|
||||
if (N < 2)
|
||||
{
|
||||
mag_inf(u);
|
||||
return;
|
||||
}
|
||||
|
||||
mag_init(TN);
|
||||
mag_init(UN);
|
||||
mag_init(t);
|
||||
|
||||
if (mag_cmp_2exp_si(z, 0) >= 0)
|
||||
{
|
||||
mag_inf(u);
|
||||
}
|
||||
else
|
||||
{
|
||||
/* Bound T(N) */
|
||||
mag_pow_ui(TN, z, N);
|
||||
|
||||
/* multiply by log(N)^d */
|
||||
if (d > 0)
|
||||
{
|
||||
mag_log_ui(t, N);
|
||||
mag_pow_ui(t, t, d);
|
||||
mag_mul(TN, TN, t);
|
||||
}
|
||||
|
||||
/* multiply by 1/k^s */
|
||||
if (sigma > 0)
|
||||
{
|
||||
mag_set_ui_lower(t, N);
|
||||
mag_pow_ui_lower(t, t, sigma);
|
||||
mag_div(TN, TN, t);
|
||||
}
|
||||
else if (sigma < 0)
|
||||
{
|
||||
mag_set_ui(t, N);
|
||||
mag_pow_ui(t, t, -sigma);
|
||||
mag_mul(TN, TN, t);
|
||||
}
|
||||
|
||||
/* Bound U(N) */
|
||||
mag_set(UN, z);
|
||||
|
||||
/* multiply by (1 + 1/N)**S */
|
||||
if (sigma < 0)
|
||||
{
|
||||
mag_binpow_uiui(t, N, -sigma);
|
||||
mag_mul(UN, UN, t);
|
||||
}
|
||||
|
||||
/* multiply by (1 + 1/(N log(N)))^d */
|
||||
if (d > 0)
|
||||
{
|
||||
ulong nl;
|
||||
|
||||
/* rounds down */
|
||||
nl = mag_d_log_lower_bound(N) * N * (1 - 1e-13);
|
||||
|
||||
mag_binpow_uiui(t, nl, d);
|
||||
mag_mul(UN, UN, t);
|
||||
}
|
||||
|
||||
/* T(N) / (1 - U(N)) */
|
||||
if (mag_cmp_2exp_si(UN, 0) >= 0)
|
||||
{
|
||||
mag_inf(u);
|
||||
}
|
||||
else
|
||||
{
|
||||
mag_one(t);
|
||||
mag_sub_lower(t, t, UN);
|
||||
mag_div(u, TN, t);
|
||||
}
|
||||
}
|
||||
|
||||
mag_clear(TN);
|
||||
mag_clear(UN);
|
||||
mag_clear(t);
|
||||
}
|
||||
|
||||
long
|
||||
polylog_choose_terms(mag_t err, long sigma, const mag_t z, long d, long prec)
|
||||
{
|
||||
long N;
|
||||
|
||||
for (N = 3; ; N = FLINT_MAX(N+3, N*1.1))
|
||||
{
|
||||
polylog_remainder_bound(err, z, sigma, d, N);
|
||||
|
||||
/* TODO: do something else when |Li_s(z)| is very small/very large? */
|
||||
if (mag_cmp_2exp_si(err, -prec) < 0)
|
||||
break;
|
||||
|
||||
if (N > 100 * prec)
|
||||
{
|
||||
N = 3;
|
||||
mag_inf(err);
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
return N;
|
||||
}
|
||||
|
||||
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
|
||||
{
|
||||
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;
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
_acb_poly_polylog_cpx_zeta(acb_ptr w, const acb_t s, const acb_t z, long len, long prec)
|
||||
{
|
||||
acb_ptr e1, e2, z1, z2, e1z1, e2z2;
|
||||
acb_t t, u, v;
|
||||
long k, len2;
|
||||
int deflate_zeta, deflate_gamma;
|
||||
|
||||
if (!acb_is_finite(s) || !acb_is_finite(z))
|
||||
{
|
||||
_acb_vec_indeterminate(w, len);
|
||||
return;
|
||||
}
|
||||
|
||||
acb_init(t);
|
||||
acb_init(u);
|
||||
acb_init(v);
|
||||
|
||||
/* v = 1-s */
|
||||
acb_one(v);
|
||||
acb_sub(v, v, s, prec);
|
||||
|
||||
/* pole of zeta */
|
||||
deflate_zeta = acb_is_one(v);
|
||||
|
||||
/* poles of gamma at nonpositive integer v */
|
||||
deflate_gamma = (arb_is_zero(acb_imagref(v)) &&
|
||||
arb_is_int(acb_realref(v)) &&
|
||||
arf_sgn(arb_midref(acb_realref(v))) <= 0);
|
||||
|
||||
len2 = len + deflate_gamma;
|
||||
|
||||
e1 = _acb_vec_init(len + 1);
|
||||
e2 = _acb_vec_init(len + 1);
|
||||
z1 = _acb_vec_init(len + 1);
|
||||
z2 = _acb_vec_init(len + 1);
|
||||
e1z1 = _acb_vec_init(len + 1);
|
||||
e2z2 = _acb_vec_init(len + 1);
|
||||
|
||||
/* u = log(-z)/(pi*i) */
|
||||
acb_neg(t, z);
|
||||
acb_log(t, t, prec);
|
||||
acb_const_pi(u, prec);
|
||||
acb_mul_onei(u, u);
|
||||
acb_div(u, t, u, prec);
|
||||
|
||||
/* z1 = zeta(v+x, 1/2 + log(-z)/(2*pi*i)) */
|
||||
acb_one(t);
|
||||
acb_add(t, t, u, prec);
|
||||
acb_mul_2exp_si(t, t, -1);
|
||||
_acb_poly_zeta_cpx_series(z1, v, t, deflate_zeta, len2, prec);
|
||||
|
||||
/* z2 = zeta(v+x, 1/2 - log(-z)/(2*pi*i)) */
|
||||
acb_one(t);
|
||||
acb_sub(t, t, u, prec);
|
||||
acb_mul_2exp_si(t, t, -1);
|
||||
_acb_poly_zeta_cpx_series(z2, v, t, deflate_zeta, len2, prec);
|
||||
|
||||
/* e1 = (i/(2pi))^(v+x) */
|
||||
acb_onei(t);
|
||||
acb_const_pi(u, prec);
|
||||
acb_div(t, t, u, prec);
|
||||
acb_mul_2exp_si(t, t, -1);
|
||||
_acb_poly_acb_pow_cpx(e1, t, v, len + (deflate_zeta || deflate_gamma), prec);
|
||||
|
||||
/* e2 = (1/(2 pi i))^(v+x) */
|
||||
acb_conj(t, t);
|
||||
_acb_poly_acb_pow_cpx(e2, t, v, len + (deflate_zeta || deflate_gamma), prec);
|
||||
|
||||
_acb_poly_mullow(e1z1, e1, len2, z1, len2, len2, prec);
|
||||
_acb_poly_mullow(e2z2, e2, len2, z2, len2, len2, prec);
|
||||
_acb_vec_add(z1, e1z1, e2z2, len2, prec);
|
||||
|
||||
if (deflate_gamma)
|
||||
{
|
||||
/* gamma(v+x) = pi/sin(pi(v+x)) * 1/gamma(1-v-x) */
|
||||
|
||||
/* TODO: write a csc function? */
|
||||
acb_zero(e1);
|
||||
acb_const_pi(e1 + 1, prec);
|
||||
acb_mul_2exp_si(e2, v, -1);
|
||||
if (!arb_is_int(acb_realref(e2)))
|
||||
acb_neg(e1 + 1, e1 + 1);
|
||||
_acb_poly_sin_series(e2, e1, 2, len2, prec);
|
||||
_acb_poly_inv_series(e1, e2 + 1, len, len, prec);
|
||||
acb_const_pi(e2, prec);
|
||||
_acb_vec_scalar_mul(e1, e1, len, e2, prec);
|
||||
|
||||
acb_set(z2, s);
|
||||
acb_set_si(z2 + 1, -1);
|
||||
_acb_poly_rgamma_series(e2, z2, 2, len, prec);
|
||||
_acb_poly_mullow(z2, e1, len, e2, len, len, prec);
|
||||
|
||||
_acb_poly_mullow(w, z1 + 1, len, z2, len, len, prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
if (deflate_zeta)
|
||||
{
|
||||
for (k = 0; k < len; k++)
|
||||
{
|
||||
arb_mul_2exp_si(acb_realref(e1 + k + 1), acb_realref(e1 + k + 1), 1);
|
||||
arb_add(acb_realref(z1 + k), acb_realref(z1 + k), acb_realref(e1 + k + 1), prec);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/* gamma(v+x) */
|
||||
acb_set(e1, v);
|
||||
if (len > 1)
|
||||
acb_one(e1 + 1);
|
||||
_acb_poly_gamma_series(z2, e1, FLINT_MIN(len, 2), len, prec);
|
||||
|
||||
_acb_poly_mullow(w, z2, len, z1, len, len, prec);
|
||||
}
|
||||
|
||||
/* correct signs (from s -> 1-s) */
|
||||
for (k = 1; k < len; k += 2)
|
||||
acb_neg(w + k, w + k);
|
||||
|
||||
_acb_vec_clear(e1, len + 1);
|
||||
_acb_vec_clear(e2, len + 1);
|
||||
_acb_vec_clear(z1, len + 1);
|
||||
_acb_vec_clear(z2, len + 1);
|
||||
_acb_vec_clear(e1z1, len + 1);
|
||||
_acb_vec_clear(e2z2, len + 1);
|
||||
|
||||
acb_clear(t);
|
||||
acb_clear(u);
|
||||
acb_clear(v);
|
||||
}
|
||||
|
||||
void
|
||||
_acb_poly_polylog_cpx_small(acb_ptr w, const acb_t s, const acb_t z, long len, long prec)
|
||||
{
|
||||
long k, N, sigma;
|
||||
int is_real;
|
||||
mag_t zmag, err, errf;
|
||||
acb_t a;
|
||||
|
||||
acb_init(a);
|
||||
mag_init(zmag);
|
||||
mag_init(err);
|
||||
mag_init(errf);
|
||||
|
||||
is_real = polylog_is_real(s, z);
|
||||
acb_get_mag(zmag, z);
|
||||
sigma = arb_get_si_lower(acb_realref(s));
|
||||
|
||||
N = polylog_choose_terms(err, sigma, zmag, len - 1, prec);
|
||||
|
||||
/* TODO: allow threading */
|
||||
acb_one(a);
|
||||
_acb_poly_powsum_series_naive(w, s, a, z, N - 1, len, prec);
|
||||
_acb_vec_scalar_mul(w, w, len, z, prec);
|
||||
|
||||
for (k = 0; k < len; k++)
|
||||
{
|
||||
polylog_remainder_bound(err, zmag, sigma, k, N);
|
||||
mag_rfac_ui(errf, k);
|
||||
mag_mul(err, err, errf);
|
||||
|
||||
if (is_real)
|
||||
arb_add_error_mag(acb_realref(w + k), err);
|
||||
else
|
||||
acb_add_error_mag(w + k, err);
|
||||
}
|
||||
|
||||
acb_clear(a);
|
||||
mag_clear(zmag);
|
||||
mag_clear(err);
|
||||
mag_clear(errf);
|
||||
}
|
||||
|
||||
void
|
||||
_acb_poly_polylog_cpx(acb_ptr w, const acb_t s, const acb_t z, long len, long prec)
|
||||
{
|
||||
mag_t zmag;
|
||||
|
||||
mag_init(zmag);
|
||||
acb_get_mag(zmag, z);
|
||||
|
||||
if (mag_cmp_2exp_si(zmag, -1) < 0)
|
||||
_acb_poly_polylog_cpx_small(w, s, z, len, prec);
|
||||
else
|
||||
_acb_poly_polylog_cpx_zeta(w, s, z, len, prec);
|
||||
|
||||
mag_clear(zmag);
|
||||
}
|
||||
|
||||
412
acb_poly/test/t-polylog_cpx.c
Normal file
412
acb_poly/test/t-polylog_cpx.c
Normal file
|
|
@ -0,0 +1,412 @@
|
|||
/*=============================================================================
|
||||
|
||||
This file is part of ARB.
|
||||
|
||||
ARB is free software; you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
ARB is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with ARB; if not, write to the Free Software
|
||||
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
|
||||
|
||||
=============================================================================*/
|
||||
/******************************************************************************
|
||||
|
||||
Copyright (C) 2013 Fredrik Johansson
|
||||
|
||||
******************************************************************************/
|
||||
|
||||
#include "acb_poly.h"
|
||||
|
||||
#define EPS 1e-13
|
||||
#define NUM_DERIVS 3
|
||||
#define NUM_TESTS 56
|
||||
|
||||
const double polylog_testdata[NUM_TESTS][10] = {
|
||||
{-2.0, 0.0, 0.0, 0.0,
|
||||
0.0, 0.0,
|
||||
0.0, 0.0,
|
||||
0.0, 0.0},
|
||||
{-2.0, 0.0, -0.5, 0.0,
|
||||
-0.0740740740740740741, 0.0,
|
||||
-0.158799035455864473, 0.0,
|
||||
0.00278425812292977954, 0.0},
|
||||
{-2.0, 0.0, 0.99609375, 0.0,
|
||||
33358080.0, 0.0,
|
||||
-215693527.719090182, 0.0,
|
||||
703924885.809675425, 0.0},
|
||||
{-2.0, 0.0, 1.00390625, 0.0,
|
||||
-33751296.0, 0.0,
|
||||
218367905.967311029, -106032823.56283664,
|
||||
-546519337.1290267, 686023009.262281684},
|
||||
{-2.0, 0.0, -4.25, 0.0,
|
||||
0.0954540546377281071, 0.0,
|
||||
-0.0297817299750221174, 0.0,
|
||||
-0.196648773789941696, 0.0},
|
||||
{-2.0, 0.0, 3.5, 0.0,
|
||||
-1.008, 0.0,
|
||||
0.691243661008350802, -3.19575256895409845,
|
||||
4.53340070033764442, 2.22882281653983385},
|
||||
{-2.0, 0.0, 0.0, -1.0,
|
||||
0.0, 0.5,
|
||||
-0.852556797635011582, -0.25189536315108886,
|
||||
0.345626503121598963, -0.46292146649778986},
|
||||
{-2.0, 0.0, 0.125, 0.125,
|
||||
0.071424, 0.280832,
|
||||
0.0646851488459364134, -0.118051747326526959,
|
||||
-0.0396751348737248634, 0.0452067417426917285},
|
||||
{-1.0, 0.0, 0.0, 0.0,
|
||||
0.0, 0.0,
|
||||
0.0, 0.0,
|
||||
0.0, 0.0},
|
||||
{-1.0, 0.0, -0.5, 0.0,
|
||||
-0.222222222222222222, 0.0,
|
||||
-0.131883201632587537, 0.0,
|
||||
0.0199225628389378777, 0.0},
|
||||
{-1.0, 0.0, 0.99609375, 0.0,
|
||||
65280.0, 0.0,
|
||||
-389461.499432291984, 0.0,
|
||||
1182814.77437627002, 0.0},
|
||||
{-1.0, 0.0, 1.00390625, 0.0,
|
||||
65792.0, 0.0,
|
||||
-392773.095775796478, 206691.925664168516,
|
||||
868954.597422132011, -1233932.55215831045},
|
||||
{-1.0, 0.0, -4.25, 0.0,
|
||||
-0.154195011337868481, 0.0,
|
||||
-0.474824989764482489, 0.0,
|
||||
-0.21428706813151403, 0.0},
|
||||
{-1.0, 0.0, 3.5, 0.0,
|
||||
0.56, 0.0,
|
||||
-0.322735710701998933, 2.00176023742981729,
|
||||
-3.09698193717149541, -0.395213225234415973},
|
||||
{-1.0, 0.0, 0.0, -1.0,
|
||||
-0.5, 0.0,
|
||||
-0.183855151549436048, -0.58312180806163756,
|
||||
0.254077251327014729, 0.0351429519611278982},
|
||||
{-1.0, 0.0, 0.125, 0.125,
|
||||
0.1088, 0.1984,
|
||||
0.0192018441183911, -0.0548368329876085499,
|
||||
-0.0115465156713130216, 0.020954643261003461},
|
||||
{0.0, 0.0, 0.0, 0.0,
|
||||
0.0, 0.0,
|
||||
0.0, 0.0,
|
||||
0.0, 0.0},
|
||||
{0.0, 0.0, -0.5, 0.0,
|
||||
-0.333333333333333333, 0.0,
|
||||
-0.0904116649848191551, 0.0,
|
||||
0.0199028637728058217, 0.0},
|
||||
{0.0, 0.0, 0.99609375, 0.0,
|
||||
255.0, 0.0,
|
||||
-1269.73106874463507, 0.0,
|
||||
3359.60296533314103, 0.0},
|
||||
{0.0, 0.0, 1.00390625, 0.0,
|
||||
-257.0, 0.0,
|
||||
1273.86116119868176, -805.817494984366455,
|
||||
-2113.9649224353639, 4004.84182305521042},
|
||||
{0.0, 0.0, -4.25, 0.0,
|
||||
-0.809523809523809524, 0.0,
|
||||
-0.795769271870623909, 0.0,
|
||||
-0.0965372637237445457, 0.0},
|
||||
{0.0, 0.0, 3.5, 0.0,
|
||||
-1.4, 0.0,
|
||||
-1.78816683065098841, -2.50773109725857055,
|
||||
1.83790031889197745, -2.01262260402527512},
|
||||
{0.0, 0.0, 0.0, -1.0,
|
||||
-0.5, -0.5,
|
||||
0.120782237635245222, -0.391594392706836776,
|
||||
0.0669079578849868665, 0.116416441749488305},
|
||||
{0.0, 0.0, 0.125, 0.125,
|
||||
0.12, 0.16,
|
||||
0.00583206536388767829, -0.0256940749384562616,
|
||||
-0.00344424777009813297, 0.00965046958017889366},
|
||||
{-0.5, 0.0, 0.0, 0.0,
|
||||
0.0, 0.0,
|
||||
0.0, 0.0,
|
||||
0.0, 0.0},
|
||||
{-0.5, 0.0, -0.5, 0.0,
|
||||
-0.283012810746506023, 0.0,
|
||||
-0.111080563146094489, 0.0,
|
||||
0.0211430885590715694, 0.0},
|
||||
{-0.5, 0.0, 0.99609375, 0.0,
|
||||
3619.14124113924805, 0.0,
|
||||
-20195.2825251289493, 0.0,
|
||||
58032.3050555790807, 0.0},
|
||||
{-0.5, 0.0, 1.00390625, 0.0,
|
||||
-0.207985517845521836, 3640.61848390552728,
|
||||
-11437.7014354502895, -20327.8206114649109,
|
||||
63861.4335397379033, 40487.2831296029499},
|
||||
{-0.5, 0.0, -4.25, 0.0,
|
||||
-0.441427843410865187, 0.0,
|
||||
-0.665497904840929966, 0.0,
|
||||
-0.162969337925229287, 0.0},
|
||||
{-0.5, 0.0, 3.5, 0.0,
|
||||
-0.232043460812781422, 0.63203565652521718,
|
||||
-2.44367789279529917, 0.119367210264842707,
|
||||
-0.763914228371074801, -2.81228485141461073},
|
||||
{-0.5, 0.0, 0.0, -1.0,
|
||||
-0.537549381115898973, -0.27517974122882025,
|
||||
0.0158915835170861709, -0.505981074388456133,
|
||||
0.147999289996985447, 0.105632732120539676},
|
||||
{-0.5, 0.0, 0.125, 0.125,
|
||||
0.116021225799463628, 0.175602658223910358,
|
||||
0.0105539353881710306, -0.0374742285737420805,
|
||||
-0.0062875342272019156, 0.0142081398192247589},
|
||||
{1.0, 0.0, 0.0, 0.0,
|
||||
0.0, 0.0,
|
||||
0.0, 0.0,
|
||||
0.0, 0.0},
|
||||
{1.0, 0.0, -0.5, 0.0,
|
||||
-0.405465108108164382, 0.0,
|
||||
-0.0556552111040376179, 0.0,
|
||||
0.0145754167941056437, 0.0},
|
||||
{1.0, 0.0, 0.99609375, 0.0,
|
||||
5.54517744447956248, 0.0,
|
||||
-13.0766609768056118, 0.0,
|
||||
24.0940545588450487, 0.0},
|
||||
{1.0, 0.0, 1.00390625, 0.0,
|
||||
5.54517744447956248, -3.14159265358979324,
|
||||
-8.16844417312495553, 15.6134381896777968,
|
||||
-0.38789687903338657, -36.2148484363837508},
|
||||
{1.0, 0.0, -4.25, 0.0,
|
||||
-1.65822807660353238, 0.0,
|
||||
-0.861140718490863904, 0.0,
|
||||
0.0240108114276114004, 0.0},
|
||||
{1.0, 0.0, 3.5, 0.0,
|
||||
-0.916290731874155065, -3.14159265358979324,
|
||||
2.57284574234941582, -2.52133906787958123,
|
||||
1.44241441233798233, 1.57208443290129947},
|
||||
{1.0, 0.0, 0.0, -1.0,
|
||||
-0.346573590279972655, -0.78539816339744831,
|
||||
0.160292055087885226, -0.192901316796912429,
|
||||
-0.0115066749133700615, 0.0770708622146679417},
|
||||
{1.0, 0.0, 0.125, 0.125,
|
||||
0.123430038965762899, 0.141897054604163923,
|
||||
0.00180689203856968544, -0.0122114794810442984,
|
||||
-0.00105066979284378462, 0.00449853805632375062},
|
||||
{2.0, 0.0, 0.0, 0.0,
|
||||
0.0, 0.0,
|
||||
0.0, 0.0,
|
||||
0.0, 0.0},
|
||||
{2.0, 0.0, -0.5, 0.0,
|
||||
-0.448414206923646202, 0.0,
|
||||
-0.0320237060779612787, 0.0,
|
||||
0.0092536571135294828, 0.0},
|
||||
{2.0, 0.0, 0.99609375, 0.0,
|
||||
1.61932072927810507, 0.0,
|
||||
-0.863010139919679163, 0.0,
|
||||
0.825515949454683347, 0.0},
|
||||
{2.0, 0.0, 1.00390625, 0.0,
|
||||
1.6704551616562229, -0.0122479400888173039,
|
||||
-0.992661633138610541, 0.0731191212424626706,
|
||||
1.04852574935488925, -0.214307793003452771},
|
||||
{2.0, 0.0, -4.25, 0.0,
|
||||
-2.46898738747224279, 0.0,
|
||||
-0.739287387997829118, 0.0,
|
||||
0.0874470644522297429, 0.0},
|
||||
{2.0, 0.0, 3.5, 0.0,
|
||||
2.1959348115757808, -3.93567093851438969,
|
||||
2.87640729249480637, 0.777030723254221356,
|
||||
-0.81939570865486996, 1.19241843763256774},
|
||||
{2.0, 0.0, 0.0, -1.0,
|
||||
-0.205616758356028305, -0.915965594177219015,
|
||||
0.117193531789480469, -0.0815807361165927951,
|
||||
-0.0255408233794369879, 0.0372076062178391175},
|
||||
{2.0, 0.0, 0.125, 0.125,
|
||||
0.124500148439393858, 0.133240721203487057,
|
||||
0.000569113907620539642, -0.00588251562210034232,
|
||||
-0.000326714635669985778, 0.00213086661076891839},
|
||||
{-2.0, 2.0, 0.0, 0.0,
|
||||
0.0, 0.0,
|
||||
0.0, 0.0,
|
||||
0.0, 0.0},
|
||||
{-2.0, 2.0, -0.5, 0.0,
|
||||
-0.132203419506655614, -0.390689361106274866,
|
||||
-0.256343261725123592, 0.106175177622871437,
|
||||
0.0749290443999118563, 0.0319207027895210173},
|
||||
{-2.0, 2.0, 0.99609375, 0.0,
|
||||
13841369.6929212372, -8339597.62269406988,
|
||||
-87251051.2148652244, 65225350.0315463319,
|
||||
273254167.909661771, -247710863.316563213},
|
||||
{-2.0, 2.0, 1.00390625, 0.0,
|
||||
-26028.6661639721745, 15961.1559569935978,
|
||||
113847.512122984734, -206358.723605581203,
|
||||
6578.40068463757013, 908971.828153761815},
|
||||
{-2.0, 2.0, -4.25, 0.0,
|
||||
1.26205464628714387, 0.274313692011971484,
|
||||
0.442242006131113713, -1.61799252235272836,
|
||||
-0.927785295033312282, -0.172205994085466244},
|
||||
{-2.0, 2.0, 3.5, 0.0,
|
||||
0.0870998538268457367, 0.0980699787750225838,
|
||||
0.152937181587789832, -0.0211963252018165386,
|
||||
0.0431921445114569281, -0.0695183434413184998},
|
||||
{-2.0, 2.0, 0.0, -1.0,
|
||||
0.199312975007554535, -0.0775592701795263651,
|
||||
-0.0210851901426345675, -0.167326332453265841,
|
||||
-0.0632997708310078762, -0.0392019115527122943},
|
||||
{-2.0, 2.0, 0.125, 0.125,
|
||||
0.315234332693229892, 0.165559827416958489,
|
||||
-0.16530664867351669, -0.0392852937947218337,
|
||||
0.0780220225006149953, 0.0215283250955342792},
|
||||
};
|
||||
|
||||
int main()
|
||||
{
|
||||
long iter;
|
||||
flint_rand_t state;
|
||||
|
||||
printf("polylog_cpx....");
|
||||
fflush(stdout);
|
||||
|
||||
flint_randinit(state);
|
||||
|
||||
/* check particular values against table */
|
||||
{
|
||||
acb_t s, z, t;
|
||||
acb_ptr w1, w2;
|
||||
long i, j, prec;
|
||||
|
||||
acb_init(s);
|
||||
acb_init(z);
|
||||
acb_init(t);
|
||||
w1 = _acb_vec_init(NUM_DERIVS);
|
||||
w2 = _acb_vec_init(NUM_DERIVS);
|
||||
|
||||
for (prec = 32; prec <= 512; prec *= 4)
|
||||
{
|
||||
for (i = 0; i < NUM_TESTS; i++)
|
||||
{
|
||||
acb_zero(s);
|
||||
arf_set_d(arb_midref(acb_realref(s)), polylog_testdata[i][0]);
|
||||
arf_set_d(arb_midref(acb_imagref(s)), polylog_testdata[i][1]);
|
||||
|
||||
acb_zero(z);
|
||||
arf_set_d(arb_midref(acb_realref(z)), polylog_testdata[i][2]);
|
||||
arf_set_d(arb_midref(acb_imagref(z)), polylog_testdata[i][3]);
|
||||
|
||||
_acb_poly_polylog_cpx_small(w1, s, z, NUM_DERIVS, prec);
|
||||
_acb_poly_polylog_cpx_zeta(w2, s, z, NUM_DERIVS, prec);
|
||||
|
||||
for (j = 0; j < NUM_DERIVS; j++)
|
||||
{
|
||||
arf_set_d(arb_midref(acb_realref(t)), polylog_testdata[i][4+2*j]);
|
||||
mag_set_d(arb_radref(acb_realref(t)), fabs(polylog_testdata[i][4+2*j]) * EPS);
|
||||
arf_set_d(arb_midref(acb_imagref(t)), polylog_testdata[i][4+2*j+1]);
|
||||
mag_set_d(arb_radref(acb_imagref(t)), fabs(polylog_testdata[i][4+2*j+1]) * EPS);
|
||||
|
||||
if (!acb_overlaps(w1 + j, t) || !acb_overlaps(w2 + j, t)
|
||||
|| !acb_overlaps(w1 + j, w2 + j))
|
||||
{
|
||||
printf("FAIL\n\n");
|
||||
printf("s = "); acb_printd(s, 15); printf("\n\n");
|
||||
printf("z = "); acb_printd(z, 15); printf("\n\n");
|
||||
printf("t = "); acb_printd(t, 15); printf("\n\n");
|
||||
printf("w1 = "); acb_printd(w1 + j, 15); printf("\n\n");
|
||||
printf("w2 = "); acb_printd(w2 + j, 15); printf("\n\n");
|
||||
abort();
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
_acb_vec_clear(w1, NUM_DERIVS);
|
||||
_acb_vec_clear(w2, NUM_DERIVS);
|
||||
acb_clear(s);
|
||||
acb_clear(z);
|
||||
acb_clear(t);
|
||||
}
|
||||
|
||||
for (iter = 0; iter < 250; iter++)
|
||||
{
|
||||
acb_t s, z;
|
||||
acb_ptr w1, w2;
|
||||
long i, len1, len2, prec1, prec2;
|
||||
|
||||
acb_init(s);
|
||||
acb_init(z);
|
||||
|
||||
if (n_randint(state, 2))
|
||||
acb_randtest(s, state, 1 + n_randint(state, 300), 3);
|
||||
else
|
||||
acb_set_si(s, n_randint(state, 20) - 10);
|
||||
|
||||
switch (n_randint(state, 3))
|
||||
{
|
||||
case 0:
|
||||
acb_randtest(z, state, 1 + n_randint(state, 300), 10);
|
||||
break;
|
||||
case 1:
|
||||
arb_randtest(acb_realref(z), state, 1 + n_randint(state, 300), 10);
|
||||
break;
|
||||
case 2:
|
||||
acb_one(z);
|
||||
break;
|
||||
}
|
||||
|
||||
prec1 = 2 + n_randint(state, 300);
|
||||
prec2 = prec1 + 30;
|
||||
len1 = 1 + n_randint(state, 20);
|
||||
len2 = 1 + n_randint(state, 20);
|
||||
|
||||
w1 = _acb_vec_init(len1);
|
||||
w2 = _acb_vec_init(len2);
|
||||
|
||||
switch (n_randint(state, 5))
|
||||
{
|
||||
case 0:
|
||||
_acb_poly_polylog_cpx_small(w1, s, z, len1, prec1);
|
||||
break;
|
||||
case 1:
|
||||
_acb_poly_polylog_cpx_zeta(w1, s, z, len1, prec1);
|
||||
break;
|
||||
default:
|
||||
_acb_poly_polylog_cpx(w1, s, z, len1, prec1);
|
||||
break;
|
||||
}
|
||||
|
||||
switch (n_randint(state, 5))
|
||||
{
|
||||
case 0:
|
||||
_acb_poly_polylog_cpx_small(w2, s, z, len2, prec2);
|
||||
break;
|
||||
case 1:
|
||||
_acb_poly_polylog_cpx_zeta(w2, s, z, len2, prec2);
|
||||
break;
|
||||
default:
|
||||
_acb_poly_polylog_cpx(w2, s, z, len2, prec2);
|
||||
break;
|
||||
}
|
||||
|
||||
for (i = 0; i < FLINT_MIN(len1, len2); i++)
|
||||
{
|
||||
if (!acb_overlaps(w1 + i, w2 + i))
|
||||
{
|
||||
printf("FAIL: overlap\n\n");
|
||||
printf("iter = %ld\n", iter);
|
||||
printf("len1 = %ld, len2 = %ld, i = %ld\n\n", len1, len2, i);
|
||||
printf("s = "); acb_printd(s, prec1 / 3.33); printf("\n\n");
|
||||
printf("z = "); acb_printd(z, prec1 / 3.33); printf("\n\n");
|
||||
printf("w1 = "); acb_printd(w1 + i, prec1 / 3.33); printf("\n\n");
|
||||
printf("w2 = "); acb_printd(w2 + i, prec2 / 3.33); printf("\n\n");
|
||||
abort();
|
||||
}
|
||||
}
|
||||
|
||||
acb_clear(s);
|
||||
acb_clear(z);
|
||||
_acb_vec_clear(w1, len1);
|
||||
_acb_vec_clear(w2, len2);
|
||||
}
|
||||
|
||||
flint_randclear(state);
|
||||
flint_cleanup();
|
||||
printf("PASS\n");
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
|
||||
|
|
@ -36,6 +36,12 @@ _acb_poly_zeta_cpx_series(acb_ptr z, const acb_t s, const acb_t a, int deflate,
|
|||
if (d < 1)
|
||||
return;
|
||||
|
||||
if (!acb_is_finite(s) || !acb_is_finite(a))
|
||||
{
|
||||
_acb_vec_indeterminate(z, d);
|
||||
return;
|
||||
}
|
||||
|
||||
arf_init(bound);
|
||||
vb = _arb_vec_init(d);
|
||||
|
||||
|
|
|
|||
7
arb.h
7
arb.h
|
|
@ -489,6 +489,13 @@ void arb_add_error_2exp_fmpz(arb_t x, const fmpz_t err);
|
|||
|
||||
void arb_add_error(arb_t x, const arb_t error);
|
||||
|
||||
/* TODO: document */
|
||||
static __inline__ void
|
||||
arb_add_error_mag(arb_t x, const mag_t err)
|
||||
{
|
||||
mag_add(arb_radref(x), arb_radref(x), err);
|
||||
}
|
||||
|
||||
int arb_contains_arf(const arb_t x, const arf_t y);
|
||||
|
||||
int arb_contains_fmpq(const arb_t x, const fmpq_t y);
|
||||
|
|
|
|||
10
arb/gamma.c
10
arb/gamma.c
|
|
@ -30,15 +30,7 @@
|
|||
#include "hypgeom.h"
|
||||
|
||||
/* TODO: move helper functions to the right modules */
|
||||
void
|
||||
acb_get_mag(mag_t z, const acb_t x)
|
||||
{
|
||||
arf_t t;
|
||||
arf_init(t);
|
||||
acb_get_abs_ubound_arf(t, x, MAG_BITS);
|
||||
arf_get_mag(z, t);
|
||||
arf_clear(t);
|
||||
}
|
||||
void acb_get_mag(mag_t z, const acb_t x);
|
||||
|
||||
void
|
||||
acb_get_mag_lower(mag_t z, const acb_t x)
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue