(28/48/50) Lambert W Function

[John Keith]

NOTE: The first and last programs have been obsoleted by the program in post #22, which should be used in most cases. That program gives accurate results over the entire complex plane for any branch, positive or negative, but accuracy is limited for values very close to the branch point -1/e.

For very accurate results near the branch point in branches 0 and -1, use the program in post #28, which is a translation of Albert Chan's HP-71 program in post # 19.

The following text refers to the (mostly obsolete) programs in this post only.

---------------------------------------------------------------

This program computes the Lambert W function W(z) for branches 0 and -1. For branch -1, the program covers the real-valued range -1/e < z < 0. For branch 0, the program covers the entire complex plane. For most values, the results are accurate to within 2 ULPs. For values very close to 1/e, the results are approximate. The program expects the branch (0 or -1) on level 2 and z on level 1.

Some notes and background information:

The gold standard for Lambert W programs on the HP 50g is Egan Ford's SPECIAL.LIB. It is very accurate for almost any input value, and can return results for any branch. It does have some problems, however. It is an HPGCC program and only works on the 49g+ and 50g. It requires two large programs on the calculators SD card. Also, it is rather slow. I believe that the slowness is caused by the need to fetch the executable from the SD card, as one would expect the HPGCC code to be quite fast.

The program described here is comparatively limited since no branches other than 0 and -1 are covered, and then only the real-valued range of branch -1. However, it is about 5 times as fast as SPECIAL.LIB and about 9 times as fast as programs using the calculator's built-in solvers. It also runs on all RPL calculators.

If anyone here has an idea of how to get an accurate estimate for other branches, or to improve accuracy around the branch point at -1/e, please chime in!

The code that provides the initial estimate for branch 0 comes from page 8 of Iacono and Boyd. The estimation code for branch -1 comes from formula 27 on page 12 of https://www.sciencedirect.com/science/ar...via%3Dihub. The recurrence used in the main loop occurs in both papers and is mentioned in Wikipedia. Many additional references here.

Code:


\<< 1. EXP ROT ROT SWAP
  IF NOT                                             @ Branch 0
  THEN DUP ABS .00000001
    IF \>=
    THEN DUP ROT * 1. + \v/ DUP 1.14956131 * 1. +    @ Estimate for
      SWAP 1. + LN .4549574 * 1. + / LN 2.036 * 1. - @ branch 0, |z| > 1E-8
    ELSE DUP ROT DROP                                @ DUP If |z| < 1E-8
    END
  ELSE                                               @ Branch -1
    IF DUP -.25 \<=                                  @ Estimate for
    THEN DUP -1. SWAP 4. ROLL * 1. + 2. * \v/ -      @ z <= 0.25
    ELSE DUP ROT DROP NEG LN DUP NEG LN -            @ For z > 0.25
    END
  END 1. 3.                                          @ Main loop
  START DUP2 / LN 1. + SWAP DUP 1. + / *
  NEXT SWAP DROP
\>>

Also, a streamlined version for the 49 and 50, 25 bytes smaller:

Code:


\<< 1. EXP UNROT SWAP
  IF NOT
  THEN DUP ABS .00000001 \>=
    { DUP ROT * 1. + \v/ DUP 1.14956131 * 1. +
      SWAP 1. + LN .4549574 * 1. + / LN 2.036 * 1. - }
    { NIP DUP } IFTE
  ELSE DUP -.25 \<=
    { DUP -1. SWAP 4. ROLL * 1. + 2. * \v/ - }
    { NIP DUP NEG LN DUP NEG LN - } IFTE
  END 1. 3.
  START DUP2 / LN 1. + SWAP DUP 1. + / *
  NEXT NIP
\>>

Now, here are a couple of experimental programs for those interested in further exploration. The first one is a variation of the recurrence from the main program using the LongFloat Library. The program expects z on level 2 and an estimate (such as the result of the main program) on level 1. It is rather slow, about 22 seconds on my 50g but reasonable on emulators. The value of c assumes the default value of DIGITS which is 50. Removing the R\<-\->F at the end will leave a LongFloat number (type 27) on the stack.

Code:


\<< R\<-\->F SWAP R\<-\->F SWAP
    1. R\<-\->F 1.E-24 R\<-\->F \-> f1 c
  \<<
    DO DUP2 FDIV FLN f1 FADD OVER DUP f1 FADD FDIV FMULT
    UNTIL DUP ROT FSUB FABS c F.LE
    END NIP R\<-\->F
  \>>
\>>

The last program implements the Halley recurrence shown at the Wikipedia link and in the Lóczi paper. The reason for including it is that the recurrence used in the main program works poorly for branches other than branch 0 and the real-valued part of branch -1. It tends to "drag" the values back to branch 0 even if given an accurate estimate for the value in another branch. The Halley recurrence converges rapidly even if the estimate is not very close, although it can be fooled if the estimate is vague. As above, the program expects z on level 2 and an estimate (user-supplied) on level 1. For HP-28 and 48, change PICK3 to 3 PICK.

Code:


\<<
  WHILE DUP2 DUP EXP
    DUP ROT *
    ROT -
    PICK3 1. + ROT *
    DUP2 / ABS .00000000005 > @ 5E-11
  REPEAT PICK3 2. + PICK3 *
    4. PICK 2. * 2. + /
    - / -
  END DROP2
\>>