(28/48/50) Lambert W Function

[John Keith]

Here is the current version of the Lambert W program for all RPL calculators. This version covers all representable numbers, real or complex, in any branch. In almost all cases the result will be accurate within 2 ULP's, except for a circle of radius .001 around -1/e.

This program replaces the first and last programs in post #1. As before, level 2 should have an integer representing the branch, and level 1 should have the real or complex number z.

Code:


\<<
  IF DUP ABS                                        @ z not 0?
  THEN SWAP
    IF DUP                                          @ k not 0?
    THEN
      IF DUP -1 SAME                                @ If branch -1
      THEN
        IF OVER RE DUP 0 < SWAP -2 > AND
        3 PICK IM DUP 0 \>= SWAP .2 < AND AND       @ Near -1/e?
        THEN DROP
          IF DUP RE -.25 >                          @ Estimate for branch -1
          THEN DUP NEG LN DUP NEG LN -              @ near -1/e
          ELSE DUP -1. SWAP 1 EXP * 1 +
            2 * \v/ -
          END 6 0                                   @ Do 6 iterations
        ELSE 1                                      @ Use complex code
        END
      ELSE 1                                        @ Use complex code
      END
    ELSE DROP DUP ABS .00000001
      IF \>=                                        @ If |z| >= 1E-8
      THEN DUP 1 EXP * 1 + \v/ DUP 1.14956131 *     @ Estimate for branch 0
        1 + SWAP 1 + LN .4549574 * 1 + /
        LN 2.036 * 1 - 3                            @ Do 3 iterations
      ELSE DUP 1                                    @ DUP if |z| < 1E-8
      END 0
    END
    IF NOT
    THEN 1 SWAP                                     @ Iteration for branch 0
      START DUP2 / LN 1 + SWAP DUP 1 + / *          @ and branch -1 near -1/e
      NEXT
    ELSE OVER SWAP (0, 6.28318530718) *             @ Estimate for
      SWAP LN + DUP LN DUP2 DUP2 - ROT ROT SWAP / + @ complex branches
      ROT ROT DUP 2 - * SWAP SQ 2 * / +
      1 3
      START DUP2 DUP EXP DUP ROT * ROT -            @ Halley iteration
        3 PICK 1 + ROT * 3 PICK 2 + 3 PICK *
        4 PICK 2 * 2 + / - / -
      NEXT
    END
  ELSE                                              @ z = 0
    IF OVER                                         @ If branch not 0, return
    THEN INV NEG                                    @ negative infinity
    END
  END SWAP DROP
\>>

Also, a shorter version for the HP 49 and 50 using UNROT, PICK3 and IFTE. The program can be used in approximate or exact mode, but the result will be approximate. Execution time on my 50g ranges from 110 to 490 ms depending on input.

Code:


\<< \->NUM
  IF DUP ABS
  THEN SWAP IP
    IF DUP
    THEN DUP -1. SAME
      { OVER RE DUP 0. < SWAP -2. > AND
        PICK3 IM DUP 0. \>= SWAP .2 < AND AND
        { DROP DUP RE -.25 >
          { DUP NEG LN DUP NEG LN - }
          { DUP -1. SWAP 1. EXP * 1. + 2. * \v/ -
          } IFTE 6. 0.
        } 1. IFTE
      } 1. IFTE
    ELSE DROP DUP ABS .00000001 \>=
      { DUP 1. EXP * 1. + \v/ DUP 1.14956131 * 1. +
        SWAP 1. + LN .4549574 * 1. + / LN 2.036 * 1. - 3.
      }
      { DUP 1. } IFTE 0.
    END NOT
    { 1. SWAP
    START DUP2 / LN 1. + SWAP DUP 1. + / *
    NEXT
    }
    { OVER SWAP (0.,6.28318530718) * SWAP LN +
      DUP LN DUP2 DUP2 - UNROT SWAP / + UNROT
      DUP 2. - * SWAP SQ 2. * / +
      1. 3.
      START DUP2 DUP EXP DUP ROT * ROT -
        PICK3 1. + ROT * PICK3 2. + PICK3 *
        4. PICK 2. * 2. + / - / -
      NEXT
     } IFTE
  ELSE OVER { INV NEG } IFT
  END SWAP DROP
\>>