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

use reflection formula in lgamma_series

Browse source
This commit is contained in:
Fredrik Johansson 2015-07-06 15:52:54 +02:00
parent 91c72f51c5
commit 250b2921dc
4 changed files with 88 additions and 10 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
acb_poly/gamma_series.c
View file

@ -45,7 +45,7 @@ bsplit(acb_ptr Q, acb_ptr T, const acb_t z, long a, long b, long num, long prec)
acb_set(Q, z);
if (num > 1) acb_one(Q + 1);
if (num > 2) acb_zero(Q + 2);
}
}
else
{ /* (z + t)^2 */
acb_mul(Q, z, z, prec); /* TODO: precompute */

55
acb_poly/lgamma_series.c
View file

@ -109,14 +109,57 @@ _acb_poly_lgamma_series(acb_ptr res, acb_srcptr h, long hlen, long len, long pre
acb_init(zr);
/* use Stirling series */
acb_gamma_stirling_choose_param(&reflect, &r, &n, h, 0, 0, wp);
acb_add_ui(zr, h, r, wp);
_acb_poly_gamma_stirling_eval(u, zr, n, len, wp);
acb_gamma_stirling_choose_param(&reflect, &r, &n, h, 1, 0, wp);
if (r != 0)
if (reflect)
{
_log_rising_ui_series(t, h, r, len, wp);
_acb_vec_sub(u, u, t, len, wp);
/* log gamma(h+x) = log rf(1-(h+x), r) - log gamma(1-(h+x)+r) - log sin(pi (h+x)) + log(pi) */
if (r != 0) /* otherwise t = 0 */
{
acb_sub_ui(u, h, 1, wp);
acb_neg(u, u);
_log_rising_ui_series(t, u, r, len, wp);
for (i = 1; i < len; i += 2)
acb_neg(t + i, t + i);
}
acb_sub_ui(u, h, 1, wp);
acb_neg(u, u);
acb_add_ui(zr, u, r, wp);
_acb_poly_gamma_stirling_eval(u, zr, n, len, wp);
for (i = 1; i < len; i += 2)
acb_neg(u + i, u + i);
_acb_vec_sub(t, t, u, len, wp);
/* log(sin) is unstable with large imaginary parts;
cot_pi is implemented in a numerically stable way */
acb_set(u, h);
acb_one(u + 1);
_acb_poly_cot_pi_series(u, u, 2, len - 1, wp);
_acb_poly_integral(u, u, len, wp);
acb_const_pi(u, wp);
_acb_vec_scalar_mul(u + 1, u + 1, len - 1, u, wp);
acb_log_sin_pi(u, h, wp);
_acb_vec_sub(u, t, u, len, wp);
acb_const_pi(t, wp); /* todo: constant for log pi */
acb_log(t, t, wp);
acb_add(u, u, t, wp);
}
else
{
/* log gamma(x) = log gamma(x+r) - log rf(x,r) */
acb_add_ui(zr, h, r, wp);
_acb_poly_gamma_stirling_eval(u, zr, n, len, wp);
if (r != 0)
{
_log_rising_ui_series(t, h, r, len, wp);
_acb_vec_sub(u, u, t, len, wp);
}
}
/* compose with nonconstant part */

35
acb_poly/test/t-lgamma_series.c
View file

@ -35,6 +35,35 @@ int main()
flint_randinit(state);
/* special accuracy test case */
{
acb_poly_t a;
acb_t c;
acb_init(c);
acb_poly_init(a);
arb_set_str(acb_realref(c), "-20.25", 53);
arb_set_str(acb_imagref(c), "1e1000", 53);
acb_poly_set_coeff_acb(a, 0, c);
acb_poly_set_coeff_si(a, 1, 1);
acb_poly_lgamma_series(a, a, 3, 53);
if (acb_rel_accuracy_bits(a->coeffs) < 40 ||
acb_rel_accuracy_bits(a->coeffs + 1) < 40 ||
acb_rel_accuracy_bits(a->coeffs + 2) < 40)
{
printf("FAIL: accuracy (reflection formula)\n\n");
acb_poly_printd(a, 15); printf("\n\n");
abort();
}
acb_poly_clear(a);
acb_clear(c);
}
for (iter = 0; iter < 500; iter++)
{
long m, n1, n2, rbits1, rbits2, rbits3;
@ -53,9 +82,9 @@ int main()
acb_poly_init(c);
acb_poly_init(d);
acb_poly_randtest(a, state, m, rbits1, 3);
acb_poly_randtest(b, state, m, rbits1, 3);
acb_poly_randtest(c, state, m, rbits1, 3);
acb_poly_randtest(a, state, m, rbits1, 10);
acb_poly_randtest(b, state, m, rbits1, 10);
acb_poly_randtest(c, state, m, rbits1, 10);
acb_poly_lgamma_series(b, a, n1, rbits2);
acb_poly_lgamma_series(c, a, n2, rbits3);

6
arb_poly/lgamma_series.c
View file

@ -61,6 +61,12 @@ _arb_poly_lgamma_series(arb_ptr res, arb_srcptr h, long hlen, long len, long pre
arb_t zr;
arb_ptr t, u;
if (!arb_is_positive(h))
{
_arb_vec_indeterminate(res, len);
return;
}
hlen = FLINT_MIN(hlen, len);
wp = prec + FLINT_BIT_COUNT(prec);