LambertW.frink

View or download LambertW.frink in plain text format


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"]
}
*/


View or download LambertW.frink in plain text format


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 18376 days, 23 hours, 33 minutes ago.