# LambertW.frink

``` use ArbitraryPrecision.frink use sqrtWayne.frink // This calculates the value of the Lambert W function, also // sometimes called the Product Log function.  This calculates the principal // branch sometimes known as W_0.  (See LambertW1 for calculating the other // branch) // // The Lambert W function is what you get when you try to solve //   z = w e^w   for w. // // You really want to be using this with Frink:The Next Generation // // See:  http://mathworld.wolfram.com/LambertW-Function.html // and   http://en.wikipedia.org/wiki/Lambert_W_function //         for a general solution. // // This algorithm was originally based on the algorithm for calculating // LambertW using Newton's method, as outlined at: // http://www.whim.org/nebula/math/lambertw.html // but it actually now uses Halley's method because it converges // faster. // // This also uses as its initial guess very good approximation to Lambert W // with a relative error of less than 0.01 in the entire complex plane. // From the paper: // "Uniform approximations for transcendental functions", Serge Winitzki // // For this, see "Numerical Evaluation of the Lambert W Function // and Application to Generation of Generalized Gaussian Noise // With Exponent 1/2" by Francois Chapeau-Blondeau and Abdelilah Monir, // IEEE Transactions on Signal Processing, Vol. 50, No. 9, // September 2002 // https://pdfs.semanticscholar.org/7a5a/76a9369586dd0dd34dda156d8f2779d1fd59.pdf // LambertW[x, digits=getPrecision[]] := { //   println["In LambertW, x = \$x, digits = \$digits"]    oldPrec = getPrecision[]    err = 10^-(digits+2) //   println["Error is \$err"]    if x == 0       return 0    // TODO:  Get cached value from E class using binary splitting.    e = arbitraryExp[1, digits+5]    // println["e is \$e"]    // Get a good first guess.    w = estimateLambertW[x]    // println["Estimate is \$w"]        setPrecision[digits+4]    retval = 1          iterCount = 0           try    {       do       {          oldw = w          iterCount = iterCount + 1          if (w == -1)             return -1          else             if (w < -.35) and (iterCount mod 3 == 2)             {                // Newton's method sometimes bounces around                // -1/e (approx -0.367879) so we sometimes                // use this alternate method to try and make it converge.                w = -1 + (w+1) sqrt[(x + 1/e) / (w arbitraryExp[w] + 1/e),                                    digits+4]                if w < -1                   w = -1             } else             {                // Normal case, Newton's method                // We normally use Halley's method (below) which is more                // complicated but has a higher order of convergence, converging                // in about half the number of steps in the usual case.  We use                // Newton's method if we haven't converged after a long time                // just in case something's going weird in Halley's method.                if (iterCount > 20) and (iterCount mod 2 == 0)                   w = (x arbitraryExp[-w] + w^2) / (w+1)                else                {                   // Halley's method, see Eq. 9 in Chapeau-Blondeau                   ew = arbitraryExp[w]                   wew = w ew                   wewmz = wew - x                   w = w - wewmz / ((w + 1) ew - ((w+2) wewmz)/(2 w + 2))                }             }          // println[w]       } while (abs[w - oldw] > abs[w err])       setPrecision[digits]       retval = 1.0 * w    }    finally    {       setPrecision[oldPrec]    } //   println["In LambertW, returning \$retval"]    return retval } /** Function to test the difference between the calculated value of LambertW     and its inverse function. */ LambertWTest[x, digits=getPrecision[]] := {    origPrec = getPrecision[]    try    {       setPrecision[digits]       w = LambertW[x, digits]       reverse = w arbitraryExp[w, digits] //      println["\$x\t\$w"]       return x - reverse    }    finally       setPrecision[origPrec] } // This procedure gives a very good approximation to Lambert W // with a relative error of less than 0.01 in the entire complex plane. // From the paper: // "Uniform approximations for transcendental functions", Serge Winitzki estimateLambertW[x] := {    if (x < 1)    {       ex = e x       return ex / (1 + ((2ex + 2)^(-1/2) + 1.2890834)^-1) //eq39    } else    {       ln1x = ln[1+x]       return ln1x (1 - (ln[1 + ln1x])/(2 + ln1x))   // Eq. 38    } } /** This calculates the secondary branch of the Lambert W function known as     W_-1.  It is defined over the domain [-1/e, 0].  */ LambertW1[x, digits=getPrecision[]] := {    oldPrec = getPrecision[]    err = 10^-(digits+2) //   println["Error is \$err"]    if x == 0       return 0    // Get a good first guess.    w = estimateLambertW1[x]    // println[w]    if w == undef or w >= -1       w = -1. - 10^-(min[oldPrec,digits]/2)    // println[w]        setPrecision[digits+4]    retval = 1          iterCount = 0           try    {       do       {          oldw = w          iterCount = iterCount + 1          // Normal case, Newton's method          // We normally use Halley's method (below) which is more          // complicated but has a higher order of convergence, converging          // in about half the number of steps in the usual case.  We use          // Newton's method if we haven't converged after a long time          // just in case something's going weird in Halley's method.          if (iterCount > 10) and (iterCount mod 3 == 0)             w = (x arbitraryExp[-w] + w^2) / (w+1)          else          {             // Halley's method, see Eq. 9 in Chapeau-Blondeau             ew = arbitraryExp[w]             wew = w ew             wewmz = wew - x             w = w - wewmz / ((w + 1) ew - ((w+2) wewmz)/(2 w + 2))             // println[w]          }       } while (abs[w - oldw] > abs[w err])       setPrecision[digits]       retval = 1.0 * w    }    finally    {       setPrecision[oldPrec]    }    return retval } /** This function gives a good initial guess for the Lambert W-1 function.     The relative error is less than 10^-4 for any value in its domain [-1/e, 0]     This should be used to seed a more accurate convergence procedure.     This algorithm is Fig. 2 in "Numerical Evaluation of the Lambert W Function     and Application to Generation of Generalized Gaussian Noise     With Exponent 1/2" by Francois Chapeau-Blondeau and Abdelilah Monir,     IEEE Transactions on Signal Processing, Vol. 50, No. 9,     September 2002      https://pdfs.semanticscholar.org/7a5a/76a9369586dd0dd34dda156d8f2779d1fd59.pd */ estimateLambertW1[z] := {    if -1/e <= z and z < -0.333    {       p = -sqrt[2 (e z + 1)]       return -1 + p - 1/3 p^2 + 11/72 p^3 - 43/540 p^4 + 769/17280 p^5 - 221/8505 p^6    } else       if -0.333 <= z and z <= -0.033       {          return (-8.0960 + 391.0025 z - 47.4252 z^2 - 4877.6330 z^3 - 5532.7760 z^4) / (1 - 82.9423 z + 433.8688 z^2 + 1515.3060 z^3)       } else          if -0.033 < z and z < 0          {             L1 = ln[-z]             L2 = ln[-ln[-z]]             return L1 - L2 + L2/L1 + ((-2+L2) L2)/(2 L1^2) +                ((6 - 9 L2 + 2 L2^2)L2)/(6 L1^3) +                ((-12 + 36 L2 - 22 L2^2 + 3 L2^3) L2)/ (12 L1^4) +                ((60 - 300 L2 + 350 L2^2 - 125 L2^3 + 12 L2^4) L2)/ (60 L1^5)          }     return undef          } /* digits = 500 for z = -1/e + 1e-15 to 20 step 0.1 {    err = LambertWTest[z, digits]    if abs[err] > 10^-(digits-1)       println["\$z\t\$err"] } */```

This is a program written in the programming language Frink.
For more information, view the Frink Documentation or see More Sample Frink Programs.

Alan Eliasen was born 18593 days, 12 hours, 37 minutes ago.