diff --git a/acb/lambertw.c b/acb/lambertw.c index 8b44bf31..fa45f469 100644 --- a/acb/lambertw.c +++ b/acb/lambertw.c @@ -263,8 +263,19 @@ void acb_lambertw_main(acb_t res, const acb_t z, else ebits = fmpz_bits(MAG_EXPREF(err)); - acb_get_mag(err, ez1); - mbits = fmpz_bits(MAG_EXPREF(err)); + if (fmpz_is_zero(k) || (fmpz_is_one(k) && arb_is_negative(acb_imagref(z))) + || (fmpz_equal_si(k, -1) && arb_is_nonnegative(acb_imagref(z)))) + { + acb_get_mag(err, ez1); + mbits = -MAG_EXP(err); + mbits = FLINT_MAX(mbits, 0); + mbits = FLINT_MIN(mbits, prec); + } + else + { + mbits = 0; + } + kbits = fmpz_bits(k); extraprec = FLINT_MAX(ebits, kbits); diff --git a/acb/lambertw_check_branch.c b/acb/lambertw_check_branch.c index 7427964f..d6e546cf 100644 --- a/acb/lambertw_check_branch.c +++ b/acb/lambertw_check_branch.c @@ -24,9 +24,25 @@ _acb_lambertw_check_branch(const acb_t w, const fmpz_t k, slong prec) arb_init(b); /* t = x sinc(y), v = -cos(y) */ - arb_sinc(t, acb_imagref(w), prec); + if (arb_is_exact(acb_imagref(w))) + { + if (arb_is_zero(acb_imagref(w))) + { + arb_one(t); + arb_one(v); + } + else + { + arb_sin_cos(t, v, acb_imagref(w), prec); + arb_div(t, t, acb_imagref(w), prec); + } + } + else + { + arb_sinc(t, acb_imagref(w), prec); + arb_cos(v, acb_imagref(w), prec); + } arb_mul(t, t, acb_realref(w), prec); - arb_cos(v, acb_imagref(w), prec); arb_neg(v, v); /* u = y / pi, with conjugate relation for k */ diff --git a/arb/lambertw.c b/arb/lambertw.c index 47c2f79e..7d5a0211 100644 --- a/arb/lambertw.c +++ b/arb/lambertw.c @@ -586,7 +586,7 @@ arb_lambertw(arb_t res, const arb_t x, int flags, slong prec) else { slong k, padding, nextstep, maxstep, *steps; - double rate; + double rate, nearm1; steps = flint_malloc(sizeof(slong) * FLINT_BITS); @@ -597,6 +597,21 @@ arb_lambertw(arb_t res, const arb_t x, int flags, slong prec) rate = 2.0 + 1.0 / (1.0 + 0.01 * ebits); padding = 6 * ebits2; + /* extra padding near -1/e */ + nearm1 = arf_get_d(arb_midref(w), ARF_RND_DOWN); + if (fabs(nearm1 + 1.0) < 0.01) + { + arf_add_ui(arb_midref(t), arb_midref(w), 1, 30, ARF_RND_DOWN); + + if (arf_is_zero(arb_midref(t))) + padding += prec; + else + { + slong ee = -ARF_EXP(arb_midref(t)); + padding += FLINT_MIN(FLINT_MAX(2 * ee, 0), prec); + } + } + maxstep = 0; steps[0] = wp;