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[]] := { // 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.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"] } */