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

incomplete gamma improvement and tests

Browse source
This commit is contained in:
p15-git-acc 2020-09-07 18:41:00 -05:00
parent 961f73ed4e
commit 4ff47db9fb
6 changed files with 65 additions and 116 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

29
acb_dirichlet/platt_c_bound.c
View file

@ -12,33 +12,6 @@
#include "acb_dirichlet.h"
#include "arb_hypgeom.h"
/* Increase precision adaptively. */
static void
_gamma_upper_workaround(arb_t res, const arb_t s, const arb_t z,
int regularized, slong prec)
{
if (!arb_is_finite(s) || !arb_is_finite(z))
{
arb_indeterminate(res);
}
else
{
arb_t x;
slong i;
arb_init(x);
for (i = 0; i < 5; i++)
{
arb_hypgeom_gamma_upper(x, s, z, regularized, prec << i);
if (arb_rel_accuracy_bits(x) > 1)
{
break;
}
}
arb_swap(res, x);
arb_clear(x);
}
}
static void
_arb_pow_si(arb_t res, const arb_t x, slong y, slong prec)
{
@ -175,7 +148,7 @@ _pre_c_p(arb_ptr res, slong sigma, const arb_t h, ulong k, slong prec)
arb_set_si(x, k + l + 1);
arb_mul_2exp_si(x, x, -1);
_gamma_upper_workaround(x, x, u, 0, prec);
arb_hypgeom_gamma_upper(x, x, u, 0, prec);
arb_mul(res + l, res + l, x, prec);
}

31
acb_dirichlet/platt_lemma_A11.c
View file

@ -12,33 +12,6 @@
#include "acb_dirichlet.h"
#include "arb_hypgeom.h"
/* Increase precision adaptively. */
static void
_gamma_upper_workaround(arb_t res, const arb_t s, const arb_t z,
int regularized, slong prec)
{
if (!arb_is_finite(s) || !arb_is_finite(z))
{
arb_indeterminate(res);
}
else
{
arb_t x;
slong i;
arb_init(x);
for (i = 0; i < 5; i++)
{
arb_hypgeom_gamma_upper(x, s, z, regularized, prec << i);
if (arb_rel_accuracy_bits(x) > 1)
{
break;
}
}
arb_swap(res, x);
arb_clear(x);
}
}
static void
_platt_lemma_A11_X(arb_t out,
const arb_t t0, const arb_t h, slong B, const arb_t beta, slong prec)
@ -90,7 +63,7 @@ _platt_lemma_A11_Y(arb_t out,
arb_sqr(x4, x4, prec);
arb_mul_2exp_si(x4, x4, -3);
_gamma_upper_workaround(x5, x3, x4, 0, prec);
arb_hypgeom_gamma_upper(x5, x3, x4, 0, prec);
arb_mul(out, x1, x2, prec);
arb_mul(out, out, x5, prec);
@ -131,7 +104,7 @@ _platt_lemma_A11_Z(arb_t out,
arb_sqr(x4, x4, prec);
arb_mul_2exp_si(x4, x4, -1);
_gamma_upper_workaround(x5, x3, x4, 0, prec);
arb_hypgeom_gamma_upper(x5, x3, x4, 0, prec);
arb_mul(out, x1, x2, prec);
arb_mul(out, out, x5, prec);

30
acb_dirichlet/platt_lemma_A5.c
View file

@ -12,34 +12,6 @@
#include "acb_dirichlet.h"
#include "arb_hypgeom.h"
/* Increase precision adaptively. */
static void
_gamma_upper_workaround(arb_t res, const arb_t s, const arb_t z,
int regularized, slong prec)
{
if (!arb_is_finite(s) || !arb_is_finite(z))
{
arb_indeterminate(res);
}
else
{
arb_t x;
slong i;
arb_init(x);
for (i = 0; i < 5; i++)
{
arb_hypgeom_gamma_upper(x, s, z, regularized, prec << i);
if (arb_rel_accuracy_bits(x) > 1)
{
break;
}
}
arb_swap(res, x);
arb_clear(x);
}
}
/* Lemma A5 requires B > h*sqrt(k) */
static int
_platt_lemma_A5_constraint(slong B, const arb_t h, slong k, slong prec)
@ -106,7 +78,7 @@ acb_dirichlet_platt_lemma_A5(arb_t out, slong B, const arb_t h,
arb_mul_2exp_si(x6, b, -1);
_gamma_upper_workaround(x6, x6, a, 0, prec);
arb_hypgeom_gamma_upper(x6, x6, a, 0, prec);
arb_mul(out, x4, x5, prec);
arb_mul(out, out, x6, prec);

31
acb_dirichlet/platt_ws_interpolation.c
View file

@ -13,33 +13,6 @@
#include "arb_hypgeom.h"
#include "arb_poly.h"
/* Increase precision adaptively. */
static void
_gamma_upper_workaround(arb_t res, const arb_t s, const arb_t z,
int regularized, slong prec)
{
if (!arb_is_finite(s) || !arb_is_finite(z))
{
arb_indeterminate(res);
}
else
{
arb_t x;
slong i;
arb_init(x);
for (i = 0; i < 5; i++)
{
arb_hypgeom_gamma_upper(x, s, z, regularized, prec << i);
if (arb_rel_accuracy_bits(x) > 1)
{
break;
}
}
arb_swap(res, x);
arb_clear(x);
}
}
static void
_arb_div_si_si(arb_t res, slong x, slong y, slong prec)
{
@ -218,7 +191,7 @@ _platt_bound_C3_Y(arb_t res, const arb_t t0, slong A, const arb_t H,
arb_div(g2, g2, H, prec);
arb_sqr(g2, g2, prec);
arb_mul_2exp_si(g2, g2, -1);
_gamma_upper_workaround(g, g1, g2, 0, prec);
arb_hypgeom_gamma_upper(g, g1, g2, 0, prec);
/* res = a*b*A*H*g */
arb_mul_si(res, H, A, prec);
@ -261,7 +234,7 @@ _platt_bound_C3_Z(arb_t res, const arb_t t0, slong A, const arb_t H,
arb_div(g2, t0, H, prec);
arb_sqr(g2, g2, prec);
arb_mul_2exp_si(g2, g2, -1);
_gamma_upper_workaround(g, g1, g2, 0, prec);
arb_hypgeom_gamma_upper(g, g1, g2, 0, prec);
/* res = a*b*A*g */
arb_mul_si(res, g, A, prec);

4
acb_hypgeom/gamma_upper.c
View file

@ -398,7 +398,9 @@ acb_hypgeom_gamma_upper(acb_t res, const acb_t s, const acb_t z, int regularized
c = 2.391097143;
ds = arf_get_d(arb_midref(acb_realref(s)), ARF_RND_DOWN);
dx = arf_get_d(arb_midref(acb_realref(z)), ARF_RND_DOWN);
if (pow(dx - c, 8) > pow(a*ds, 8) + pow(b*prec, 8))
if (dx - c > a*ds + b*prec ||
(dx - c > FLINT_MAX(a*ds, b*prec) &&
pow(dx - c, 8) > pow(a*ds, 8) + pow(b*prec, 8)))
{
acb_hypgeom_gamma_upper_asymp(res, s, z, regularized, prec);
return;

56
acb_hypgeom/test/t-gamma_upper.c
View file

@ -11,6 +11,28 @@
#include "acb_hypgeom.h"
static void
_accuracy_regression_test(const acb_t s, const acb_t z,
int regularized, slong prec, slong issue, slong accuracy)
{
acb_t g;
acb_init(g);
acb_hypgeom_gamma_upper(g, s, z, regularized, prec);
if (acb_rel_accuracy_bits(g) < accuracy)
{
flint_printf("FAIL: accuracy regression in issue #%ld\n\n", issue);
flint_printf("prec = %d\n\n", prec);
flint_printf("regularized = %d\n\n", regularized);
flint_printf("s = "); acb_printd(s, 30); flint_printf("\n\n");
flint_printf("z = "); acb_printd(z, 30); flint_printf("\n\n");
flint_printf("g = "); acb_printd(g, 30); flint_printf("\n\n");
flint_abort();
}
acb_clear(g);
}
int main()
{
slong iter;
@ -211,6 +233,40 @@ int main()
acb_clear(u);
}
/* Accuracy regression tests. */
{
acb_t s, z;
slong prec, issue, accuracy;
acb_init(s);
acb_init(z);
issue = 166;
prec = 165;
accuracy = 100;
acb_zero(s);
acb_set_si(z, 110);
_accuracy_regression_test(s, z, 2, prec, issue, accuracy);
issue = 276;
prec = 300;
accuracy = 100;
acb_set_ui(s, 357);
acb_set_ui(z, 356);
_accuracy_regression_test(s, z, 0, prec, issue, accuracy);
arb_set_str(acb_realref(s), "356.123", prec);
arb_set_str(acb_realref(z), "356.456", prec);
_accuracy_regression_test(s, z, 0, prec, issue, accuracy);
arb_set_str(acb_realref(s), "357.123", prec);
arb_set_str(acb_realref(z), "356.456", prec);
_accuracy_regression_test(s, z, 0, prec, issue, accuracy);
arb_set_str(acb_realref(s), "357.456", prec);
arb_set_str(acb_realref(z), "356.123", prec);
_accuracy_regression_test(s, z, 0, prec, issue, accuracy);
acb_clear(s);
acb_clear(z);
}
flint_randclear(state);
flint_cleanup();
flint_printf("PASS\n");