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.