arb_hypgeom_lgamma, acb_hypgeom_lgamma

This commit is contained in:
fredrik 2021-09-01 09:52:11 +02:00
parent 83b9ae5b58
commit 97f2325896
8 changed files with 592 additions and 4 deletions

View file

@ -41,6 +41,8 @@ int acb_hypgeom_gamma_taylor(acb_t res, const acb_t x, int reciprocal, slong pre
void acb_hypgeom_gamma(acb_t y, const acb_t x, slong prec);
void acb_hypgeom_rgamma(acb_t y, const acb_t x, slong prec);
void acb_hypgeom_lgamma(acb_t y, const acb_t x, slong prec);
void acb_hypgeom_pfq_bound_factor(mag_t C,
acb_srcptr a, slong p, acb_srcptr b, slong q, const acb_t z, ulong n);

262
acb_hypgeom/lgamma.c Normal file
View file

@ -0,0 +1,262 @@
/*
Copyright (C) 2014, 2015, 2021 Fredrik Johansson
This file is part of Arb.
Arb is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version. See <http://www.gnu.org/licenses/>.
*/
#include "acb_hypgeom.h"
#include "arb_hypgeom.h"
void acb_hypgeom_gamma_stirling_choose_param(int * reflect, slong * r, slong * n,
const acb_t z, int use_reflect, int digamma, slong prec);
void acb_hypgeom_gamma_stirling_inner(acb_t s, const acb_t z, slong N, slong prec);
static double
want_taylor(double x, double y, slong prec)
{
if (y < 0.0) y = -y;
if (x < 0.0) x = -2.0 * x;
if ((prec < 128 && y > 4.0) || (prec < 256 && y > 5.0) ||
(prec < 512 && y > 8.0) || (prec < 1024 && y > 9.0) || y > 10.0)
{
return 0;
}
if (x * (1.0 + 0.75 * y) > 8 + 0.15 * prec)
{
return 0;
}
return 1;
}
/* Linear fit on [0.5, 1.5] for
lambda x: findroot(lambda y: im(loggamma(x+1j*y)) - (n+0.5)*pi */
static const double Atab[] = {
4.5835631239879990091,
6.4037921417161376741,
7.9938623618272375768,
9.4449131928216797873,
10.802608819487725856,
12.0918817314347272,
};
static const double Btab[] = {
-1.1432582881376479127,
-0.86248117216701645437,
-0.75778990135448922722,
-0.69734688055939976228,
-0.65626499937495627271,
-0.62578331900739100617,
};
void
_arb_const_log_pi(arb_t t, slong prec)
{
arb_const_pi(t, prec + 2);
arb_log(t, t, prec);
}
ARB_DEF_CACHED_CONSTANT(arb_const_log_pi, _arb_const_log_pi)
int
acb_hypgeom_lgamma_taylor(acb_t res, const acb_t z, slong prec)
{
double x, y, acc;
slong k, r, wp;
acb_t t, u;
int reflect;
/* Assume xerr, yerr <= 1/16 */
if (mag_cmp_2exp_si(arb_radref(acb_realref(z)), -4) > 0)
return 0;
if (mag_cmp_2exp_si(arb_radref(acb_imagref(z)), -4) > 0)
return 0;
acc = acb_rel_accuracy_bits(z);
acc = FLINT_MAX(acc, 0);
wp = FLINT_MIN(prec, acc + 20);
wp = FLINT_MAX(wp, 2);
/* x, y plus eventual rounding error */
x = arf_get_d(arb_midref(acb_realref(z)), ARF_RND_NEAR);
y = arf_get_d(arb_midref(acb_imagref(z)), ARF_RND_NEAR);
if (!want_taylor(x, y, wp))
return 0;
acb_init(t);
acb_init(u);
/* Reduce real part to (approximately) [0.5, 1.5]. */
r = floor(x - 0.5);
/* Reflection formula is slower but improves accuracy. */
reflect = (x < -3.0);
if (reflect)
{
acb_neg(u, z);
acb_add_si(u, u, 2 + r, 2 * prec + 10);
x = 2.0 + r - x;
y = -y;
}
else
{
acb_sub_si(u, z, r, 2 * prec + 10);
x = x - r;
}
for (k = 0; k < 6; k++)
{
if (fabs(y) <= Atab[k] + Btab[k] * x)
{
if (!acb_hypgeom_gamma_taylor(t, u, 1, wp))
{
acb_clear(t);
acb_clear(u);
return 0;
}
if (k % 2 == 0)
{
acb_log(t, t, wp);
acb_neg(t, t);
}
else
{
acb_neg(t, t);
acb_log(t, t, wp);
acb_neg(t, t);
}
if (k != 0)
{
arb_t pi;
arb_init(pi);
arb_const_pi(pi, wp);
arb_addmul_si(acb_imagref(t), pi, (y > 0) ? k : -k, wp);
arb_clear(pi);
}
if (reflect)
{
acb_t v;
acb_init(v);
/* loggamma(x) = log(pi) - lsin(x) - loggamma(2+r-x) - logrf(2+r-x, -r-1) */
acb_hypgeom_log_rising_ui(v, u, -r-1, wp);
acb_log_sin_pi(res, z, wp);
acb_add(res, res, v, wp);
acb_add(res, res, t, wp);
acb_neg(res, res);
arb_const_log_pi(acb_realref(t), wp);
arb_zero(acb_imagref(t));
acb_add(res, res, t, prec);
acb_clear(v);
}
else if (r == 0)
{
acb_set_round(res, t, prec);
}
else if (r > 0)
{
acb_hypgeom_log_rising_ui(res, u, r, wp);
acb_add(res, res, t, prec);
}
else
{
acb_hypgeom_log_rising_ui(res, z, -r, wp);
acb_sub(res, t, res, prec);
}
acb_clear(t);
acb_clear(u);
return 1;
}
}
acb_clear(t);
acb_clear(u);
return 0;
}
void
acb_hypgeom_lgamma(acb_t y, const acb_t x, slong prec)
{
int reflect;
slong r, n, wp;
acb_t t, u, v;
double acc;
if (acb_is_real(x) && arb_is_positive(acb_realref(x)))
{
arb_hypgeom_lgamma(acb_realref(y), acb_realref(x), prec);
arb_zero(acb_imagref(y));
return;
}
if (acb_hypgeom_lgamma_taylor(y, x, prec))
return;
acc = acb_rel_accuracy_bits(x);
acc = FLINT_MAX(acc, 0);
wp = FLINT_MIN(prec, acc + 20);
wp = FLINT_MAX(wp, 2);
wp = wp + FLINT_BIT_COUNT(wp);
acb_hypgeom_gamma_stirling_choose_param(&reflect, &r, &n, x, 1, 0, wp);
acb_init(t);
acb_init(u);
acb_init(v);
if (reflect)
{
/* log gamma(x) = log rf(1-x, r) - log gamma(1-x+r) - log sin(pi x) + log(pi) */
acb_sub_ui(u, x, 1, wp);
acb_neg(u, u);
acb_hypgeom_log_rising_ui(t, u, r, wp);
acb_add_ui(u, u, r, wp);
acb_hypgeom_gamma_stirling_inner(v, u, n, wp);
acb_sub(t, t, v, wp);
acb_log_sin_pi(u, x, wp);
acb_sub(t, t, u, wp);
arb_const_log_pi(acb_realref(u), wp);
arb_zero(acb_imagref(u));
acb_add(y, t, u, wp);
}
else
{
/* log gamma(x) = log gamma(x+r) - log rf(x,r) */
acb_add_ui(t, x, r, wp);
acb_hypgeom_gamma_stirling_inner(u, t, n, wp);
acb_hypgeom_log_rising_ui(t, x, r, wp);
acb_sub(y, u, t, prec);
}
if (!acb_is_finite(y))
acb_indeterminate(y);
acb_clear(t);
acb_clear(u);
acb_clear(v);
}

120
acb_hypgeom/test/t-lgamma.c Normal file
View file

@ -0,0 +1,120 @@
/*
Copyright (C) 2021 Fredrik Johansson
This file is part of Arb.
Arb is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version. See <http://www.gnu.org/licenses/>.
*/
#include "acb_hypgeom.h"
int main()
{
slong iter;
flint_rand_t state;
flint_printf("lgamma....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 10000 * arb_test_multiplier(); iter++)
{
acb_t z, s1, s2, a, b;
slong prec, ebits, prec2;
prec = 2 + n_randint(state, 200);
if (n_randint(state, 10) == 0)
prec = 2 + n_randint(state, 1000);
if (n_randint(state, 10) == 0)
ebits = 100;
else
ebits = 10;
ebits = 2;
prec2 = 2 + n_randint(state, 200);
acb_init(z);
acb_init(s1);
acb_init(s2);
acb_init(a);
acb_init(b);
acb_randtest(z, state, prec, ebits);
acb_randtest(s1, state, prec, 10);
acb_randtest(s2, state, prec, 10);
if (n_randint(state, 2))
{
acb_hypgeom_lgamma(s1, z, prec);
}
else
{
acb_set(s1, z);
acb_hypgeom_lgamma(s1, s1, prec);
}
acb_add_ui(s2, z, 1, prec2);
acb_hypgeom_lgamma(s2, s2, prec2);
acb_log(a, z, prec2);
acb_sub(s2, s2, a, prec2);
if (!acb_overlaps(s1, s2))
{
flint_printf("FAIL\n\n");
flint_printf("prec = %wd\n\n", prec);
flint_printf("z = "); acb_printn(z, 1000, 0); flint_printf("\n\n");
flint_printf("s1 = "); acb_printn(s1, 1000, 0); flint_printf("\n\n");
flint_printf("s2 = "); acb_printn(s2, 1000, 0); flint_printf("\n\n");
acb_sub(s1, s1, s2, prec2);
flint_printf("s1 - s2 = "); acb_printd(s1, 1000); flint_printf("\n\n");
flint_abort();
}
acb_get_mid(a, z);
if (n_randint(state, 2))
{
arf_set_mag(arb_midref(acb_realref(b)), arb_radref(acb_realref(z)));
arf_set_mag(arb_midref(acb_imagref(b)), arb_radref(acb_imagref(z)));
if (n_randint(state, 2))
acb_neg(b, b);
if (n_randint(state, 2))
acb_conj(b, b);
acb_add(a, a, b, prec);
}
acb_hypgeom_lgamma(s2, a, prec);
if (!acb_overlaps(s1, s2))
{
flint_printf("FAIL (2)\n\n");
flint_printf("prec = %wd\n\n", prec);
flint_printf("z = "); acb_printn(z, 1000, 0); flint_printf("\n\n");
flint_printf("a = "); acb_printn(a, 1000, 0); flint_printf("\n\n");
flint_printf("s1 = "); acb_printn(s1, 1000, 0); flint_printf("\n\n");
flint_printf("s2 = "); acb_printn(s2, 1000, 0); flint_printf("\n\n");
acb_sub(s1, s1, s2, prec2);
flint_printf("s1 - s2 = "); acb_printd(s1, 1000); flint_printf("\n\n");
flint_abort();
}
acb_clear(z);
acb_clear(s1);
acb_clear(s2);
acb_clear(a);
acb_clear(b);
}
flint_randclear(state);
flint_cleanup();
flint_printf("PASS\n");
return EXIT_SUCCESS;
}

View file

@ -59,6 +59,8 @@ int arb_hypgeom_gamma_taylor(arb_t res, const arb_t x, int reciprocal, slong pre
void arb_hypgeom_gamma(arb_t y, const arb_t x, slong prec);
void arb_hypgeom_rgamma(arb_t y, const arb_t x, slong prec);
void arb_hypgeom_lgamma(arb_t y, const arb_t x, slong prec);
void arb_hypgeom_pfq(arb_t res, arb_srcptr a, slong p, arb_srcptr b, slong q,
const arb_t z, int regularized, slong prec);

75
arb_hypgeom/lgamma.c Normal file
View file

@ -0,0 +1,75 @@
/*
Copyright (C) 2021 Fredrik Johansson
This file is part of Arb.
Arb is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version. See <http://www.gnu.org/licenses/>.
*/
#include "arb_hypgeom.h"
void arb_hypgeom_gamma_stirling_choose_param(int * reflect, slong * r, slong * n, const arb_t x, int use_reflect, int digamma, slong prec);
int arb_hypgeom_gamma_exact(arb_t res, const arb_t x, int reciprocal, slong prec);
void arb_hypgeom_gamma_stirling_inner(arb_t s, const arb_t z, slong N, slong prec);
void
arb_hypgeom_lgamma_stirling(arb_t y, const arb_t x, slong prec)
{
int reflect;
slong r, n, wp;
arb_t t, u;
double acc;
/* todo: for large x (if exact or accurate enough), increase precision */
acc = arb_rel_accuracy_bits(x);
acc = FLINT_MAX(acc, 0);
wp = FLINT_MIN(prec, acc + 20);
wp = FLINT_MAX(wp, 2);
wp = wp + FLINT_BIT_COUNT(wp);
arb_hypgeom_gamma_stirling_choose_param(&reflect, &r, &n, x, 0, 0, wp);
arb_init(t);
arb_init(u);
/* log(gamma(x)) = log(gamma(x+r)) - log(rf(x,r)) */
arb_init(t);
arb_init(u);
arb_add_ui(t, x, r, wp);
arb_hypgeom_gamma_stirling_inner(u, t, n, wp);
arb_hypgeom_rising_ui_rec(t, x, r, wp);
arb_log(t, t, wp);
arb_sub(y, u, t, prec);
arb_clear(t);
arb_clear(u);
}
void
arb_hypgeom_lgamma(arb_t res, const arb_t x, slong prec)
{
if (!arb_is_positive(x) || !arb_is_finite(x))
{
arb_indeterminate(res);
return;
}
if (arb_hypgeom_gamma_exact(res, x, 0, prec))
{
arb_log(res, res, prec);
return;
}
if (arb_hypgeom_gamma_taylor(res, x, 0, prec))
{
arb_log(res, res, prec);
return;
}
arb_hypgeom_lgamma_stirling(res, x, prec);
}

118
arb_hypgeom/test/t-lgamma.c Normal file
View file

@ -0,0 +1,118 @@
/*
Copyright (C) 2021 Fredrik Johansson
This file is part of Arb.
Arb is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version. See <http://www.gnu.org/licenses/>.
*/
#include "arb_hypgeom.h"
int main()
{
slong iter;
flint_rand_t state;
flint_printf("lgamma....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 10000 * arb_test_multiplier(); iter++)
{
arb_t z, s1, s2, a, b;
slong prec, ebits, prec2;
prec = 2 + n_randint(state, 200);
if (n_randint(state, 10) == 0)
prec = 2 + n_randint(state, 1000);
if (n_randint(state, 10) == 0)
ebits = 100;
else
ebits = 10;
ebits = 2;
prec2 = 2 + n_randint(state, 200);
arb_init(z);
arb_init(s1);
arb_init(s2);
arb_init(a);
arb_init(b);
arb_randtest(z, state, prec, ebits);
arb_randtest(s1, state, prec, 10);
arb_randtest(s2, state, prec, 10);
if (n_randint(state, 2))
{
arb_hypgeom_lgamma(s1, z, prec);
}
else
{
arb_set(s1, z);
arb_hypgeom_lgamma(s1, s1, prec);
}
arb_add_ui(s2, z, 1, prec2);
arb_hypgeom_lgamma(s2, s2, prec2);
arb_log(a, z, prec2);
arb_sub(s2, s2, a, prec2);
if (!arb_overlaps(s1, s2))
{
flint_printf("FAIL\n\n");
flint_printf("prec = %wd\n\n", prec);
flint_printf("z = "); arb_printn(z, 1000, 0); flint_printf("\n\n");
flint_printf("s1 = "); arb_printn(s1, 1000, 0); flint_printf("\n\n");
flint_printf("s2 = "); arb_printn(s2, 1000, 0); flint_printf("\n\n");
arb_sub(s1, s1, s2, prec2);
flint_printf("s1 - s2 = "); arb_printd(s1, 1000); flint_printf("\n\n");
flint_abort();
}
arb_set(a, z);
mag_zero(arb_radref(a));
if (n_randint(state, 2))
{
arf_set_mag(arb_midref(b), arb_radref(z));
if (n_randint(state, 2))
arb_neg(b, b);
arb_add(a, a, b, prec);
}
arb_hypgeom_lgamma(s2, a, prec);
if (!arb_overlaps(s1, s2))
{
flint_printf("FAIL (2)\n\n");
flint_printf("prec = %wd\n\n", prec);
flint_printf("z = "); arb_printn(z, 1000, 0); flint_printf("\n\n");
flint_printf("a = "); arb_printn(a, 1000, 0); flint_printf("\n\n");
flint_printf("s1 = "); arb_printn(s1, 1000, 0); flint_printf("\n\n");
flint_printf("s2 = "); arb_printn(s2, 1000, 0); flint_printf("\n\n");
arb_sub(s1, s1, s2, prec2);
flint_printf("s1 - s2 = "); arb_printd(s1, 1000); flint_printf("\n\n");
flint_abort();
}
arb_clear(z);
arb_clear(s1);
arb_clear(s2);
arb_clear(a);
arb_clear(b);
}
flint_randclear(state);
flint_cleanup();
flint_printf("PASS\n");
return EXIT_SUCCESS;
}

View file

@ -99,16 +99,21 @@ Gamma function
and returns 0. If *reciprocal* is set, the reciprocal gamma function is
computed instead.
.. function:: void acb_hypgeom_gamma(acb_t y, const acb_t x, slong prec)
.. function:: void acb_hypgeom_gamma(acb_t res, const acb_t x, slong prec)
Sets *res* to the gamma function of *x* computed using a default
algorithm choice.
.. function:: void acb_hypgeom_rgamma(acb_t y, const acb_t x, slong prec)
.. function:: void acb_hypgeom_rgamma(acb_t res, const acb_t x, slong prec)
Sets *res* to the reciprocal gamma function of *x* computed using a default
algorithm choice.
.. function:: void acb_hypgeom_lgamma(acb_t res, const acb_t x, slong prec)
Sets *res* to the principal branch of the log-gamma function of *x*
computed using a default algorithm choice.
Convergent series
-------------------------------------------------------------------------------

View file

@ -103,16 +103,20 @@ Gamma function
and returns 0. If *reciprocal* is set, the reciprocal gamma function is
computed instead.
.. function:: void arb_hypgeom_gamma(arb_t y, const arb_t x, slong prec)
.. function:: void arb_hypgeom_gamma(arb_t res, const arb_t x, slong prec)
Sets *res* to the gamma function of *x* computed using a default
algorithm choice.
.. function:: void arb_hypgeom_rgamma(arb_t y, const arb_t x, slong prec)
.. function:: void arb_hypgeom_rgamma(arb_t res, const arb_t x, slong prec)
Sets *res* to the reciprocal gamma function of *x* computed using a default
algorithm choice.
.. function:: void arb_hypgeom_lgamma(arb_t res, const arb_t x, slong prec)
Sets *res* to the log-gamma function of *x* computed using a default
algorithm choice.
Binomial coefficients
-------------------------------------------------------------------------------