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 renaming

Browse source
This commit is contained in:
Fredrik Johansson 2013-07-30 15:21:17 +02:00
parent 6f62daaccf
commit d6ef6f70b2
11 changed files with 28 additions and 135 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

18
doc/source/gamma.rst
View file

@ -82,9 +82,9 @@ Evaluation using the Stirling series
Sets `b = B_{2k} / (2k (2k-1))`, rounded to *prec* bits, or if *digamma*
is nonzero, sets `b = B_{2k} / (2k)`.
.. function :: void gamma_stirling_eval_series_fmprb(fmprb_t s, const fmprb_t z, long n, int digamma, long prec)
.. function :: void gamma_stirling_eval_fmprb(fmprb_t s, const fmprb_t z, long n, int digamma, long prec)
.. function :: void gamma_stirling_eval_series_fmpcb(fmpcb_t s, const fmpcb_t z, long n, int digamma, long prec)
.. function :: void gamma_stirling_eval_fmpcb(fmpcb_t s, const fmpcb_t z, long n, int digamma, long prec)
Evaluates the Stirling series
@ -101,17 +101,7 @@ Evaluation using the Stirling series
If *digamma* is nonzero, the derivative of this series (i.e. the
expansion for the digamma function) is evaluated.
The error bound
.. math ::
|R(n,z)| \le \frac{|t_n|}{\cos(0.5 \arg(z))^{2n}}
is included in the output (when evaluating the digamma function, the
expression is the same except that `t_n` changes to the
differentiated term and the exponent `2n` changes to `2n+1`
(section 5.11 in [NIST2012]_).
The error bound for the tail `R(n,z)` is included in the output.
.. function :: void gamma_stirling_bound_phase(fmpr_t bound, const fmpcb_t z, long prec)
@ -138,7 +128,7 @@ Evaluation using the Stirling series
.. math ::
R_n(z) = \int_0^{\infty} \frac{B_{2n} - {\tilde B}_{2n}(x)}{2n(x+z)^{2n}} dt.
R_n(z) = \int_0^{\infty} \frac{B_{2n} - {\tilde B}_{2n}(x)}{2n(x+z)^{2n}} dx.
We optionally evaluate the bound for several terms in the
Taylor series: considering `R_n(z+t) \in \mathbb{C}[[t]]`, we

4
fmpcb/digamma.c
View file

@ -52,13 +52,13 @@ fmpcb_digamma(fmpcb_t y, const fmpcb_t x, long prec)
gamma_harmonic_sum_fmpcb_ui_bsplit(u, t, r, wp);
fmpcb_add(v, v, u, wp);
fmpcb_add_ui(t, t, r, wp);
gamma_stirling_eval_series_fmpcb(u, t, n, 1, wp);
gamma_stirling_eval_fmpcb(u, t, n, 1, wp);
fmpcb_sub(y, u, v, wp);
}
else
{
fmpcb_add_ui(t, x, r, wp);
gamma_stirling_eval_series_fmpcb(u, t, n, 1, wp);
gamma_stirling_eval_fmpcb(u, t, n, 1, wp);
gamma_harmonic_sum_fmpcb_ui_bsplit(t, x, r, wp);
fmpcb_sub(y, u, t, prec);
}

6
fmpcb/gamma.c
View file

@ -50,7 +50,7 @@ _fmpcb_gamma(fmpcb_t y, const fmpcb_t x, long prec, int inverse)
fmprb_const_pi(fmpcb_realref(v), wp);
fmpcb_mul_fmprb(u, u, fmpcb_realref(v), wp);
fmpcb_add_ui(t, t, r, wp);
gamma_stirling_eval_series_fmpcb(v, t, n, 0, wp);
gamma_stirling_eval_fmpcb(v, t, n, 0, wp);
fmpcb_exp(v, v, wp);
fmpcb_sin_pi(t, x, wp);
fmpcb_mul(v, v, t, wp);
@ -59,7 +59,7 @@ _fmpcb_gamma(fmpcb_t y, const fmpcb_t x, long prec, int inverse)
{
/* gamma(x) = gamma(x+r) / rf(x,r) */
fmpcb_add_ui(t, x, r, wp);
gamma_stirling_eval_series_fmpcb(u, t, n, 0, wp);
gamma_stirling_eval_fmpcb(u, t, n, 0, wp);
fmpcb_exp(u, u, prec);
gamma_rising_fmpcb_ui_bsplit(v, x, r, wp);
}
@ -126,7 +126,7 @@ fmpcb_lgamma(fmpcb_t y, const fmpcb_t x, long prec)
fmpcb_init(u);
fmpcb_add_ui(t, x, r, wp);
gamma_stirling_eval_series_fmpcb(u, t, n, 0, wp);
gamma_stirling_eval_fmpcb(u, t, n, 0, wp);
fmpcb_log_rfac_ui(t, x, r, wp);
fmpcb_sub(y, u, t, prec);

4
fmprb/digamma.c
View file

@ -52,13 +52,13 @@ fmprb_digamma(fmprb_t y, const fmprb_t x, long prec)
gamma_harmonic_sum_fmprb_ui_bsplit(u, t, r, wp);
fmprb_add(v, v, u, wp);
fmprb_add_ui(t, t, r, wp);
gamma_stirling_eval_series_fmprb(u, t, n, 1, wp);
gamma_stirling_eval_fmprb(u, t, n, 1, wp);
fmprb_sub(y, u, v, wp);
}
else
{
fmprb_add_ui(t, x, r, wp);
gamma_stirling_eval_series_fmprb(u, t, n, 1, wp);
gamma_stirling_eval_fmprb(u, t, n, 1, wp);
gamma_harmonic_sum_fmprb_ui_bsplit(t, x, r, wp);
fmprb_sub(y, u, t, prec);
}

6
fmprb/gamma.c
View file

@ -92,7 +92,7 @@ _fmprb_gamma(fmprb_t y, const fmprb_t x, long prec, int inverse)
fmprb_const_pi(v, wp);
fmprb_mul(u, u, v, wp);
fmprb_add_ui(t, t, r, wp);
gamma_stirling_eval_series_fmprb(v, t, n, 0, wp);
gamma_stirling_eval_fmprb(v, t, n, 0, wp);
fmprb_exp(v, v, wp);
fmprb_sin_pi(t, x, wp);
fmprb_mul(v, v, t, wp);
@ -101,7 +101,7 @@ _fmprb_gamma(fmprb_t y, const fmprb_t x, long prec, int inverse)
{
/* gamma(x) = gamma(x+r) / rf(x,r) */
fmprb_add_ui(t, x, r, wp);
gamma_stirling_eval_series_fmprb(u, t, n, 0, wp);
gamma_stirling_eval_fmprb(u, t, n, 0, wp);
fmprb_exp(u, u, prec);
gamma_rising_fmprb_ui_bsplit(v, x, r, wp);
}
@ -144,7 +144,7 @@ fmprb_lgamma(fmprb_t y, const fmprb_t x, long prec)
fmprb_init(u);
fmprb_add_ui(t, x, r, wp);
gamma_stirling_eval_series_fmprb(u, t, n, 0, wp);
gamma_stirling_eval_fmprb(u, t, n, 0, wp);
gamma_rising_fmprb_ui_bsplit(t, x, r, wp);
fmprb_log(t, t, wp);
fmprb_sub(y, u, t, prec);

7
gamma.h
View file

@ -64,7 +64,7 @@ gamma_taylor_coeffs_for_prec(long prec)
}
void gamma_taylor_precompute(long num, long prec);
void gamma_taylor_eval_series_fmprb(fmprb_t y, const fmprb_t x, long prec);
void gamma_taylor_eval_fmprb(fmprb_t y, const fmprb_t x, long prec);
void gamma_taylor_fmprb(fmprb_t y, const fmprb_t x, long prec);
@ -76,8 +76,9 @@ void gamma_stirling_bound_fmprb(fmpr_struct * err, const fmprb_t x, long k0, lon
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);
void gamma_stirling_eval_fmprb(fmprb_t s, const fmprb_t z, long nterms, int digamma, long prec);
void gamma_stirling_eval_fmpcb(fmpcb_t s, const fmpcb_t z, long nterms, int digamma, long prec);
void gamma_stirling_eval_fmprb_series(fmprb_ptr res, const fmprb_t z, long n, long num, long prec);
void gamma_stirling_eval_fmpcb_series(fmpcb_ptr res, const fmpcb_t z, long n, long num, long prec);

4
gamma/fmpq_stirling.c
View file

@ -56,7 +56,7 @@ gamma_fmpq_stirling(fmprb_t y, const fmpq_t a, long prec)
fmprb_sub_ui(t, x, 1, wp);
fmprb_neg(t, t);
fmprb_add_ui(t, t, r, wp);
gamma_stirling_eval_series_fmprb(v, t, n, 0, wp);
gamma_stirling_eval_fmprb(v, t, n, 0, wp);
fmprb_exp(v, v, wp);
fmprb_sin_pi_fmpq(t, a, wp);
fmprb_mul(v, v, t, wp);
@ -65,7 +65,7 @@ gamma_fmpq_stirling(fmprb_t y, const fmpq_t a, long prec)
{
/* gamma(x) = gamma(x+r) / rf(x,r) */
fmprb_add_ui(t, x, r, wp);
gamma_stirling_eval_series_fmprb(u, t, n, 0, wp);
gamma_stirling_eval_fmprb(u, t, n, 0, wp);
fmprb_exp(u, u, prec);
gamma_rising_fmprb_fmpq_ui_bsplit(v, a, r, wp);
}

61
gamma/stirling_eval_series_fmpcb.c → gamma/stirling_eval_fmpcb.c
View file

@ -27,65 +27,8 @@
#include "gamma.h"
#include "bernoulli.h"
/*
B_(2n) / (2n (2n-1)) / |z|^(2n-1) * (1/cos(0.5*arg(z))^(2n))
*/
void
gamma_stirling_bound_remainder_fmpcb(fmpr_t err, const fmpcb_t z, int digamma, long n)
{
fmpr_t t;
fmprb_t b;
if (fmprb_contains_zero(fmpcb_imagref(z)) &&
fmprb_contains_nonpositive(fmpcb_realref(z)))
{
fmpr_pos_inf(err);
return;
}
/* bound 1 / |z|^(2n-1) */
fmpcb_get_abs_lbound_fmpr(err, z, FMPRB_RAD_PREC);
if (fmpr_is_zero(err))
{
fmpr_pos_inf(err);
return;
}
fmpr_ui_div(err, 1, err, FMPRB_RAD_PREC, FMPR_RND_UP);
if (digamma)
fmpr_pow_sloppy_ui(err, err, 2 * n, FMPRB_RAD_PREC, FMPR_RND_UP);
else
fmpr_pow_sloppy_ui(err, err, 2 * n - 1, FMPRB_RAD_PREC, FMPR_RND_UP);
/* bound coefficient */
fmprb_init(b);
fmpr_init(t);
gamma_stirling_coeff(b, n, digamma, FMPRB_RAD_PREC);
fmprb_get_abs_ubound_fmpr(t, b, FMPRB_RAD_PREC);
fmpr_mul(err, err, t, FMPRB_RAD_PREC, FMPR_RND_UP);
/* bound 1/cos(0.5*arg(z))^(2n) */
fmpcb_arg(b, z, FMPRB_RAD_PREC);
fmprb_mul_2exp_si(b, b, -1);
fmprb_cos(b, b, FMPRB_RAD_PREC);
fmprb_get_abs_lbound_fmpr(t, b, FMPRB_RAD_PREC);
fmpr_ui_div(t, 1, t, FMPRB_RAD_PREC, FMPR_RND_UP);
if (digamma)
fmpr_pow_sloppy_ui(t, t, 2 * n + 1, FMPRB_RAD_PREC, FMPR_RND_UP);
else
fmpr_pow_sloppy_ui(t, t, 2 * n, FMPRB_RAD_PREC, FMPR_RND_UP);
fmpr_mul(err, err, t, FMPRB_RAD_PREC, FMPR_RND_UP);
fmprb_clear(b);
fmpr_clear(t);
}
void
gamma_stirling_eval_series_fmpcb(fmpcb_t s, const fmpcb_t z, long nterms, int digamma, long prec)
gamma_stirling_eval_fmpcb(fmpcb_t s, const fmpcb_t z, long nterms, int digamma, long prec)
{
fmpcb_t t, logz, zinv, zinv2;
fmprb_t b;
@ -141,7 +84,7 @@ gamma_stirling_eval_series_fmpcb(fmpcb_t s, const fmpcb_t z, long nterms, int di
/* remainder bound */
fmpr_init(err);
gamma_stirling_bound_remainder_fmpcb(err, z, digamma, nterms);
gamma_stirling_bound_fmpcb(err, z, digamma ? 1 : 0, 1, nterms);
fmprb_add_error_fmpr(fmpcb_realref(s), err);
fmprb_add_error_fmpr(fmpcb_imagref(s), err);
fmpr_clear(err);

45
gamma/stirling_eval_series_fmprb.c → gamma/stirling_eval_fmprb.c
View file

@ -28,48 +28,7 @@
#include "bernoulli.h"
void
gamma_stirling_bound_remainder(fmpr_t err, const fmprb_t z, int digamma, long n)
{
if (fmprb_contains_nonpositive(z))
{
fmpr_pos_inf(err);
}
else
{
fmpr_t t;
fmprb_t b;
fmpr_init(t);
fmprb_init(b);
fmpr_sub(t, fmprb_midref(z), fmprb_radref(z), FMPRB_RAD_PREC, FMPR_RND_FLOOR);
if (fmpr_sgn(t) <= 0)
{
fmpr_pos_inf(err);
}
else
{
fmpr_ui_div(t, 1, t, FMPRB_RAD_PREC, FMPR_RND_UP);
if (digamma)
fmpr_pow_sloppy_ui(t, t, 2 * n, FMPRB_RAD_PREC, FMPR_RND_UP);
else
fmpr_pow_sloppy_ui(t, t, 2 * n - 1, FMPRB_RAD_PREC, FMPR_RND_UP);
gamma_stirling_coeff(b, n, digamma, FMPRB_RAD_PREC);
fmprb_get_abs_ubound_fmpr(err, b, FMPRB_RAD_PREC);
fmpr_mul(err, err, t, FMPRB_RAD_PREC, FMPR_RND_UP);
}
fmpr_clear(t);
fmprb_clear(b);
}
}
void
gamma_stirling_eval_series_fmprb(fmprb_t s, const fmprb_t z, long nterms, int digamma, long prec)
gamma_stirling_eval_fmprb(fmprb_t s, const fmprb_t z, long nterms, int digamma, long prec)
{
fmprb_t b, t, logz, zinv, zinv2;
fmpr_t err;
@ -124,7 +83,7 @@ gamma_stirling_eval_series_fmprb(fmprb_t s, const fmprb_t z, long nterms, int di
/* remainder bound */
fmpr_init(err);
gamma_stirling_bound_remainder(err, z, digamma, nterms);
gamma_stirling_bound_fmprb(err, z, digamma ? 1 : 0, 1, nterms);
fmprb_add_error_fmpr(s, err);
fmpr_clear(err);

2
gamma/taylor_eval_series_fmprb.c
View file

@ -28,7 +28,7 @@
/* evaluate sum_{k=0}^{n-1} c_{k+1} x^k
precision estimate assumes x <= 1/2 */
void
gamma_taylor_eval_series_fmprb(fmprb_t y, const fmprb_t x, long prec)
gamma_taylor_eval_fmprb(fmprb_t y, const fmprb_t x, long prec)
{
long i, n, wp;
fmprb_t t, u, v;

6
gamma/taylor_fmprb.c
View file

@ -40,14 +40,14 @@ gamma_taylor_fmprb(fmprb_t y, const fmprb_t x, long prec)
if (v == 0)
{
gamma_taylor_eval_series_fmprb(u, x, prec);
gamma_taylor_eval_fmprb(u, x, prec);
fmprb_mul(u, u, x, prec);
fmprb_ui_div(y, 1, u, prec);
}
else if (v > 0)
{
fmprb_sub_si(t, x, v, prec);
gamma_taylor_eval_series_fmprb(t, t, prec);
gamma_taylor_eval_fmprb(t, t, prec);
fmprb_sub_si(u, x, v - 1, prec);
gamma_rising_fmprb_ui_bsplit(u, u, v - 1, prec);
fmprb_div(y, u, t, prec);
@ -55,7 +55,7 @@ gamma_taylor_fmprb(fmprb_t y, const fmprb_t x, long prec)
else
{
fmprb_add_si(t, x, (-v), prec);
gamma_taylor_eval_series_fmprb(t, t, prec);
gamma_taylor_eval_fmprb(t, t, prec);
gamma_rising_fmprb_ui_bsplit(u, x, (-v) + 1, prec);
fmprb_mul(y, u, t, prec);
fmprb_ui_div(y, 1, y, prec);