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

some further lambertw precision tweaks

Browse source
This commit is contained in:
Fredrik Johansson 2017-04-22 13:59:37 +02:00
parent 8dcf5b383a
commit 919e6179a5
3 changed files with 47 additions and 5 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

13
acb/lambertw.c
View file

@ -263,8 +263,19 @@ void acb_lambertw_main(acb_t res, const acb_t z,
else else
ebits = fmpz_bits(MAG_EXPREF(err)); ebits = 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); acb_get_mag(err, ez1);
mbits = fmpz_bits(MAG_EXPREF(err)); mbits = -MAG_EXP(err);
mbits = FLINT_MAX(mbits, 0);
mbits = FLINT_MIN(mbits, prec);
}
else
{
mbits = 0;
}
kbits = fmpz_bits(k); kbits = fmpz_bits(k);
extraprec = FLINT_MAX(ebits, kbits); extraprec = FLINT_MAX(ebits, kbits);

18
acb/lambertw_check_branch.c
View file

@ -24,9 +24,25 @@ _acb_lambertw_check_branch(const acb_t w, const fmpz_t k, slong prec)
arb_init(b); arb_init(b);
/* t = x sinc(y), v = -cos(y) */ /* t = x sinc(y), v = -cos(y) */
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_sinc(t, acb_imagref(w), prec);
arb_mul(t, t, acb_realref(w), prec);
arb_cos(v, acb_imagref(w), prec); arb_cos(v, acb_imagref(w), prec);
}
arb_mul(t, t, acb_realref(w), prec);
arb_neg(v, v); arb_neg(v, v);
/* u = y / pi, with conjugate relation for k */ /* u = y / pi, with conjugate relation for k */

17
arb/lambertw.c
View file

@ -586,7 +586,7 @@ arb_lambertw(arb_t res, const arb_t x, int flags, slong prec)
else else
{ {
slong k, padding, nextstep, maxstep, *steps; slong k, padding, nextstep, maxstep, *steps;
double rate; double rate, nearm1;
steps = flint_malloc(sizeof(slong) * FLINT_BITS); 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); rate = 2.0 + 1.0 / (1.0 + 0.01 * ebits);
padding = 6 * ebits2; 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; maxstep = 0;
steps[0] = wp; steps[0] = wp;