use ArbitraryPrecision.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.pdf
*/
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"]
}
*/