implement Lambert W function by jlapeyre · Pull Request #84 · JuliaMath/SpecialFunctions.jl

jlapeyre · 2018-03-28T18:33:52Z

This copies the implementation of the Lambert W function and
related functions and data from jlapeyre/LambertW.jl to
SpecialFunctions.jl

This copies the implementation of the Lambert W function and related functions and data from jlapeyre/LambertW.jl to SpecialFunctions.jl

jlapeyre · 2018-03-28T18:40:27Z

Including the Lambert W function was first mentioned in #27

jlapeyre · 2018-04-18T10:27:01Z

Anyone home ?

ararslan

I've left some comments but I can't speak to the accuracy of the algorithms.

ararslan · 2018-04-18T15:12:09Z

src/lambertw.jl

+#export lambertw, lambertwbp
+using Compat
+
+const euler =


You can just do

using Compat.MathConstants: e

then use e instead of euler.

Thanks. I couldn't figure out the compat syntax from the docs, so I just did what I could.

ararslan · 2018-04-18T15:13:30Z

src/lambertw.jl

+
+function __init__()
+    omega_const_bf_[] =
+        parse(BigFloat,"0.5671432904097838729999686622103555497538157871865125081351310792230457930866845666932194")


Why is it necessary to rerun this every time the module is loaded? Seems like this should just be

const omega_const_bf = parse(BigFloat, "...")

IIUC, some storage is malloced when the big float is constructed and Julia keeps the pointer. The next time the module runs, the memory has been freed (of course), but the address Julia stored has not changed. BANG!... If you don't precompile, then the heap allocation is done every time the module is loaded. But, I want to precompile. AFAICT, this is the recommended way to a) allocate the storage every time the module is used. b) precompile the module.

ararslan · 2018-04-18T15:15:25Z

src/lambertw.jl

+    two_t = convert(T,2)
+    lastx = x
+    lastdiff = zero(T)
+    for i in 1:100


Where does 100 come from? Seems like a rather arbitrary "magic number"; perhaps we could have that be configurable, e.g.

function _lambertw(z::T, x::T, maxiter::Int=100) where T<:Number

Two things. 1) Some of this code was taken from mpmath. Maybe this number was, although I don't think a reason was given there either. I should acknowledge the source somewhere, it has a MIT compatible license. I will look at the python source. 2) Yes, I think changing it is ok. How ? Maybe.... pick a number far larger than the number of iterations required to converge in most, or all cases. But, make it small enough to finish the loop in a short time. We could have it error or warn about lack of precision if the entire loop is completed. I'm pretty sure I had something like this in there at one point.

ararslan · 2018-04-18T15:15:52Z

src/lambertw.jl

+        ex = exp(x)
+        xexz = x * ex - z
+        x1 = x + 1
+        x = x - xexz / (ex * x1 - (x + two_t) * xexz / (two_t * x1 ) )


x -= instead of x = x -

ararslan · 2018-04-18T15:18:16Z

src/lambertw.jl

+
+# The fancy initial condition selection does not seem to help speed, but we leave it for now.
+function lambertwk0(x::T)::T where T<:AbstractFloat
+    x == Inf && return Inf


Perhaps isfinite(x) || return x? That will also handle the case where x is NaN, or if T is anything other than Float64 (Inf is a Float64 literal).

I don't recall my reasoning here. It is highly likely that I compared to python and Mathematica. But, maybe there is a Julia convention. Some clues for what to do:

julia> sin(NaN) NaN julia> sin(Inf) ERROR: DomainError with Inf: sin(x) is only defined for finite x.

But sin does not tend to a limit with increasing x. So Inf is in no sense the correct result. lambertw tends to infinity, so Inf in this sense is the correct answer.

In [5]: lambertw(inf) # python mpmath Out[5]: mpf('+inf') In [9]: lambertw(nan) Out[9]: mpf('nan')

Another issue is:

julia> convert(BigFloat,Inf) Inf julia> typeof(convert(BigFloat,Inf)) BigFloat

Yikes, no visible difference. In any case, at a minimum, this should be convert(T,Inf).

ararslan · 2018-04-18T15:20:14Z

src/lambertw.jl

+
+# maybe compute higher precision. converges very quickly
+function omega_const(::Type{BigFloat})
+  @compat  precision(BigFloat) <= 256 && return omega_const_bf_[]


Weird indentation here, and why is @compat necessary?

I think compatibility with v0.4. There is already a deprecation in v0.5:

julia> get_bigfloat_precision() WARNING: get_bigfloat_precision() is deprecated

In other words, this can be removed. (Although I am interested what you/community think is the proper indentation, were @compat necessary here.

The Julia standard for indentation is 4 spaces. That's what's used in the Julia repository as well as the majority of packages.

We don't need compatibility with 0.4 or even 0.5 here, since the package requires 0.6 at a minimum.

ararslan · 2018-04-18T15:21:29Z

src/lambertw.jl

+
+const LAMWMU_FLOAT64 = lamwcoeff(Float64,500)
+
+function horner(x, p::AbstractArray,n)


Why not use Base.Math.@horner?

I knew about Base.Math.@horner when I wrote this. So, the answer is I don't know. I'll look into it.

ararslan · 2018-04-18T15:22:30Z

src/lambertw.jl

+function lambertwbp(x::Number,k::Integer)
+    k == 0 && return _lambertw0(x)
+    k == -1 && return _lambertwm1(x)
+    error("expansion about branch point only implemented for k = 0 and -1")


An error regarding the arguments should be throw(ArgumentError("...")) instead of error("...").

Right. will do.

ararslan · 2018-04-18T15:23:51Z

src/lambertw.jl

+
+#### Lambert W function ####
+
+const LAMBERTW_USE_NAN = false


I must say, I don't really like this. I think we should pick a convention that matches the other functions here and stick to it.

Ok. IIRC, there was some debate around the time I wrote the module about how to handle this in general. I also wanted to see which I liked better. I think the issue is settled, so sure, I'll try to find the convention. It can only be changed by editing the source. In SpecialFunctions this certainly shouldn't be advertised as something to do.

jlapeyre · 2018-04-18T17:00:42Z

@ararslan thanks for the review. I'll do these things.

jlapeyre · 2018-04-18T18:51:43Z

How do you prefer updates to the PR ? I could rebase and force push, but I think we lose link between comments and code this way.

@stevengj I suggest putting off reviewing until this round of changes are made to avoid duplicating effort / wasting time.

ararslan · 2018-04-18T19:01:35Z

You can just push more commits to your branch and they'll show up here. We can squash on merge.

jlapeyre · 2018-04-19T18:48:55Z

The latest commit addresses all comments made by @ararslan. I did a bit of other refactoring and tidying.

jlapeyre · 2018-07-06T10:25:52Z

I'd like to close this and issue a new PR after v1.0. In the upstream source LambertW.jl, I dropped v0.6 support, and have improved clarity and compliance with style standards. The upstream source already incorporates changes made during the code review here. (I have not yet pushed all these changes to github.)

stevengj · 2018-07-12T02:07:52Z

Just do a forced push to this branch to update the existing PR.

cossio · 2018-09-26T16:08:11Z

I hope this gets merged soon. :)

cossio · 2018-09-27T22:00:05Z

src/SpecialFunctions.jl

+function __init__()
+    # allocate storage for this BigFloat constant each time this module is loaded
+    omega_const_bf_[] =
+        parse(BigFloat,"0.5671432904097838729999686622103555497538157871865125081351310792230457930866845666932194")


Isn't big"0.56..." better than parse(...)?

If you check, for instance with @less big"1.0", you see that the macro does more work.

You're right. @code_llvm big"0.56..." gives a LOT more code.

Also, big"..." doesn't respect current set precision:

julia> setprecision(1024) do precision(big"1.0") end 256

Is that intended behavior of big"..."? Or should one raise an issue?

giordano · 2019-05-20T16:45:59Z

Bump 🙂 I'd also love to have this one merged

chriscoey · 2020-01-05T13:53:45Z

bump. need this one quite often!

longemen3000 · 2020-02-05T22:32:40Z

bumping

alyst

Thanks for this PR, it would be a very nice addition to SpecialFunctions.jl!
Also merging it here would ensure that the implementation would be maintained by the community.

alyst · 2021-11-15T14:27:35Z

src/lambertw.jl

+    two_t = convert(T,2)
+    lastx = x
+    lastdiff = zero(T)
+    converged::Bool = false


Suggested change

converged::Bool = false

converged = false

I guess type-assert is not required here

alyst · 2021-11-15T14:28:46Z

src/lambertw.jl

+        ex = exp(x)
+        xexz = x * ex - z
+        x1 = x + 1
+        x -= xexz / (ex * x1 - (x + two_t) * xexz / (two_t * x1 ) )


Suggested change

x -= xexz / (ex * x1 - (x + two_t) * xexz / (two_t * x1 ) )

x -= xexz / (ex * x1 - (x + two_t) * xexz / (two_t * x1))

Also, does it really help to define two_t = convert(T, 2)? Maybe the compiler can figure out the best way to calculate x + 2 (like it's done for x + 1 above)?

Suggested change

x -= xexz / (ex * x1 - (x + two_t) * xexz / (two_t * x1 ) )

x -= xexz / (ex * x1 - (x + 2) * xexz / (2 * x1))

alyst · 2021-11-15T14:36:08Z

src/lambertw.jl

+        x1 = x + 1
+        x -= xexz / (ex * x1 - (x + two_t) * xexz / (two_t * x1 ) )
+        xdiff = abs(lastx - x)
+        if xdiff <= 3*eps(abs(lastx)) || lastdiff == xdiff  # second condition catches two-value cycle


Suggested change

if xdiff <= 3*eps(abs(lastx)) || lastdiff == xdiff # second condition catches two-value cycle

if xdiff <= 3*eps(lastx) || lastdiff == xdiff # second condition catches two-value cycle

alyst · 2021-11-15T14:37:36Z

src/lambertw.jl

+        lastx = x
+        lastdiff = xdiff
+    end
+    converged || warn("lambertw with z=", z, " did not converge in ", maxits, " iterations.")


Suggested change

converged || warn("lambertw with z=", z, " did not converge in ", maxits, " iterations.")

converged || warn("lambertw(", z, ") did not converge in ", maxits, " iterations.")

alyst · 2021-11-15T14:44:07Z

src/lambertw.jl

+
+# Use Halley's root-finding method to find
+# x = lambertw(z) with initial point x.
+function _lambertw(z::T, x::T, maxits) where T <: Number


Suggested change

function _lambertw(z::T, x::T, maxits) where T <: Number

function _lambertw(z::T, x::T, maxits::Integer) where T <: Number

alyst · 2021-11-15T15:15:16Z

src/lambertw.jl

+# solve for α₂. We get α₂ = 0.
+# compute array of coefficients μ in (4.22).
+# m[1] is μ₀
+function lamwcoeff(T::DataType, n::Int)


Suggested change

function lamwcoeff(T::DataType, n::Int)

function lambertw_coefficients(T::DataType, n::Int)

alyst · 2021-11-15T15:19:34Z

src/lambertw.jl

+    return m
+end
+
+const LAMWMU_FLOAT64 = lamwcoeff(Float64,500)


Suggested change

const LAMWMU_FLOAT64 = lamwcoeff(Float64,500)

const LambertW_μ = lambertw_coefficients(Float64, 500)

alyst · 2021-11-15T15:41:40Z

src/lambertw.jl

+    x < 4e-11 && return wser3(p)
+    x < 1e-5 && return wser7(p)
+    x < 1e-3 && return wser12(p)
+    x < 1e-2 && return wser19(p)
+    x < 3e-2 && return wser26(p)
+    x < 5e-2 && return wser32(p)
+    x < 1e-1 && return wser50(p)
+    x < 1.9e-1 && return wser100(p)
+    x > 1/MathConstants.e && throw(DomainError(x))  # radius of convergence
+    return wser290(p)  # good for x approx .32


Suggested change

x < 4e-11 && return wser3(p)

x < 1e-5 && return wser7(p)

x < 1e-3 && return wser12(p)

x < 1e-2 && return wser19(p)

x < 3e-2 && return wser26(p)

x < 5e-2 && return wser32(p)

x < 1e-1 && return wser50(p)

x < 1.9e-1 && return wser100(p)

x > 1/MathConstants.e && throw(DomainError(x)) # radius of convergence

return wser290(p) # good for x approx .32

x < 4e-11 && return Base.Math.@horner p LambertW_μ[1:3]...

x < 1e-5 && return Base.Math.@horner p LambertW_μ[1:7]...

n = x < 1e-3 ? 12 :

x < 1e-2 ? 19 :

x < 3e-2 ? 26 :

x < 5e-2 ? 32 :

x < 1e-1 ? 50 :

x < 1.9e-1 ? 100 :

x <= 1/MathConstants.e ? 290 : # good for x approx .32

throw(DomainError(x)) # outside of convergence range

return evalpoly(p, view(LambertW_μ, 1:n))

alyst · 2021-11-15T15:42:19Z

src/lambertw.jl

+
+# Converges to Float64 precision
+# We could get finer tuning by separating k=0,-1 branches.
+function wser(p,x)


Suggested change

function wser(p,x)

function lambertw_series(p, x)

alyst · 2021-11-15T15:52:01Z

src/lambertw.jl

+# Base.Math.@horner requires literal coefficients
+# But, we have an array `p` of computed coefficients
+function horner(x, p::AbstractArray, n)
+    n += 1
+    ex = p[n]
+    for i = n-1:-1:2
+        ex = :(muladd(t, $ex, $(p[i])))
+    end
+    ex = :( t * $ex)
+    return Expr(:block, :(t = $x), ex)
+end
+
+function mkwser(name, n)
+    iex = horner(:x,LAMWMU_FLOAT64,n)
+    return :(function ($name)(x) $iex  end)
+end
+
+eval(mkwser(:wser3, 3))
+eval(mkwser(:wser5, 5))
+eval(mkwser(:wser7, 7))
+eval(mkwser(:wser12, 12))
+eval(mkwser(:wser19, 19))
+eval(mkwser(:wser26, 26))
+eval(mkwser(:wser32, 32))
+eval(mkwser(:wser50, 50))
+eval(mkwser(:wser100, 100))
+eval(mkwser(:wser290, 290))


I think it's possible to use @horner/evalpoly() directly, so these functions are not needed.

Suggested change

# Base.Math.@horner requires literal coefficients

# But, we have an array `p` of computed coefficients

function horner(x, p::AbstractArray, n)

n += 1

ex = p[n]

for i = n-1:-1:2

ex = :(muladd(t, $ex, $(p[i])))

end

ex = :( t * $ex)

return Expr(:block, :(t = $x), ex)

end

function mkwser(name, n)

iex = horner(:x,LAMWMU_FLOAT64,n)

return :(function ($name)(x) $iex end)

end

eval(mkwser(:wser3, 3))

eval(mkwser(:wser5, 5))

eval(mkwser(:wser7, 7))

eval(mkwser(:wser12, 12))

eval(mkwser(:wser19, 19))

eval(mkwser(:wser26, 26))

eval(mkwser(:wser32, 32))

eval(mkwser(:wser50, 50))

eval(mkwser(:wser100, 100))

eval(mkwser(:wser290, 290))

alyst · 2021-11-15T18:07:36Z

src/lambertw.jl

+### omega constant ###
+
+const omega_const_ = 0.567143290409783872999968662210355
+# The BigFloat `omega_const_bf_` is set via a literal in the function __init__ to prevent a segfault


I wonder if the issue (segfault) that required omega_const_bf initialization in __init__() is resolved, and it could be declared here:

Suggested change

### omega constant ###

const omega_const_ = 0.567143290409783872999968662210355

# The BigFloat `omega_const_bf_` is set via a literal in the function __init__ to prevent a segfault

### Omega constant: Ω exp(Ω) = 1 (Ω = lambertw(1)) ###

const Omega_constant = 0.567143290409783872999968662210355

const Omega_constant_BigFloat = BigFloat("0.5671432904097838729999686622103555497538157871865125081351310792230457930866845666932194")

Or maybe Omega constant could be moved to IrrationalConstants.jl?

The "workaround" for the segfault was the standard way to deal with this problem. AFAIK, it has not changed. This is work that has to be done after the dynamic library (GMP, or something like that in this case) is linked.

I see. The potential alternative that does not implicate __init__() could be

const Omega_BigFloat256 = Ref{BigFloat}() @irrational Omega 0.567143290409783872999968662210355 Omega_constant() # compute BigFloat Omega constant at arbitrary precision function Omega_constant() isassigned(Omega_BigFloat256) || (Omega_BigFloat256[] = BigFloat("0.5671432904097838729999686622103555497538157871865125081351310792230457930866845666932194")) oc = Omega_BigFloat256[] precision(BigFloat) <= 256 && return oc myeps = eps(BigFloat) for i in 1:100 nextoc = (1 + oc) / (1 + exp(oc)) abs(oc - nextoc) <= myeps && break oc = nextoc end return oc end

jlapeyre · 2021-11-15T18:57:40Z

@alyst, this PR is out of date. LambertW has changed in the meantime. I (or we) need to merge the changes into this PR branch.

alyst · 2021-11-15T19:25:53Z

@alyst, this PR is out of date. LambertW has changed in the meantime. I (or we) need to merge the changes into this PR branch.

Let me know how you want to do it.

ViralBShah · 2021-12-29T08:27:16Z

Seems like lots of people want this, but maybe it is easiest to just refresh the package that @jlapeyre has. Is it registered? Should we move it to JuliaMath?

jlapeyre/LambertW.jl

jlapeyre · 2021-12-31T15:41:55Z

It's not clear to me the advantage of moving lambertw to SpecialFunctions. A few people want it moved, but I didn't see an argument. Maybe they use it frequently and already import SpecialFunctions a lot, and don't want to be required to do import LambertW.

@ViralBShah : Yes, this package is registered. I very much support moving this to JuliaMath.

It seems like a bit of work to get this moved to SpecialFunctions. One issue is the big, flat, namespace in SpecialFunctions. What to do with LambertW.omega? I think SpecialFunctions.omega is not a good choice.

I'm not completely against moving LambertW.jl, but I'd like to see arguments for why.

alyst · 2021-12-31T16:34:48Z

I'm not completely against moving LambertW.jl, but I'd like to see arguments for why.

it's more likely that packages like SpecialFunctions.jl will have (co-)maintainers in the long run
even if (co-)maintainers of "fine-grained" packages are available, it would be an extra burden for them to do the tedious, but necessary work of keeping the API, documentation, unit tests and CI scripts up-to-date
Labert W function does not depend on any other packages, and the code base is sufficiently lean and straightforward, so adding it to SpecialFunctions.jl would not make the latter heavier
it's a matter of taste, but I find it's easier to discover/evaluate methods/tools when they are organized into dedicated registries like SpecialFunctions.jl

cossio · 2021-12-31T17:25:56Z

What to do with LambertW.omega? I think SpecialFunctions.omega is not a good choice.

How about using lambertw(1) as the public API to access the value of the Omega constant? No need for an additional name. Just guarantee that lambertw(1) has a fast evaluation path.

I think SpecialFunctions.omega is not a good choice.

Another option is SpecialFunctions.lambert_omega, and it doesn't have to be exported.

simonbyrne · 2022-01-02T00:26:34Z

it's more likely that packages like SpecialFunctions.jl will have (co-)maintainers in the long run

My experience so far has been the opposite: large packages like SpecialFunctions.jl and Distributions.jl struggle to keep maintainers. Smaller packages tend to have an easier time keeping maintainers (I presume due to a lower maintenance burden, and the code author feeling an ownership of the package).

ViralBShah · 2022-01-02T09:41:15Z

However, lots of small packages are a burden to the user.

originalsouth · 2022-01-02T09:45:46Z

Perhaps, as a middle ground solution, the functions can be reexported then?

devmotion · 2022-01-02T15:32:28Z

I agree with @simonbyrne, I don't think it's necessarily better maintained if it is part of SpecialFunctions.

IMO separate packages are fine in this case. In my experience they are annoying mainly if one always has to load a specific set of packages to achieve some basic functionality but I don't think this is the case here. One can use LambertW just fine without loading SpecialFunctions, and similarly SpecialFunctions without LambertW.

One problem with reexporting LambertW is that any breaking release of LambertW requires a breaking release of SpecialFunctions. Due to the large number of dependencies (1625 packages), it takes quite a while until a breaking release of SpecialFunctions is propagated through the Julia ecosystem and it requires work and time from many downstream developers, whereas a breaking release of LambertW (if SpecialFunctions does not reexport it) would have less severe effects since it only has a few (4 to be precise) dependencies.

alyst · 2022-01-04T11:01:12Z

How about using lambertw(1) as the public API to access the value of the Omega constant? No need for an additional name. Just guarantee that lambertw(1) has a fast evaluation path.

The problem is that Omega should be of type Base.Irrational{:XXX} (as the other constants), which could then be explicitly converted to a numeric type the user specifies, so it cannot be returned by lambertw(1) directly for the sake of type stability.
But the fast path to calculate lambertw(1)/lambertw(big(1))/lambertw(complex(1,0)) is already there.

Smaller packages tend to have an easier time keeping maintainers (I presume due to a lower maintenance burden, and the code author feeling an ownership of the package).

"Feeling an ownership" is a double-edged sword. In my experience the code quality of community-maintained packages is better -- more eyes to look at the code, more chances to keep up with the evolving Julia standards, more predictable PR review and acceptance policy.
Anyway, what I am saying about community maintenance only concerns the functionality that is "adjacent" to the core of the language. Clearly, SpecialFunctions.jl belongs to this category. I believe lambertw() belongs there too, given its numerous applications in different branches of science. But then again, it is up to the core Julia "stakeholders" to define and implement the strategy that would keep such functionality in good shape.
For my part I have just implemented in #371 what I have thought is need for this PR to be accepted.

cossio · 2022-01-04T11:09:28Z

In my experience the code quality of community-maintained packages is better -- more eyes to look at the code, more chances to keep up with the evolving Julia standards, more predictable PR review and acceptance policy.

It seems that just transferring to LambertW.jl to JuliaMath would also achieve this? It could then be seen as a community package, no longer with a sole owner. Also more people will have commit access, code reviews, etc. ...
And that way it remains independent of SpecialFunctions.jl, which has a lot more ~~dependencies~~dependants.

Looked at it this way, it would seem that transferring LambertW.jl to JuliaMath is the happy solution for everybody.

alyst · 2022-01-04T11:17:19Z

Looked at it this way, it would seem that transferring LambertW.jl to JuliaMath is the happy solution for everybody.

Authorship aside, it's yet another package to make sure that e.g. .github/workflows or docs/make.jl scripts are up-to-date.
So far lambertw() doesn't have any dependencies, so I don't expect that in the future SpecialFunctions.jl would need minor version updates because of lambertw()
(Well, for my fork I had to add Compat to support evalpoly() on Julia 1.3, but if that is critical, I can switch to @evalpoly macro that should work everywhere).

jlapeyre · 2022-01-09T20:19:08Z

I thumbs-upped the comments supporting not putting LambertW in SpecialFunctions. I propose moving it to JuliaMath. This does not preclude, or make it any more difficult, to later move it to SpecialFunctions. I know there are arguments both ways. The problem is you can't put yourself in the shoes of every user and know what their workflow is like.

cossio · 2023-03-07T22:03:32Z

@simonbyrne Can I ask why was this closed? Is it still the plan to transfer lambertw to this package?

ararslan · 2023-03-08T00:27:19Z

I assume it was closed in favor of #371.

simonbyrne · 2023-03-08T19:05:40Z

That, and it seems there was a consensus in keeping it as a separate package?

ararslan · 2023-03-08T19:38:20Z

I didn't interpret there to be a strong consensus one way or another

implement Lambert W function

5c16731

This copies the implementation of the Lambert W function and related functions and data from jlapeyre/LambertW.jl to SpecialFunctions.jl

jlapeyre mentioned this pull request Apr 9, 2018

Add (and unify) Julia Lambert W function? #27

Open

ararslan reviewed Apr 18, 2018

View reviewed changes

ararslan requested a review from stevengj April 18, 2018 15:24

made changes request in PR review

9ef0205

cossio reviewed Sep 27, 2018

View reviewed changes

This was referenced Sep 4, 2021

Move LambertW to SpecialFunctions JuliaMath/LambertW.jl#8

Open

Add ZeroInflatedPoisson distribution JuliaStats/Distributions.jl#1393

Open

alyst reviewed Nov 15, 2021

View reviewed changes

alyst mentioned this pull request Nov 19, 2021

Lambert's Omega constant JuliaMath/IrrationalConstants.jl#12

Closed

alyst mentioned this pull request Dec 31, 2021

Updated Lambert W function implementation #371

Open

simonbyrne closed this Mar 7, 2023


		const LAMWMU_FLOAT64 = lamwcoeff(Float64,500)

		function horner(x, p::AbstractArray,n)

	x -= xexz / (ex * x1 - (x + two_t) * xexz / (two_t * x1 ) )
	x -= xexz / (ex * x1 - (x + two_t) * xexz / (two_t * x1))

	x -= xexz / (ex * x1 - (x + two_t) * xexz / (two_t * x1 ) )
	x -= xexz / (ex * x1 - (x + 2) * xexz / (2 * x1))

	if xdiff <= 3*eps(abs(lastx)) \|\| lastdiff == xdiff # second condition catches two-value cycle
	if xdiff <= 3*eps(lastx) \|\| lastdiff == xdiff # second condition catches two-value cycle

	converged \|\| warn("lambertw with z=", z, " did not converge in ", maxits, " iterations.")
	converged \|\| warn("lambertw(", z, ") did not converge in ", maxits, " iterations.")

	function _lambertw(z::T, x::T, maxits) where T <: Number
	function _lambertw(z::T, x::T, maxits::Integer) where T <: Number

	function lamwcoeff(T::DataType, n::Int)
	function lambertw_coefficients(T::DataType, n::Int)

	const LAMWMU_FLOAT64 = lamwcoeff(Float64,500)
	const LambertW_μ = lambertw_coefficients(Float64, 500)

-    x < 4e-11 && return wser3(p)
-    x < 1e-5 && return wser7(p)
-    x < 1e-3 && return wser12(p)
-    x < 1e-2 && return wser19(p)
-    x < 3e-2 && return wser26(p)
-    x < 5e-2 && return wser32(p)
-    x < 1e-1 && return wser50(p)
-    x < 1.9e-1 && return wser100(p)
-    x > 1/MathConstants.e && throw(DomainError(x))  # radius of convergence
-    return wser290(p)  # good for x approx .32
+    x < 4e-11 && return Base.Math.@horner p LambertW_μ[1:3]...
+    x < 1e-5 && return Base.Math.@horner p LambertW_μ[1:7]...
+    n = x < 1e-3 ? 12 :
+        x < 1e-2 ? 19 :
+        x < 3e-2 ? 26 :
+        x < 5e-2 ? 32 :
+        x < 1e-1 ? 50 :
+        x < 1.9e-1 ? 100 :
+        x <= 1/MathConstants.e ? 290 : # good for x approx .32
+        throw(DomainError(x)) # outside of convergence range
+    return evalpoly(p, view(LambertW_μ, 1:n))

Uh oh!

Conversation

jlapeyre commented Mar 28, 2018

Uh oh!

jlapeyre commented Mar 28, 2018

Uh oh!

jlapeyre commented Apr 18, 2018

Uh oh!

ararslan left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

jlapeyre Apr 18, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

jlapeyre commented Apr 18, 2018

Uh oh!

jlapeyre commented Apr 18, 2018

Uh oh!

ararslan commented Apr 18, 2018

Uh oh!

jlapeyre commented Apr 19, 2018

Uh oh!

jlapeyre commented Jul 6, 2018

Uh oh!

stevengj commented Jul 12, 2018

Uh oh!

cossio commented Sep 26, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

giordano May 20, 2019 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

giordano commented May 20, 2019

Uh oh!

jlapeyre Apr 18, 2018 •

edited

Loading

cossio commented Sep 26, 2018 •

edited

Loading

giordano May 20, 2019 •

edited

Loading