mirror of
https://github.com/vale981/arb
synced 2025-03-04 17:01:40 -05:00
complex taylor gamma; bugfixes
This commit is contained in:
fredrik
parent
5eeb61ca2e
commit
01c219f876
5 changed files with 445 additions and 2 deletions
|
|
@ -167,6 +167,9 @@ acb_hypgeom_gamma(acb_t y, const acb_t x, slong prec)
|
|||
return;
|
||||
}
|
||||
|
||||
if (acb_hypgeom_gamma_taylor(y, x, 0, prec))
|
||||
return;
|
||||
|
||||
acb_hypgeom_gamma_stirling(y, x, 0, prec);
|
||||
}
|
||||
|
||||
|
|
@ -180,6 +183,9 @@ acb_hypgeom_rgamma(acb_t y, const acb_t x, slong prec)
|
|||
return;
|
||||
}
|
||||
|
||||
if (acb_hypgeom_gamma_taylor(y, x, 1, prec))
|
||||
return;
|
||||
|
||||
acb_hypgeom_gamma_stirling(y, x, 1, prec);
|
||||
}
|
||||
|
||||
|
|
|
|||
295
acb_hypgeom/gamma_taylor.c
Normal file
295
acb_hypgeom/gamma_taylor.c
Normal file
|
|
@ -0,0 +1,295 @@
|
|||
/*
|
||||
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"
|
||||
#include "acb_hypgeom.h"
|
||||
|
||||
static void
|
||||
evaluate_rect(acb_t res, const short * term_prec, slong len, const acb_t x, slong prec)
|
||||
{
|
||||
slong i, j, m, r, n1, n2;
|
||||
acb_ptr xs;
|
||||
acb_t s, t;
|
||||
arb_struct c[17];
|
||||
|
||||
m = n_sqrt(len) + 1;
|
||||
m = FLINT_MIN(m, 16);
|
||||
r = (len + m - 1) / m;
|
||||
|
||||
xs = _acb_vec_init(m + 1);
|
||||
acb_init(s);
|
||||
acb_init(t);
|
||||
_acb_vec_set_powers(xs, x, m + 1, prec);
|
||||
|
||||
acb_zero(res);
|
||||
|
||||
for (i = r - 1; i >= 0; i--)
|
||||
{
|
||||
n1 = m * i;
|
||||
n2 = FLINT_MIN(len, n1 + m);
|
||||
|
||||
for (j = n1; j < n2; j++)
|
||||
{
|
||||
if (j == 0)
|
||||
{
|
||||
arb_init(c);
|
||||
arb_one(c);
|
||||
}
|
||||
else
|
||||
{
|
||||
if (!_arb_hypgeom_gamma_coeff_shallow(arb_midref(c + j - n1), arb_radref(c + j - n1), j, term_prec[j]))
|
||||
flint_abort();
|
||||
}
|
||||
}
|
||||
|
||||
arb_dot(acb_realref(s), NULL, 0, acb_realref(xs), 2, c, 1, n2 - n1, prec);
|
||||
arb_dot(acb_imagref(s), NULL, 0, acb_imagref(xs), 2, c, 1, n2 - n1, prec);
|
||||
|
||||
#if 0
|
||||
acb_set_round(t, xs + m, term_prec[n1]);
|
||||
acb_mul(res, res, t, term_prec[n1]);
|
||||
acb_add(res, res, s, term_prec[n1]);
|
||||
#else
|
||||
acb_mul(res, res, xs + m, term_prec[n1]);
|
||||
acb_add(res, res, s, term_prec[n1]);
|
||||
#endif
|
||||
}
|
||||
|
||||
_acb_vec_clear(xs, m + 1);
|
||||
acb_clear(s);
|
||||
acb_clear(t);
|
||||
}
|
||||
|
||||
/* Bound requires: |u| <= 20, N <= 10000, N != (1443, 2005, 9891). */
|
||||
static void
|
||||
error_bound(mag_t err, const acb_t u, slong N)
|
||||
{
|
||||
mag_t t;
|
||||
mag_init(t);
|
||||
|
||||
acb_get_mag(t, u);
|
||||
|
||||
if (N >= 1443 || mag_cmp_2exp_si(t, 4) > 0)
|
||||
{
|
||||
mag_inf(err);
|
||||
}
|
||||
else
|
||||
{
|
||||
mag_pow_ui(err, t, N);
|
||||
mag_mul_2exp_si(err, err, arb_hypgeom_gamma_coeffs[N].exp);
|
||||
|
||||
if (mag_cmp_2exp_si(t, -1) > 0)
|
||||
mag_mul(err, err, t);
|
||||
else
|
||||
mag_mul_2exp_si(err, err, -1);
|
||||
|
||||
mag_mul_2exp_si(err, err, 3);
|
||||
|
||||
if (mag_cmp_2exp_si(err, -8) > 0)
|
||||
mag_inf(err);
|
||||
}
|
||||
|
||||
mag_clear(t);
|
||||
}
|
||||
|
||||
static double
|
||||
want_taylor(double x, double y, slong prec)
|
||||
{
|
||||
if (y < 0.0) y = -y;
|
||||
if (x < 0.0) x = -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;
|
||||
}
|
||||
|
||||
int
|
||||
acb_hypgeom_gamma_taylor(acb_t res, const acb_t z, int reciprocal, slong prec)
|
||||
{
|
||||
acb_t s, u;
|
||||
int success;
|
||||
double dua, dub, du2, log2u;
|
||||
slong i, r, n, wp, tail_bound, goal;
|
||||
short term_prec[ARB_HYPGEOM_GAMMA_TAB_NUM];
|
||||
mag_t err;
|
||||
|
||||
if (!acb_is_finite(z) ||
|
||||
arf_cmp_2exp_si(arb_midref(acb_imagref(z)), 4) >= 0 ||
|
||||
arf_cmp_2exp_si(arb_midref(acb_realref(z)), 10) >= 0)
|
||||
{
|
||||
return 0;
|
||||
}
|
||||
|
||||
dua = arf_get_d(arb_midref(acb_realref(z)), ARF_RND_UP);
|
||||
dub = arf_get_d(arb_midref(acb_imagref(z)), ARF_RND_UP);
|
||||
dub = fabs(dub);
|
||||
|
||||
if (!want_taylor(dua, dub, prec))
|
||||
return 0;
|
||||
|
||||
if (dua >= 0.0)
|
||||
r = (slong) (dua + 0.5);
|
||||
else
|
||||
r = -(slong) (-dua + 0.5);
|
||||
|
||||
acb_init(s);
|
||||
acb_init(u);
|
||||
mag_init(err);
|
||||
|
||||
success = 0;
|
||||
|
||||
/* Argument reduction: u = z - r */
|
||||
acb_sub_si(u, z, r, 2 * prec + 10);
|
||||
dua -= r;
|
||||
|
||||
goal = acb_rel_accuracy_bits(u);
|
||||
|
||||
/* not designed for wide intervals (yet) */
|
||||
if (goal < 8)
|
||||
{
|
||||
success = 0;
|
||||
goto cleanup;
|
||||
}
|
||||
|
||||
goal = FLINT_MIN(goal, prec - MAG_BITS) + MAG_BITS;
|
||||
goal = FLINT_MAX(goal, 5);
|
||||
goal = goal + 5;
|
||||
wp = goal + 4 + FLINT_BIT_COUNT(FLINT_ABS(r));
|
||||
|
||||
if (wp > ARB_HYPGEOM_GAMMA_TAB_PREC)
|
||||
{
|
||||
success = 0;
|
||||
goto cleanup;
|
||||
}
|
||||
|
||||
if (!want_taylor(r, dub, goal))
|
||||
{
|
||||
success = 0;
|
||||
goto cleanup;
|
||||
}
|
||||
|
||||
du2 = dua * dua + dub * dub;
|
||||
|
||||
if (du2 > 1e-8)
|
||||
{
|
||||
log2u = 0.5 * mag_d_log_upper_bound(du2) * 1.4426950408889634074 * (1 + 1e-14);
|
||||
}
|
||||
else
|
||||
{
|
||||
slong aexp, bexp;
|
||||
|
||||
aexp = arf_cmpabs_2exp_si(arb_midref(acb_realref(u)), -wp) >= 0 ? ARF_EXP(arb_midref(acb_realref(u))) : -wp;
|
||||
bexp = arf_cmpabs_2exp_si(arb_midref(acb_imagref(u)), -wp) >= 0 ? ARF_EXP(arb_midref(acb_imagref(u))) : -wp;
|
||||
log2u = FLINT_MAX(aexp, bexp) + 1;
|
||||
}
|
||||
|
||||
term_prec[0] = wp;
|
||||
n = 0;
|
||||
|
||||
for (i = 1; i < ARB_HYPGEOM_GAMMA_TAB_NUM; i++)
|
||||
{
|
||||
tail_bound = arb_hypgeom_gamma_coeffs[i].exp + i * log2u + 5;
|
||||
|
||||
if (tail_bound <= -goal)
|
||||
{
|
||||
n = i;
|
||||
break;
|
||||
}
|
||||
|
||||
term_prec[i] = FLINT_MIN(FLINT_MAX(wp + tail_bound, 2), wp);
|
||||
|
||||
if (term_prec[i] > arb_hypgeom_gamma_coeffs[i].nlimbs * FLINT_BITS)
|
||||
{
|
||||
success = 0;
|
||||
goto cleanup;
|
||||
}
|
||||
}
|
||||
|
||||
if (n != 0)
|
||||
error_bound(err, u, n);
|
||||
|
||||
if (n == 0 || mag_is_inf(err))
|
||||
{
|
||||
success = 0;
|
||||
goto cleanup;
|
||||
}
|
||||
|
||||
evaluate_rect(s, term_prec, n, u, wp);
|
||||
acb_add_error_mag(s, err);
|
||||
|
||||
if (r == 0 || r == 1)
|
||||
{
|
||||
if (r == 0)
|
||||
acb_mul(s, s, u, wp);
|
||||
|
||||
if (reciprocal)
|
||||
{
|
||||
acb_set_round(res, s, prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
acb_one(u);
|
||||
acb_div(res, u, s, prec);
|
||||
}
|
||||
}
|
||||
else if (r >= 2)
|
||||
{
|
||||
acb_add_ui(u, u, 1, wp);
|
||||
acb_hypgeom_rising_ui_rec(u, u, r - 1, wp);
|
||||
|
||||
if (reciprocal)
|
||||
acb_div(res, s, u, prec);
|
||||
else
|
||||
acb_div(res, u, s, prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
/* gamma(x) = (-1)^r / (rgamma(1+x-r)*rf(1+r-x,-r)*(x-r)) */
|
||||
/* 1/gamma(x) = (-1)^r * rgamma(1+x-r) * rf(1+r-x,-r) * (x-r) */
|
||||
|
||||
acb_neg(res, z);
|
||||
acb_add_si(res, res, 1 + r, wp);
|
||||
acb_hypgeom_rising_ui_rec(res, res, -r, wp);
|
||||
acb_mul(u, res, u, wp);
|
||||
|
||||
if (reciprocal)
|
||||
{
|
||||
acb_mul(res, s, u, prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
acb_mul(u, s, u, wp);
|
||||
acb_inv(res, u, prec);
|
||||
}
|
||||
|
||||
if (r % 2)
|
||||
acb_neg(res, res);
|
||||
}
|
||||
|
||||
success = 1;
|
||||
|
||||
cleanup:
|
||||
acb_clear(s);
|
||||
acb_clear(u);
|
||||
mag_clear(err);
|
||||
|
||||
return success;
|
||||
}
|
||||
|
||||
139
acb_hypgeom/test/t-gamma_taylor.c
Normal file
139
acb_hypgeom/test/t-gamma_taylor.c
Normal file
|
|
@ -0,0 +1,139 @@
|
|||
/*
|
||||
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("gamma_taylor....");
|
||||
fflush(stdout);
|
||||
|
||||
flint_randinit(state);
|
||||
|
||||
for (iter = 0; iter < 10000 * arb_test_multiplier(); iter++)
|
||||
{
|
||||
acb_t x, s1, s2, a, b;
|
||||
slong prec, ebits, prec2;
|
||||
int success, success2, alias, reciprocal;
|
||||
|
||||
if (n_randint(state, 10) == 0)
|
||||
prec = 2 + n_randint(state, 4000);
|
||||
else
|
||||
prec = 2 + n_randint(state, 300);
|
||||
|
||||
if (n_randint(state, 10) == 0)
|
||||
ebits = 100;
|
||||
else
|
||||
ebits = 10;
|
||||
|
||||
prec2 = prec + 1 + n_randint(state, 30);
|
||||
|
||||
acb_init(x);
|
||||
acb_init(s1);
|
||||
acb_init(s2);
|
||||
acb_init(a);
|
||||
acb_init(b);
|
||||
|
||||
acb_randtest(x, state, prec, ebits);
|
||||
acb_randtest(s1, state, prec, 10);
|
||||
acb_randtest(s2, state, prec, 10);
|
||||
alias = n_randint(state, 2);
|
||||
reciprocal = n_randint(state, 2);
|
||||
|
||||
if (alias)
|
||||
{
|
||||
success = acb_hypgeom_gamma_taylor(s1, x, reciprocal, prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
acb_set(s1, x);
|
||||
success = acb_hypgeom_gamma_taylor(s1, s1, reciprocal, prec);
|
||||
}
|
||||
|
||||
if (success)
|
||||
{
|
||||
/* printf("%ld\n", iter); */
|
||||
|
||||
/* Compare with Stirling series algorithm. */
|
||||
acb_hypgeom_gamma_stirling(s2, x, reciprocal, prec);
|
||||
|
||||
if (!acb_overlaps(s1, s2))
|
||||
{
|
||||
flint_printf("FAIL\n\n");
|
||||
flint_printf("prec = %wd\n\n", prec);
|
||||
flint_printf("x = "); acb_printn(x, 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();
|
||||
}
|
||||
|
||||
/* Compare with different level of precision. */
|
||||
success2 = acb_hypgeom_gamma_taylor(s2, x, reciprocal, prec2);
|
||||
|
||||
if (success2 && !acb_overlaps(s1, s2))
|
||||
{
|
||||
flint_printf("FAIL (2)\n\n");
|
||||
flint_printf("prec = %wd\n\n", prec);
|
||||
flint_printf("x = "); acb_printn(x, 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_printn(s1, 1000, 0); flint_printf("\n\n");
|
||||
flint_abort();
|
||||
}
|
||||
|
||||
acb_get_mid(a, x);
|
||||
|
||||
if (n_randint(state, 2))
|
||||
{
|
||||
arf_set_mag(arb_midref(acb_realref(b)), arb_radref(acb_realref(x)));
|
||||
arf_set_mag(arb_midref(acb_imagref(b)), arb_radref(acb_imagref(x)));
|
||||
|
||||
if (n_randint(state, 2))
|
||||
acb_neg(b, b);
|
||||
if (n_randint(state, 2))
|
||||
acb_conj(b, b);
|
||||
|
||||
acb_add(a, a, b, prec2);
|
||||
}
|
||||
|
||||
success2 = acb_hypgeom_gamma_taylor(s2, a, reciprocal, prec2);
|
||||
|
||||
if (success2 && !acb_overlaps(s1, s2))
|
||||
{
|
||||
flint_printf("FAIL (3)\n\n");
|
||||
flint_printf("prec = %wd\n\n", prec);
|
||||
flint_printf("x = "); acb_printn(x, 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_printn(s1, 1000, 0); flint_printf("\n\n");
|
||||
flint_abort();
|
||||
}
|
||||
}
|
||||
|
||||
acb_clear(x);
|
||||
acb_clear(s1);
|
||||
acb_clear(s2);
|
||||
acb_clear(a);
|
||||
acb_clear(b);
|
||||
}
|
||||
|
||||
flint_randclear(state);
|
||||
flint_cleanup();
|
||||
flint_printf("PASS\n");
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
|
|
@ -4232,7 +4232,7 @@ _arb_hypgeom_gamma_coeff_shallow(arf_t c, mag_t err, slong i, slong prec)
|
|||
|
||||
if (err != NULL)
|
||||
{
|
||||
MAG_EXP(err) = exp - term_limbs * FLINT_BITS;
|
||||
MAG_EXP(err) = exp - term_limbs * FLINT_BITS + 1;
|
||||
MAG_MAN(err) = MAG_ONE_HALF;
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -340,7 +340,10 @@ arb_hypgeom_gamma_taylor(arb_t res, const arb_t x, int reciprocal, slong prec)
|
|||
/* Nearest (roughly) integer to x, to use as shift for argument reduction
|
||||
to move to the interval [-0.5,0.5]. It's OK that dx is approximate so
|
||||
that the reduced argument will actually lie in [-0.5-eps,0.5+eps]. */
|
||||
r = (slong) (dx + 0.5);
|
||||
if (dx >= 0.0)
|
||||
r = (slong) (dx + 0.5);
|
||||
else
|
||||
r = -(slong) (-dx + 0.5);
|
||||
|
||||
/* Tuning cutoff. */
|
||||
if (prec >= 40)
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue