diff --git a/doc/source/gamma.rst b/doc/source/gamma.rst index 1c6dd3ac..f4cab2d3 100644 --- a/doc/source/gamma.rst +++ b/doc/source/gamma.rst @@ -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 diff --git a/fmpcb/digamma.c b/fmpcb/digamma.c index 420bfc01..27fbef55 100644 --- a/fmpcb/digamma.c +++ b/fmpcb/digamma.c @@ -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); } diff --git a/fmpcb/gamma.c b/fmpcb/gamma.c index a3965458..55f110db 100644 --- a/fmpcb/gamma.c +++ b/fmpcb/gamma.c @@ -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); diff --git a/fmprb/digamma.c b/fmprb/digamma.c index 5b611373..4e1747ac 100644 --- a/fmprb/digamma.c +++ b/fmprb/digamma.c @@ -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); } diff --git a/fmprb/gamma.c b/fmprb/gamma.c index 27527a1a..52d11f07 100644 --- a/fmprb/gamma.c +++ b/fmprb/gamma.c @@ -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); diff --git a/gamma.h b/gamma.h index 5306301b..dd08dd71 100644 --- a/gamma.h +++ b/gamma.h @@ -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); diff --git a/gamma/fmpq_stirling.c b/gamma/fmpq_stirling.c index 00883376..62298167 100644 --- a/gamma/fmpq_stirling.c +++ b/gamma/fmpq_stirling.c @@ -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); } diff --git a/gamma/stirling_eval_series_fmpcb.c b/gamma/stirling_eval_fmpcb.c similarity index 66% rename from gamma/stirling_eval_series_fmpcb.c rename to gamma/stirling_eval_fmpcb.c index 521a8505..13adb4db 100644 --- a/gamma/stirling_eval_series_fmpcb.c +++ b/gamma/stirling_eval_fmpcb.c @@ -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); diff --git a/gamma/stirling_eval_series_fmprb.c b/gamma/stirling_eval_fmprb.c similarity index 73% rename from gamma/stirling_eval_series_fmprb.c rename to gamma/stirling_eval_fmprb.c index ae1a5f6c..d00e35e4 100644 --- a/gamma/stirling_eval_series_fmprb.c +++ b/gamma/stirling_eval_fmprb.c @@ -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); diff --git a/gamma/taylor_eval_series_fmprb.c b/gamma/taylor_eval_series_fmprb.c index 9ad7afb2..febe419e 100644 --- a/gamma/taylor_eval_series_fmprb.c +++ b/gamma/taylor_eval_series_fmprb.c @@ -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; diff --git a/gamma/taylor_fmprb.c b/gamma/taylor_fmprb.c index 088c111b..488314c3 100644 --- a/gamma/taylor_fmprb.c +++ b/gamma/taylor_fmprb.c @@ -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);