From 9916182253deeb400d9b022e30056a4ac9a01b18 Mon Sep 17 00:00:00 2001 From: fredrik Date: Sun, 19 Sep 2021 12:56:03 +0200 Subject: [PATCH] use table for const_euler at low prec --- arb/const_euler.c | 29 ++++++++++++++++++++++++++++- 1 file changed, 28 insertions(+), 1 deletion(-) diff --git a/arb/const_euler.c b/arb/const_euler.c index bbfcc2fb..90fa4279 100644 --- a/arb/const_euler.c +++ b/arb/const_euler.c @@ -11,6 +11,7 @@ #include "arb.h" #include "hypgeom.h" +#include "arb_hypgeom.h" typedef struct { @@ -333,5 +334,31 @@ arb_const_euler_eval(arb_t res, slong prec) euler_bsplit_clear(sum); } -ARB_DEF_CACHED_CONSTANT(arb_const_euler, arb_const_euler_eval) +ARB_DEF_CACHED_CONSTANT(arb_const_euler_brent_mcmillan, arb_const_euler_eval) + +extern const mp_limb_t arb_hypgeom_gamma_tab_limbs[]; + +void +arb_const_euler(arb_t res, slong prec) +{ + if (prec < ARB_HYPGEOM_GAMMA_TAB_PREC - 16) + { + slong exp; + mp_size_t n; + + n = ARB_HYPGEOM_GAMMA_TAB_PREC / FLINT_BITS; + + /* just reading the table is known to give the correct rounding */ + _arf_set_round_mpn(arb_midref(res), &exp, arb_hypgeom_gamma_tab_limbs + n, n, 0, prec, ARF_RND_NEAR); + _fmpz_set_si_small(ARF_EXPREF(arb_midref(res)), exp); + + /* 1/2 ulp error */ + _fmpz_set_si_small(MAG_EXPREF(arb_radref(res)), exp - prec); + MAG_MAN(arb_radref(res)) = MAG_ONE_HALF; + } + else + { + arb_const_euler_brent_mcmillan(res, prec); + } +}