mirror of
https://github.com/vale981/arb
synced 2025-03-05 09:21:38 -05:00
factor out mag_polylog_tail
This commit is contained in:
Fredrik Johansson
parent
dd7cdca99f
commit
a22075c52e
6 changed files with 298 additions and 134 deletions
|
|
@ -43,92 +43,6 @@ arb_get_si_lower(const arb_t x)
|
|||
return v;
|
||||
}
|
||||
|
||||
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)
|
||||
{
|
||||
|
|
@ -136,7 +50,7 @@ polylog_choose_terms(mag_t err, long sigma, const mag_t z, long d, long prec)
|
|||
|
||||
for (N = 3; ; N = FLINT_MAX(N+3, N*1.1))
|
||||
{
|
||||
polylog_remainder_bound(err, z, sigma, d, N);
|
||||
mag_polylog_tail(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)
|
||||
|
|
@ -341,7 +255,7 @@ _acb_poly_polylog_cpx_small(acb_ptr w, const acb_t s, const acb_t z, long len, l
|
|||
|
||||
for (k = 0; k < len; k++)
|
||||
{
|
||||
polylog_remainder_bound(err, zmag, sigma, k, N);
|
||||
mag_polylog_tail(err, zmag, sigma, k, N);
|
||||
mag_rfac_ui(errf, k);
|
||||
mag_mul(err, err, errf);
|
||||
|
||||
|
|
|
|||
|
|
@ -327,3 +327,52 @@ Special functions
|
|||
Sets *z* to an upper bound for `|B_n| / n!` where `B_n` denotes
|
||||
a Bernoulli number.
|
||||
|
||||
.. function:: void mag_polylog_tail(mag_t u, const mag_t z, long s, ulong d, ulong N)
|
||||
|
||||
Sets *u* to an upper bound for
|
||||
|
||||
.. math ::
|
||||
|
||||
\sum_{k=N}^{\infty} \frac{z^k \log^d(k)}{k^s}.
|
||||
|
||||
Note: in applications where `s` in this formula may be
|
||||
real or complex, the user can simply
|
||||
substitute any convenient integer `s'` such that `s' \le \operatorname{Re}(s)`.
|
||||
|
||||
Denote the terms by `T(k)`. We pick a nonincreasing function `U(k)` such that
|
||||
|
||||
.. math ::
|
||||
|
||||
\frac{T(k+1)}{T(k)} = z \left(\frac{k}{k+1}\right)^s
|
||||
\left( \frac{\log(k+1)}{\log(k)} \right)^d \le U(k).
|
||||
|
||||
Then, as soon as `U(N) < 1`,
|
||||
|
||||
.. math ::
|
||||
|
||||
\sum_{k=N}^{\infty} T(k)
|
||||
\le T(N) \sum_{k=0}^{\infty} U(N)^k = \frac{T(N)}{1 - U(N)}.
|
||||
|
||||
In particular, we take
|
||||
|
||||
.. math ::
|
||||
|
||||
U(k) = z \; B(k, \max(0, -s)) \; B(k \log(k), d)
|
||||
|
||||
where `B(m,n) = (1 + 1/m)^n`. This follows from the bounds
|
||||
|
||||
.. math ::
|
||||
|
||||
\left(\frac{k}{k+1}\right)^{s}
|
||||
\le \begin{cases}
|
||||
1 & \text{if } s \ge 0 \\
|
||||
(1 + 1/k)^{-s} & \text{if } s < 0.
|
||||
\end{cases}
|
||||
|
||||
and
|
||||
|
||||
.. math ::
|
||||
|
||||
\left( \frac{\log(k+1)}{\log(k)} \right)^d \le
|
||||
\left(1 + \frac{1}{k \log(k)}\right)^d.
|
||||
|
||||
|
|
|
|||
|
|
@ -30,52 +30,8 @@ where
|
|||
|
||||
T(k) = \frac{z^k \log^d(k)}{k^s}.
|
||||
|
||||
Let `U(k)` be a nonincreasing bound for the magnitude of the term ratio
|
||||
|
||||
.. math ::
|
||||
|
||||
\frac{T(k+1)}{T(k)} = z \left(\frac{k}{k+1}\right)^s
|
||||
\left( \frac{\log(k+1)}{\log(k)} \right)^d.
|
||||
|
||||
Then the remainder after the `k = N-1` term is bounded by
|
||||
|
||||
.. math ::
|
||||
|
||||
\left| \sum_{k=N}^{\infty} T(k) \right|
|
||||
\le |T(N)| \sum_{k=0}^{\infty} U(N)^k = \frac{|T(N)|}{1 - U(N)}
|
||||
|
||||
whenever `U(N) < 1`.
|
||||
|
||||
If `s = \sigma + i \tau` where `\sigma, \tau \in \mathbb{R}`, we can take
|
||||
|
||||
.. math ::
|
||||
|
||||
U(k) = |z| \; B(k, \max(0, -\sigma)) \; B(k \log(k), d)
|
||||
|
||||
wherein `B(m,n) = (1 + 1/m)^n`. This follows from the bounds
|
||||
|
||||
.. math ::
|
||||
|
||||
\left| \left(\frac{k}{k+1}\right)^s \right|
|
||||
= \left(\frac{k}{k+1}\right)^{\sigma}
|
||||
\le \begin{cases}
|
||||
1 & \text{if } \sigma \ge 0 \\
|
||||
(1 + 1/k)^{-\sigma} & \text{if } \sigma < 0.
|
||||
\end{cases}
|
||||
|
||||
and
|
||||
|
||||
.. math ::
|
||||
|
||||
\left( \frac{\log(k+1)}{\log(k)} \right)^d \le
|
||||
\left(1 + \frac{1}{k \log(k)}\right)^d.
|
||||
|
||||
We can replace `\sigma` with any lower bound, conveniently
|
||||
the nearest integer in the direction of `-\infty`, and `|z|` with any
|
||||
upper bound.
|
||||
|
||||
To bound `B(m,n)` when `m` is large, it is useful to note that
|
||||
`B(m,n) = \exp(n \log(1+1/m)) \le \exp(n/m)`).
|
||||
The remainder term `\left| \sum_{k=N}^{\infty} T(k) \right|` is bounded
|
||||
via :func:`mag_polylog_tail`.
|
||||
|
||||
Expansion for general z
|
||||
-------------------------------------------------------------------------------
|
||||
|
|
|
|||
2
mag.h
2
mag.h
|
|
@ -585,6 +585,8 @@ void mag_hypot(mag_t z, const mag_t x, const mag_t y);
|
|||
|
||||
void mag_binpow_uiui(mag_t b, ulong m, ulong n);
|
||||
|
||||
void mag_polylog_tail(mag_t u, const mag_t z, long sigma, ulong d, ulong N);
|
||||
|
||||
void mag_set_ui(mag_t z, ulong x);
|
||||
void mag_set_ui_lower(mag_t z, ulong x);
|
||||
|
||||
|
|
|
|||
113
mag/polylog_tail.c
Normal file
113
mag/polylog_tail.c
Normal file
|
|
@ -0,0 +1,113 @@
|
|||
/*=============================================================================
|
||||
|
||||
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 "mag.h"
|
||||
|
||||
void
|
||||
mag_polylog_tail(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);
|
||||
}
|
||||
|
||||
130
mag/test/t-polylog_tail.c
Normal file
130
mag/test/t-polylog_tail.c
Normal file
|
|
@ -0,0 +1,130 @@
|
|||
/*=============================================================================
|
||||
|
||||
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 "mag.h"
|
||||
#include "arb.h"
|
||||
|
||||
int main()
|
||||
{
|
||||
long iter;
|
||||
flint_rand_t state;
|
||||
|
||||
printf("polylog_tail....");
|
||||
fflush(stdout);
|
||||
|
||||
flint_randinit(state);
|
||||
|
||||
for (iter = 0; iter < 2000; iter++)
|
||||
{
|
||||
mag_t xb, yb;
|
||||
ulong N, k, d;
|
||||
long s, prec;
|
||||
arb_t z, t, u, r;
|
||||
|
||||
mag_init(xb);
|
||||
mag_init(yb);
|
||||
arb_init(z);
|
||||
arb_init(t);
|
||||
arb_init(u);
|
||||
arb_init(r);
|
||||
|
||||
mag_randtest_special(xb, state, 6);
|
||||
mag_randtest_special(yb, state, 6);
|
||||
N = n_randint(state, 100);
|
||||
d = n_randint(state, 100);
|
||||
s = n_randint(state, 100) - 50;
|
||||
prec = 4 * MAG_BITS;
|
||||
|
||||
mag_polylog_tail(yb, xb, s, d, N);
|
||||
|
||||
arb_zero(z);
|
||||
arf_set_mag(arb_midref(z), xb);
|
||||
arb_zero(r);
|
||||
|
||||
for (k = N; k < N + 100; k++)
|
||||
{
|
||||
arb_pow_ui(t, z, k, prec);
|
||||
arb_log_ui(u, k, prec);
|
||||
arb_pow_ui(u, u, d, prec);
|
||||
arb_mul(t, t, u, prec);
|
||||
|
||||
arb_set_ui(u, k);
|
||||
if (s >= 0)
|
||||
{
|
||||
arb_pow_ui(u, u, s, prec);
|
||||
arb_div(t, t, u, prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
arb_pow_ui(u, u, -s, prec);
|
||||
arb_mul(t, t, u, prec);
|
||||
}
|
||||
|
||||
arb_add(r, r, t, prec);
|
||||
}
|
||||
|
||||
MAG_CHECK_BITS(xb)
|
||||
MAG_CHECK_BITS(yb)
|
||||
|
||||
arb_zero(z);
|
||||
mag_set(arb_radref(z), yb);
|
||||
if (!arb_is_finite(z))
|
||||
arb_indeterminate(z);
|
||||
|
||||
if (!arb_contains(z, r))
|
||||
{
|
||||
printf("FAIL\n\n");
|
||||
printf("N = %lu\n\n", N);
|
||||
printf("d = %lu\n\n", d);
|
||||
printf("s = %ld\n\n", s);
|
||||
printf("xb = "); mag_printd(xb, 15); printf("\n\n");
|
||||
printf("yb = "); mag_printd(yb, 15); printf("\n\n");
|
||||
printf("z = "); arb_printd(z, 15); printf("\n\n");
|
||||
printf("r = "); arb_printd(r, 15); printf("\n\n");
|
||||
abort();
|
||||
}
|
||||
|
||||
mag_polylog_tail(xb, xb, s, d, N);
|
||||
|
||||
if (!mag_equal(xb, yb))
|
||||
{
|
||||
printf("FAIL (aliasing)\n\n");
|
||||
abort();
|
||||
}
|
||||
|
||||
mag_clear(xb);
|
||||
mag_clear(yb);
|
||||
arb_clear(z);
|
||||
arb_clear(t);
|
||||
arb_clear(u);
|
||||
arb_clear(r);
|
||||
}
|
||||
|
||||
flint_randclear(state);
|
||||
flint_cleanup();
|
||||
printf("PASS\n");
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
|
||||
Loading…
Add table
Reference in a new issue