LambertW.frink

View or download LambertW.frink in plain text format


use ArbitraryPrecision.frink
use root.frink

// This calculates the value of the Lambert W function.
//
// See:  http://mathworld.wolfram.com/LambertW-Function.html
// and   http://en.wikipedia.org/wiki/Lambert_W_function
//         for a general solution.
//
// This algorithm is based on the algorithm for calculating
// LambertW using Newton's method, as outlined at:
// http://www.whim.org/nebula/math/lambertw.html
//
// TODO:  Calculate non-principal value.
// TODO:  Improve convergence rate around -1/e
LambertW[x, digits=getPrecision[]] :=
{
   oldPrec = getPrecision[]
   err = 10^-(digits+1)
   println["Error is $err"]
   if x == 0
      return 0

   if -1/e <= x and x <= 10
      w = 0
   else
      if x > 10
         w = ln[x] - ln[ln[x]]
      else
         return undef

   setPrecision[digits+2]
   retval = 1      

   try
   {
      do
      {
         oldw = w
         if (digits > 15)
            exp = arbitraryExp[-w]
         else
            exp = e^-w
         
         if (w == -1)
            return -1
         else
            if (w < -.35)  // Slow to converge around -1/e (approx -0.367879)
               w = -1 + (w+1) sqrt[(x + 1/e) / (w arbitraryExp[w] + 1/e)]
            else
               w = (x exp + w^2) / (w+1)  // Normal case
         println[w]
      } while (abs[w - oldw] > abs[w err])

      setPrecision[digits]
      retval = 1.0 * w
   }
   finally
   {
      setPrecision[oldPrec]
   }

   return retval
}


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 17778 days, 20 hours, 28 minutes ago.