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

some more general error bounding code for the Stirling series (work in progress)

Browse source
This commit is contained in:
Fredrik Johansson 2013-07-19 14:38:38 +02:00
parent 558e0bc8a8
commit 6231061911
6 changed files with 291 additions and 18 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
doc/source/credits.rst
View file

@ -68,6 +68,8 @@ Bibliography
.. [NIST2012] National Institute of Standards and Technology, *Digital Library of Mathematical Functions* (2012), http://dlmf.nist.gov/
.. [Olv1997] \F. Olver, *Asymptotics and special functions*, AKP Classics, AK Peters Ltd., Wellesley, MA, 1997. Reprint of the 1974 original.
.. [PS1973] \M. S. Paterson and L. J. Stockmeyer, "On the number of nonscalar multiplications necessary to evaluate polynomials", SIAM J. Comput (1973)
.. [Smi2001] \D. M. Smith, "Algorithm: Fortran 90 Software for Floating-Point Multiple Precision Arithmetic, Gamma and Related Functions", Transactions on Mathematical Software 27 (2001) 377-387, http://myweb.lmu.edu/dmsmith/toms2001.pdf

44
doc/source/gamma.rst
View file

@ -113,6 +113,50 @@ Evaluation using the Stirling series
differentiated term and the exponent `2n` changes to `2n+1`
(section 5.11 in [NIST2012]_).
.. function :: void gamma_stirling_bound_phase(fmpr_t bound, const fmpcb_t z, long prec)
Sets *bound* to an upper bound for the phase factor
`b = 1/\cos(\operatorname{arg}(z)/2)` which appears in the error bound
for the Stirling series. By trigonometric identities, assuming
that `z = x+yi`, we have `b = \sqrt{1 + t^2}` where
.. math ::
t = \frac{y}{\sqrt{x^2 + y^2} + x} = \frac{\sqrt{x^2 + y^2} - x}{y}
We bound `x` and `y` such that `|\operatorname{arg}(x+yi)|`
is maximized, and then evaluate `t` with the choice of square root
expression that avoids cancellation, using directional rounding throughout.
.. function :: void gamma_stirling_bound_fmprb(fmpr_struct * err, const fmprb_t z, long k0, long knum, long n)
.. function :: void gamma_stirling_bound_fmpcb(fmpr_struct * err, const fmpcb_t z, long k0, long knum, long n)
Computes bounds for the truncation error in the Stirling series
when summed up to term `n - 1` inclusive. An exact expression for the
truncation error is given (see [Olv1997]_ pp. 293-295) by
.. math ::
R_n(z) = \int_0^{\infty} \frac{B_{2n} - {\tilde B}_{2n}(x)}{2n(x+z)^{2n}} dt.
We optionally evaluate the bound for several terms in the
Taylor series: considering `R_n(z+t) \in \mathbb{C}[[t]]`, we
compute bounds for the coefficient of `t^k` for *knum* consecutive
values of *k* starting with *k0*.
Using the fact that the numerator of the integrand is bounded in
absolute value by `2 |B_{2n}|`, and using the bound for `|x+z|`
given by [Olv1997]_, we obtain
.. math ::
|[t^k] R_n(z+t)| \le 2 |B_{2n}|
\frac{\Gamma(2n+k-1)}{\Gamma(k+1) \Gamma(2n+1)}
\; |z| \; (b / |z|)^{2n+k}
where `b` is the phase factor implemented by
:func:`gamma_stirling_bound_phase`.
Rising factorials
--------------------------------------------------------------------------------

58
fmpcb.h
View file

@ -241,35 +241,57 @@ fmpcb_set_round_fmprb(fmpcb_t z, const fmprb_t x, long prec)
static __inline__ void
fmpcb_get_abs_ubound_fmpr(fmpr_t u, const fmpcb_t z, long prec)
{
fmpr_t v;
fmpr_init(v);
if (fmprb_is_zero(fmpcb_imagref(z)))
{
fmprb_get_abs_ubound_fmpr(u, fmpcb_realref(z), prec);
}
else if (fmprb_is_zero(fmpcb_realref(z)))
{
fmprb_get_abs_ubound_fmpr(u, fmpcb_imagref(z), prec);
}
else
{
fmpr_t v;
fmpr_init(v);
fmprb_get_abs_ubound_fmpr(u, fmpcb_realref(z), prec);
fmprb_get_abs_ubound_fmpr(v, fmpcb_imagref(z), prec);
fmprb_get_abs_ubound_fmpr(u, fmpcb_realref(z), prec);
fmprb_get_abs_ubound_fmpr(v, fmpcb_imagref(z), prec);
fmpr_mul(u, u, u, prec, FMPR_RND_UP);
fmpr_mul(v, v, v, prec, FMPR_RND_UP);
fmpr_add(u, u, v, prec, FMPR_RND_UP);
fmpr_sqrt(u, u, prec, FMPR_RND_UP);
fmpr_mul(u, u, u, prec, FMPR_RND_UP);
fmpr_mul(v, v, v, prec, FMPR_RND_UP);
fmpr_add(u, u, v, prec, FMPR_RND_UP);
fmpr_sqrt(u, u, prec, FMPR_RND_UP);
fmpr_clear(v);
fmpr_clear(v);
}
}
static __inline__ void
fmpcb_get_abs_lbound_fmpr(fmpr_t u, const fmpcb_t z, long prec)
{
fmpr_t v;
fmpr_init(v);
if (fmprb_is_zero(fmpcb_imagref(z)))
{
fmprb_get_abs_lbound_fmpr(u, fmpcb_realref(z), prec);
}
else if (fmprb_is_zero(fmpcb_realref(z)))
{
fmprb_get_abs_lbound_fmpr(u, fmpcb_imagref(z), prec);
}
else
{
fmpr_t v;
fmpr_init(v);
fmprb_get_abs_lbound_fmpr(u, fmpcb_realref(z), prec);
fmprb_get_abs_lbound_fmpr(v, fmpcb_imagref(z), prec);
fmprb_get_abs_lbound_fmpr(u, fmpcb_realref(z), prec);
fmprb_get_abs_lbound_fmpr(v, fmpcb_imagref(z), prec);
fmpr_mul(u, u, u, prec, FMPR_RND_DOWN);
fmpr_mul(v, v, v, prec, FMPR_RND_DOWN);
fmpr_add(u, u, v, prec, FMPR_RND_DOWN);
fmpr_sqrt(u, u, prec, FMPR_RND_DOWN);
fmpr_mul(u, u, u, prec, FMPR_RND_DOWN);
fmpr_mul(v, v, v, prec, FMPR_RND_DOWN);
fmpr_add(u, u, v, prec, FMPR_RND_DOWN);
fmpr_sqrt(u, u, prec, FMPR_RND_DOWN);
fmpr_clear(v);
fmpr_clear(v);
}
}
static __inline__ void

4
gamma.h
View file

@ -71,6 +71,10 @@ void gamma_taylor_fmprb(fmprb_t y, const fmprb_t x, long prec);
void gamma_stirling_choose_param_fmprb(int * reflect, long * r, long * n, const fmprb_t z, int use_reflect, int digamma, long prec);
void gamma_stirling_choose_param_fmpcb(int * reflect, long * r, long * n, const fmpcb_t z, int use_reflect, int digamma, long prec);
void gamma_stirling_bound_phase(fmpr_t bound, const fmpcb_t z, long prec);
void gamma_stirling_bound_fmprb(fmpr_struct * err, const fmprb_t x, long k0, long knum, long n);
void gamma_stirling_bound_fmpcb(fmpr_struct * err, const fmpcb_t z, long k0, long knum, long n);
void gamma_stirling_coeff(fmprb_t b, ulong k, int digamma, long prec);
void gamma_stirling_eval_series_fmprb(fmprb_t s, const fmprb_t z, long nterms, int digamma, long prec);
void gamma_stirling_eval_series_fmpcb(fmpcb_t s, const fmpcb_t z, long nterms, int digamma, long prec);

108
gamma/stirling_bound.c Normal file
View file

@ -0,0 +1,108 @@
/*=============================================================================
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) 2012 Fredrik Johansson
******************************************************************************/
#include <math.h>
#include "gamma.h"
#include "bernoulli.h"
void fmpr_gamma_ui_lbound(fmpr_t x, ulong n, long prec);
void fmpr_gamma_ui_ubound(fmpr_t x, ulong n, long prec);
/*
2 |B_{2n}| G(2n+k-1) / (G(k+1) G(2n+1)) |z| (T |z|^{-1})^(2n+k)
*/
void
gamma_stirling_bound_fmpcb(fmpr_struct * err, const fmpcb_t z, long k0, long knum, long n)
{
fmpr_t c, t, u, v;
long i, k, prec = FMPRB_RAD_PREC;
if (fmprb_contains_zero(fmpcb_imagref(z)) &&
fmprb_contains_nonpositive(fmpcb_realref(z)))
{
for (i = 0; i < knum; i++)
fmpr_pos_inf(err + i);
return;
}
fmpr_init(c);
fmpr_init(t);
fmpr_init(u);
fmpr_init(v);
/* t = lower bound for |z| */
fmpcb_get_abs_lbound_fmpr(t, z, prec);
fmpcb_get_abs_ubound_fmpr(v, z, prec);
/* c = upper bound for 1/(cos(arg(z)/2) |z|) */
gamma_stirling_bound_phase(c, z, prec);
fmpr_div(c, c, t, prec, FMPR_RND_UP);
/* numerator: 2 B_{2n} gamma(2n+k-1) |z| */
BERNOULLI_ENSURE_CACHED(2 * n);
fmpr_set_round_fmpz(err, fmpq_numref(bernoulli_cache + 2 * n), prec, FMPR_RND_UP);
fmpr_abs(err, err);
fmpr_div_fmpz(err, err, fmpq_denref(bernoulli_cache + 2 * n), prec, FMPR_RND_UP);
fmpr_mul_2exp_si(err, err, 1);
fmpr_gamma_ui_ubound(u, 2 * n + k0 - 1, prec);
fmpr_mul(err, err, u, prec, FMPR_RND_UP);
fmpr_mul(err, err, v, prec, FMPR_RND_UP);
/* denominator gamma(k+1) gamma(2n+1) */
fmpr_gamma_ui_lbound(t, 2 * n + 1, prec);
fmpr_gamma_ui_lbound(u, k0 + 1, prec);
fmpr_mul(t, t, u, prec, FMPR_RND_DOWN);
fmpr_div(err, err, t, prec, FMPR_RND_UP);
/* multiply by c^(2n+k) */
fmpr_pow_sloppy_ui(t, c, 2 * n + k0, prec, FMPR_RND_UP);
fmpr_mul(err, err, t, prec, FMPR_RND_UP);
for (i = 1; i < knum; i++)
{
/* recurrence factor: c * (2n+k-2) / k */
k = k0 + i;
fmpr_mul(err + i, err + i - 1, c, prec, FMPR_RND_UP);
fmpr_mul_ui(err + i, err + i, 2 * n + k - 2, prec, FMPR_RND_UP);
fmpr_div_ui(err + i, err + i, k, prec, FMPR_RND_UP);
}
fmpr_clear(c);
fmpr_clear(t);
fmpr_clear(u);
fmpr_clear(v);
}
void
gamma_stirling_bound_fmprb(fmpr_struct * err, const fmprb_t x, long k0, long knum, long n)
{
fmpcb_t z;
fmpcb_init(z);
fmpcb_set_fmprb(z, x);
gamma_stirling_bound_fmpcb(err, z, k0, knum, n);
fmpcb_clear(z);
}

93
gamma/stirling_bound_phase.c Normal file
View file

@ -0,0 +1,93 @@
/*=============================================================================
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 "gamma.h"
void
gamma_stirling_bound_phase(fmpr_t bound, const fmpcb_t z, long prec)
{
fmpr_t x, y, t, u;
int xsign;
fmpr_init(x);
fmpr_init(y);
fmpr_init(t);
fmpr_init(u);
/* first compute x, y such that |arg(z)| <= arg(x+yi) */
/* argument increases with smaller real parts */
fmpr_sub(x, fmprb_midref(fmpcb_realref(z)),
fmprb_radref(fmpcb_realref(z)), prec, FMPR_RND_FLOOR);
xsign = fmpr_sgn(x);
if (xsign >= 0) /* argument increases away from the real axis */
fmprb_get_abs_ubound_fmpr(y, fmpcb_imagref(z), prec);
else /* argument increases closer to the real axis */
fmprb_get_abs_lbound_fmpr(y, fmpcb_imagref(z), prec);
if (fmpr_is_zero(y))
{
if (xsign > 0)
fmpr_one(bound);
else
fmpr_pos_inf(bound);
}
else
{
if (xsign >= 0)
{
/* compute upper bound for t = y / (sqrt(x^2 + y^2) + x) */
fmpr_mul(t, x, x, prec, FMPR_RND_DOWN);
fmpr_mul(u, y, y, prec, FMPR_RND_DOWN);
fmpr_add(t, t, u, prec, FMPR_RND_DOWN);
fmpr_sqrt(t, t, prec, FMPR_RND_DOWN);
fmpr_add(t, t, x, prec, FMPR_RND_DOWN);
fmpr_div(t, y, t, prec, FMPR_RND_UP);
}
else
{
/* compute upper bound for t = (sqrt(x^2 + y^2) - x) / y */
fmpr_mul(t, x, x, prec, FMPR_RND_UP);
fmpr_mul(u, y, y, prec, FMPR_RND_UP);
fmpr_add(t, t, u, prec, FMPR_RND_UP);
fmpr_sqrt(t, t, prec, FMPR_RND_UP);
fmpr_sub(t, t, x, prec, FMPR_RND_UP);
fmpr_div(t, t, y, prec, FMPR_RND_UP);
}
/* compute upper bound for sqrt(1 + t^2) */
fmpr_mul(t, t, t, prec, FMPR_RND_UP);
fmpr_add_ui(t, t, 1, prec, FMPR_RND_UP);
fmpr_sqrt(bound, t, prec, FMPR_RND_UP);
}
fmpr_clear(x);
fmpr_clear(y);
fmpr_clear(t);
fmpr_clear(u);
}