Threaded Mode | Linear Mode

lyuka · (This post was last modified: 11-07-2020 02:41 PM by lyuka.)

I modified LW0H and LW0 for CAS and used them.

In CAS, both LW0 and LW0H show straightforward convergence,
and in this case shown below, correct to the last digit.
but in HOME, that is, in numerical calculation with 12 decimal digits,
convergence from around the point where the correction value reaches
several tens of ULP (unit in last place) is unnaturally slow.

In the following cases of x = -1.78, the calculation results were all
correct up to the 12th digit in CAS, but not in HOME.

LW0 Rev.1.38 convergence example at x = -1.78 at HOME
[Image: LW0_num-1.38-convergence-at--1.78.png]

[Image: LW0_num-1.38-convergence-at--1.78.png]

LW0 Rev.1.38 convergence example at x = -1.78 in CAS
[Image: LW0_CAS-1.38-convergence-at--1.78.png]

LW0H Rev.1.37 convergence example at x = -1.78 at HOME
[Image: LW0H_num-1.37-convergence-at--1.78.png]

LW0H Rev.1.37 convergence example at x = -1.78 in CAS
[Image: LW0H_CAS-1.37-convergence-at--1.78.png]

Calculation result by WolframAlpha --
W0(-1.78) ≈ 0.089218049856209317716109060765294251305387407793167833503833064… + i1.6256236744277681331520891017303249073577948841857860157186500…

This convergence behavior in HOME may be due to the fact that 12 digits are not large enough, but,
it's similar to the behavior when there is a problem with normalize-and-rounding for each numerical calculation.

Here is the LW0 code for CAS used for testing.

Code:

#CAS

//LW0 - Principal branch of the Lambert W function

//Rev.1.38 (Nov. 7, 2020) (c) Takayuki HOSODA (aka Lyuka)

LW0_CAS_test(x):=

BEGIN

  LOCAL y, p, s, t; 

  LOCAL q, r; 

  LOCAL n;

  n:=0;

  HComplex := 1; // Enable complex result from a real input.

  r := 1 / e;

  IF simplify(x + r) == 0 THEN return -1; END;

  q := e - sqrt(2) - 1;

  s := x + r;

  y := ln(r + sqrt((r + r) * s) + q * s); // approximation near x=-1/e

  PRINT;  // clear terminal

  PRINT("LW0_CAS("+x +")");

  PRINT("y" + n + "=" + y);

  r := 9;

  REPEAT

    n := n + 1;

    p := y;

    q := r;

    t := exp(y);

    s := y * t;

    t := (1 + 0.5 * y) * (s + x) + t;

    IF (0 == t) THEN break; END;

    y := y - (y + 1) * (s - x) / t; // Halley's method

    r := abs(y - p); // correction radius;

    PRINT("r" + n + "=" + r);

  UNTIL 0 == r OR (q < 1 AND q <= r); // convergence check

  PRINT("y=" + p);

  return p;

END;

#END

I'm new to programming in CAS, so I'd appreciate it if you could tell me anything wrong.

Please refer to the page below for details.
Lambert W function and exponent of Lambert W function for HP Prime