View or download LambertW.frink in plain text format
use ArbitraryPrecision.frink
use sqrtWayne.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[]] :=
{
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]
// Get a good first guess.
w = estimateLambertW[x]
// println[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]
}
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.pd
*/
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"]
}
*/
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 18234 days, 8 hours, 6 minutes ago.