Download or view LambertW.frink in plain text format
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"]
}
*/
Download or view 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 19966 days, 13 hours, 59 minutes ago.