The arb_lambertw anomaly near -1/e

[advanpix]

The arb_lambertw is not capable to provide accurate result near x=-1/e point, even with adaptive refinement of working precision.

Code below computes x = -1/e with 100 bits (~30 decimal digits) of precision. Then it calls arb_lambertw, checks relative accuracy of the result and increases working precision if it was not enough.

The anomaly is that arb_lambertw is not able to provide result with accuracy higher than 15 decimal digits even if use thousands of bits as working precision.

    /* LambertW accuracy test near branch cut point -1/e */
    slong digits = 30; /* decimal */
    slong prec   = slong(ceil(digits * log2(10.0))); /* target precision in bits */
    slong i = 0;

    /* x = -1/e */
    arb_t x;
    arb_init(x);
    arb_const_e(x, prec);
    arb_ui_div (x,1,x,prec);
    arb_neg(x,x);

    arb_t w;
    arb_init(w);

    bool fconverged = false;
    bool fstopped   = false;
    for(int wp = 1*prec; !(fconverged || fstopped) ; wp *= 2)
    {
        arb_lambertw(w, x, 0, wp);

        ++i;
        fconverged = (arb_rel_accuracy_bits(w) >= prec);
        fstopped   = (i > 10);
    }

    printf("target precision = %d\narb_rel_accuracy_bits(w) = %d (refinement steps = %d)\n",prec,arb_rel_accuracy_bits(w),i);

    printf("x = ");
    arb_printd(x,digits);
    printf("\n");

    printf("w = ");
    arb_printd(w,digits);
    printf("\n");

    arb_clear(w);
    arb_clear(x);
    flint_cleanup();

Code output:

target precision = 100
arb_rel_accuracy_bits(w) = 1 (refinement steps = 11)
x = -0.367879441171442321595523770161 +/- 8.2148e-31
w = -0.999999999999999701721040753804 +/- inf

Interestingly, if we use true complex-oriented code from acb_lambertw (without switching to arb_lambertw), then accuracy converges quickly and result is good.

I suspect heuristic in arb_lambert function for handling near -1/e arguments is failing somehow.

I use the latest 2.17.0 (on Windows + Intel compiler).

[fredrik-johansson]

(Note: depending on what you intended, you might want to recompute x = -1/e to the new precision inside the loop.)

The non-finite output is expected behavior. The ball for -1/e contains points to the left of -1/e, where W becomes nonreal, so the arb_lambertw result should actually be a nan-containing ball. If you want the complex extension, you have to use acb_lambertw. If you want the image of W(x) evaluated at (x intersected with [-1/e,inf)), you can either use acb_lambertw and extract the real part, or in this case perhaps exploit monotonicity and evaluate at the right endpoint of the interval.

You have the same behavior with arb_sqrt evaluated on [-eps,eps], for example. In that case, there is a workaround: arb_sqrtpos.

I suppose if people find it useful, we could add an option to flags of arb_lambertw for taking the intersection.

[advanpix]

Thank you for quick reply.

@"If you want the image of W(x) evaluated at (x intersected with [-1/e,inf)), you can either use acb_lambertw and extract the real part...."

The acb_lambertw switches to arb_lambertw near the -1/e automatically, so that pure complex code cannot be used now. Is there official way to use pure complex code without automatic switching to arb_lambertw in acb_lambertw?

Of course, I can do this by changing acb_lambertw code itself, but I think it would have sense to have this in Arb officially.

Follow up questions, I see that acb_lambertw switches to arb_lambert by checking the sign of z*e+1:

    /* precompute z*e + 1 */
    arb_const_e(acb_realref(ez1), prec);  
    acb_mul(ez1, ez1, z, prec);
    acb_add_ui(ez1, ez1, 1, prec);    <-- cancellation for z ~ -1/e

    /* use real code when possible */
    if (acb_is_real(z) && arb_is_positive(acb_realref(ez1)) &&
        (fmpz_is_zero(k) ||
        (fmpz_equal_si(k, -1) && arb_is_negative(acb_realref(z)))))
    {
        arb_lambertw(acb_realref(res), acb_realref(z), !fmpz_is_zero(k), prec);
        arb_zero(acb_imagref(res));
    }
    else
    {
        _acb_lambertw(res, z, ez1, k, flags, prec);
    }

So that the code is using arb_lambertw for any points arbitrarily close to -1/e from the right. As you explained above, this results in non-finite and low accuracy result. Maybe it has sense to relax this condition and use arb_lambertw only for sufficiently distant points on the right (so that whole ball is on the right)?

Moreover, do we need to do some extra accuracy check for ze+1 to detect/handle cancellation near z=-1/e? I think, for different prec (e.g. in adaptive refinement) we might get different sign of ze+1 for the same argument z near -1/e.

[fredrik-johansson]

The acb_lambertw switches to arb_lambertw near the -1/e automatically

Yes, but it should do so safely. Do you have a counterexample?

Because of the square root-type singularity, you lose about half the digits close to it; for example the image of W_0 evaluated at [-1/e +/- 10^-n] has radius about 10^(-n/2). This is a mathematical property of the function and not a bug in the implementation.

[advanpix]

Quick summary of the most important points of our discussion (I think we are loosing them):

me: The arb_lambertw gives high error near -1/e (from the right). you: yes, it is normal, since ball includes nan. Use acb_lambertw and cut imaginary part of the result if you want higher accuracy. me: Ok, but this is impossible since acb_lambertw automatically calls arb_lambertw for any points on the right of -1/e.

So, to use pure-complex code to compute lambertw accurately near -1/e, I have to change acb_lambertw (disable it from calling the arb_lambertw). Am I right?

[fredrik-johansson]

No. acb_lambertw should never give a nan-containing ball near -1/e. Here are several examples of acb_lambertw working as intended near -1/e:

>>> ctx.dps = 15
>>> (-acb(-1).exp() + acb("1e-15")).lambertw()
[-0.99999992 +/- 7.07e-9]
>>> (-acb(-1).exp() + acb("+/- 1e-10")).lambertw()
[-1.000 +/- 2.34e-5] + [+/- 2.34e-5]j
>>> ctx.dps = 30
>>> (-acb(-1).exp() + acb("1e-30")).lambertw()
[-0.999999999999998 +/- 5.14e-16]
>>> (-acb(-1).exp() + acb("1e-15 +/- 1e-10")).lambertw()
[-1.000 +/- 2.34e-5] + [+/- 2.34e-5]j
>>> (-acb(-1).exp() + acb("1e-15 +/- 1e-20")).lambertw()
[-0.999999926267 +/- 4.39e-13]
>>> ctx.dps = 50
>>> (-acb(-1).exp() + acb("1e-30")).lambertw()
[-0.99999999999999766835601840287760882 +/- 5.11e-36]
>>> (-acb(-1).exp() + acb("1e-30 +/- 1e-40")).lambertw()
[-0.999999999999997668356018 +/- 5.25e-25]
>>> ctx.dps = 1000
>>> a = -acb(-1).exp() + arb("1e-800")
>>> ctx.dps = 15
>>> a.lambertw()
[-1.000000 +/- 2.78e-8] + [+/- 2.35e-8]j    # not optimal, but correct
>>> a.lambertw()
[-1.000000000000000000000000 +/- 9.64e-26] + [+/- 8.13e-26]j
>>> ctx.dps = 1000
>>> a.lambertw()
[-0.999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999766835601840287579663646393783159912361976370081241157699190355222398995058734265049737000820452224575166497368792663755163357022208938301544460451399444692583003117861458257305665825912607687786636316 +/- 2.21e-601]

Do you have a failing test case?

[fredrik-johansson]

acb_lambertw does a low-precision comparison with -1/e to determine whether to call arb_lambertw, and therefore does not call it with any point to the right of -1/e; only for points isolated from -1/e by about 2^-prec.

[advanpix]

Oh, see now. I was confused by your first reply: "The non-finite output is expected behavior. The ball for -1/e contains points to the left of -1/e, where W becomes non-real, so the arb_lambertw result should actually be a nan-containing ball. "

Now it is clear, thank you for clarification.