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 18010 days, 7 hours, 50 minutes ago.