# 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 18491 days, 4 hours, 26 minutes ago.