diff --git a/LICENSE b/LICENSE index 28e07032..b9a9755c 100644 --- a/LICENSE +++ b/LICENSE @@ -22,3 +22,35 @@ LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. +Portions of this code are derived from SciPy and are licensed under +the Scipy License: + +> Copyright (c) 2001, 2002 Enthought, Inc. +> All rights reserved. + +> Copyright (c) 2003-2012 SciPy Developers. +> All rights reserved. + +> Redistribution and use in source and binary forms, with or without +> modification, are permitted provided that the following conditions are met: + +> a. Redistributions of source code must retain the above copyright notice, +> this list of conditions and the following disclaimer. +> b. Redistributions in binary form must reproduce the above copyright +> notice, this list of conditions and the following disclaimer in the +> documentation and/or other materials provided with the distribution. +> c. Neither the name of Enthought nor the names of the SciPy Developers +> may be used to endorse or promote products derived from this software +> without specific prior written permission. +> +> THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +> AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +> IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +> ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS +> BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, +> OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +> SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +> INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +> CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +> ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF +> THE POSSIBILITY OF SUCH DAMAGE. diff --git a/Project.toml b/Project.toml index cf84beb2..a0bde5eb 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ version = "2.2.0" [deps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" +Compat = "34da2185-b29b-5c13-b0c7-acf172513d20" IrrationalConstants = "92d709cd-6900-40b7-9082-c6be49f344b6" LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688" OpenLibm_jll = "05823500-19ac-5b8b-9628-191a04bc5112" @@ -18,6 +19,7 @@ SpecialFunctionsChainRulesCoreExt = "ChainRulesCore" [compat] ChainRulesCore = "0.9.44, 0.10, 1" ChainRulesTestUtils = "0.6.8, 0.7, 1" +Compat = "3.7, 4" IrrationalConstants = "0.1, 0.2" LogExpFunctions = "0.3.2" OpenLibm_jll = "0.7, 0.8" diff --git a/README.md b/README.md index 206679b8..d4e5ef42 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,7 @@ # SpecialFunctions.jl Special mathematical functions in Julia, including Bessel, Hankel, Airy, error, Dawson, exponential (or sine and cosine) integrals, -eta, zeta, digamma, inverse digamma, trigamma, and polygamma functions. +eta, zeta, Lambert's W, digamma, inverse digamma, trigamma, and polygamma functions. Most of these functions were formerly part of Base in early versions of Julia. CI (Linux, macOS, FreeBSD, Windows): diff --git a/docs/src/functions_list.md b/docs/src/functions_list.md index 56d36c45..8757b66c 100644 --- a/docs/src/functions_list.md +++ b/docs/src/functions_list.md @@ -69,4 +69,7 @@ SpecialFunctions.beta SpecialFunctions.logbeta SpecialFunctions.logabsbeta SpecialFunctions.logabsbinomial +SpecialFunctions.lambertw +SpecialFunctions.lambertwbp +SpecialFunctions.LambertW.Ω ``` diff --git a/docs/src/functions_overview.md b/docs/src/functions_overview.md index 01c6089b..350fb671 100644 --- a/docs/src/functions_overview.md +++ b/docs/src/functions_overview.md @@ -95,3 +95,9 @@ Here the *Special Functions* are listed according to the structure of [NIST Digi |:-------- |:----------- | | [`eta(x)`](@ref SpecialFunctions.eta) | [Dirichlet eta function](https://en.wikipedia.org/wiki/Dirichlet_eta_function) at `x` | | [`zeta(x)`](@ref SpecialFunctions.zeta) | [Riemann zeta function](https://en.wikipedia.org/wiki/Riemann_zeta_function) at `x` | + +## [Lambert's W Function](https://dlmf.nist.gov/4.13) +| Function | Description | +|:-------- |:----------- | +| [`lambertw(x, [k=0])`](@ref SpecialFunctions.lambertw) | [Lambert's W function](https://en.wikipedia.org/wiki/Lambert_W_function) at `x` for `k`-th branch | +| [`lambertwbp(x, [k=0])`](@ref SpecialFunctions.lambertwbp) | Accurate value of ``1 + W_k(-\frac{1}{\mathrm{e}} + x)`` for small `x` | \ No newline at end of file diff --git a/src/SpecialFunctions.jl b/src/SpecialFunctions.jl index 9c482350..2725a946 100644 --- a/src/SpecialFunctions.jl +++ b/src/SpecialFunctions.jl @@ -1,5 +1,7 @@ module SpecialFunctions +using Compat + using IrrationalConstants: twoπ, halfπ, @@ -13,6 +15,9 @@ using IrrationalConstants: logπ, log2π +# FIXME temporary until the fate of inve is decided +Base.@irrational inve 0.367879441171442321595 inv(big(ℯ)) + import LogExpFunctions using OpenLibm_jll @@ -77,7 +82,10 @@ export expintx, sinint, cosint, - lbinomial + lbinomial, + lambertw, + lambertwbp, + LambertW include("bessel.jl") include("erf.jl") @@ -90,6 +98,7 @@ include("gamma.jl") include("gamma_inc.jl") include("betanc.jl") include("beta_inc.jl") +include("lambertw.jl") if !isdefined(Base, :get_extension) include("../ext/SpecialFunctionsChainRulesCoreExt.jl") end diff --git a/src/lambertw.jl b/src/lambertw.jl new file mode 100644 index 00000000..a739f0cb --- /dev/null +++ b/src/lambertw.jl @@ -0,0 +1,287 @@ +#### Lambert W function #### + +""" + lambertw(z::Number, k::Integer=0; [maxiter=1000]) + +Compute the `k`th branch of the [Lambert W function](https://en.wikipedia.org/wiki/Lambert_W_function) of `z`. If `z` is real, `k` must be +either `0` or `-1`. For `Real` `z`, the domain of the branch `k = -1` is `[-1/e, 0]` and the +domain of the branch `k = 0` is `[-1/e, Inf]`. For `Complex` `z`, and all `k`, the domain is +the complex plane. + +```jldoctest; setup=:(using SpecialFunctions) +julia> lambertw(-1/ℯ, -1) +-1.0 + +julia> lambertw(-1/ℯ, 0) +-1.0 + +julia> lambertw(0, 0) +0.0 + +julia> lambertw(0, -1) +-Inf + +julia> lambertw(Complex(-10.0, 3.0), 4) +-0.9274337508660128 + 26.37693445371142im +``` +""" +lambertw(z::Number, k::Integer=0; maxiter::Integer=1000) = _lambertw(z, k, maxiter) + +# lambertw(e + 0im, k) is ok for all k +# Maybe this should return a float. But, this should cause no type instability in any case +function _lambertw(::typeof(MathConstants.e), k::Integer, maxits::Integer) + k == 0 && return 1 + throw(DomainError(k)) +end +_lambertw(x::Irrational, k::Integer, maxits::Integer) = _lambertw(float(x), k, maxits) +function _lambertw(x::Union{Integer, Rational}, k::Integer, maxits::Integer) + if k == 0 + x == 0 && return float(zero(x)) + x == 1 && return convert(typeof(float(x)), LambertW.Omega) # must be a more efficient way + end + return _lambertw(float(x), k, maxits) +end + +### Real x + +_lambertw(x::Real, k::Integer, maxits::Integer) = _lambertw(x, Val(Int(k)), maxits) + +# Real x, k = 0 +# There is a magic number here. It could be noted, or possibly removed. +# In particular, the fancy initial condition selection does not seem to help speed. +function _lambertw(x::T, ::Val{0}, maxits::Integer) where T<:Real + isfinite(x) || return x + one_t = one(T) + oneoe = -T(inve) # The branch point + x == oneoe && return -one_t + oneoe < x || throw(DomainError(x)) + itwo_t = 1 / convert(T, 2) + if x > one_t + lx = log(x) + llx = log(lx) + x0 = lx - llx - log(one_t - llx / lx) * itwo_t + else + x0 = (567//1000) * x + end + return lambertw_root_finding(x, x0, maxits) +end + +# Real x, k = -1 +function _lambertw(x::T, ::Val{-1}, maxits::Integer) where T<:Real + oneoe = -T(inve) + x == oneoe && return -one(T) # W approaches -1 as x -> -1/e from above + oneoe < x || throw(DomainError(x)) # branch domain exludes x < -1/e + x == zero(T) && return -convert(T, Inf) # W decreases w/o bound as x -> 0 from below + x < zero(T) || throw(DomainError(x)) + return lambertw_root_finding(x, log(-x), maxits) +end + +_lambertw(x::Real, k::Val, maxits::Integer) = + throw(DomainError(x, "lambertw: for branch k=$k not defined, real x must have branch k == 0 or k == -1")) + +### Complex z + +_lambertw(z::Complex{<:Integer}, k::Integer, maxits::Integer) = _lambertw(float(z), k, maxits) +# choose initial value inside correct branch for root finding +function _lambertw(z::Complex{T}, k::Integer, maxits::Integer) where T<:Real + local w::Complex{T} + pointseven = 7//10 + if abs(z) <= T(inve) + if z == 0 + k == 0 && return z + return complex(-convert(T, Inf), zero(T)) + end + if k == 0 + w = z + elseif k == -1 && imag(z) == 0 && real(z) < 0 + w = complex(log(-real(z)), 1//10^7) # need offset for z ≈ -1/e. + else + w = log(z) + k != 0 ? w += complex(0, k * 2 * pi) : nothing + end + elseif k == 0 && imag(z) <= pointseven && abs(z) <= pointseven + w = abs(z+ 1//2) < 1//10 ? imag(z) > 0 ? + complex(pointseven, pointseven) : + complex(pointseven, -pointseven) : z + else + if real(z) == convert(T, Inf) + k == 0 && return z + return z + complex(0, 2*k*pi) + end + real(z) == -convert(T, Inf) && return -z + complex(0, (2*k+1)*pi) + w = log(z) + k != 0 ? w += complex(0, 2*k*pi) : nothing + end + return lambertw_root_finding(z, w, maxits) +end + +### root finding, iterative solution + +# Use Halley's root-finding method to find +# x = lambertw(z) with initial point x0. +function lambertw_root_finding(z::T, x0::T, maxits::Integer) where T <: Number + x = x0 + lastx = x + lastdiff = zero(real(T)) + converged = false + for _ in 1:maxits + ex = exp(x) + xexz = x * ex - z + x1 = x + 1 + x -= xexz / (ex * x1 - (x + 2) * xexz / (2 * x1)) + xdiff = abs(lastx - x) + if xdiff <= 3 * eps(lastdiff) || lastdiff == xdiff # second condition catches two-value cycle + converged = true + break + end + lastx = x + lastdiff = xdiff + end + converged || @warn "lambertw($z) did not converge in $maxits iterations." + return x +end + +### Lambert's Omega constant + +# compute BigFloat Omega constant at arbitrary precision +function compute_lambertw_Omega() + oc = BigFloat("0.5671432904097838729999686622103555497538157871865125081351310792230457930866845666932194") + precision(oc) <= 256 && return oc + # iteratively improve the precision + # see https://en.wikipedia.org/wiki/Omega_constant#Computation + myeps = eps(BigFloat) + for _ in 1:1000 + nextoc = (1 + oc) / (1 + exp(oc)) + abs(oc - nextoc) <= myeps && return oc + oc = nextoc + end + @warn "Omega precision is less than current BigFloat precision ($(precision(BigFloat)))" + return oc +end + +# "private" declaration of Omega constant +Base.@irrational lambertw_Omega 0.567143290409783872999968662210355 compute_lambertw_Omega() + +module LambertW + +""" +Lambert's Omega (*Ω*) constant. + +Lambert's *Ω* is the solution to *W(Ω) = 1* equation, +where *W(t) = t exp(t)* is the +[Lambert's *W* function](https://en.wikipedia.org/wiki/Lambert_W_function). + +# See also + * https://en.wikipedia.org/wiki/Omega_constant + * [`lambertw()`][@ref SpecialFunctions.lambertw] +""" +const Ω = Irrational{:lambertw_Omega}() +const Omega = Ω # ASCII alias + +end + +### Expansion about branch point x = -1/e + +# Refer to the paper "On the Lambert W function". In (4.22) +# coefficients μ₀ through μ₃ are given explicitly. Recursion relations +# (4.23) and (4.24) for all μ are also given. This code implements the +# recursion relations. + +# We plug the known value μ₂ == -1//3 for (4.22) into (4.23) and +# solve for α₂. We get α₂ = 0. +# Compute array of coefficients μ in (4.22). +# m[1] is μ₀ +function lambertw_coeffs(T::Type{<:Number}, n::Integer) + a = Vector{T}(undef, n) + m = Vector{T}(undef, n) + + a[1] = 2 # α₀ literal in paper + a[2] = -1 # α₁ literal in paper + a[3] = 0 # α₂ get this by solving (4.23) for alpha_2 with values printed in paper + m[1] = -1 # μ₀ literal in paper + m[2] = 1 # μ₁ literal in paper + m[3] = -1//3 # μ₂ literal in paper, but only in (4.22) + + for i in 4:n + # (4.24) + msum = zero(T) + @inbounds for j in 2:(i - 2) + msum += m[j + 1] * m[i + 1 - j] + end + a[i] = msum + + # (4.23) + it = convert(T, i) + m[i] = (it - 2) / it *(m[i - 2] / 2 + a[i - 2] / 4) - + a[i] / 2 - m[i - 1] / it + end + return m +end + +const LAMBERTW_COEFFS_FLOAT64 = lambertw_coeffs(Float64, 500) + +(lambertwbp_evalpoly(x::T, ::Val{N})::T) where {T<:Number, N} = + # assume that Julia compiler is smart to decide for which N to unroll at compile time + # note that we skip μ₀=-1 + evalpoly(x, ntuple(i -> LAMBERTW_COEFFS_FLOAT64[i+1], N-1))*x + +# how many coefficients of the series to use +# to converge to Float64 precision for given x +# We could get finer tuning by separating k=0, -1 branches. +function lambertwbp_series_length(x::Real) + x < 4e-11 && return 3 + # Why N = 5 is omitted? + x < 1e-5 && return 7 + x < 1e-3 && return 12 + x < 1e-2 && return 19 + x < 3e-2 && return 26 + x < 5e-2 && return 32 + x < 1e-1 && return 50 + x < 1.9e-1 && return 100 + x > typeof(x)(inve) && throw(DomainError(x)) # radius of convergence + return 290 # good for x approx .32 +end + +# These may need tuning. +lambertwbp_series_length(z::Complex) = lambertwbp_series_length(abs(z)) + +# p is the argument to the series which is computed from x, +# see `_lambertwbp()`. +lambertwbp_series(p::Number, x::Number) = + lambertwbp_evalpoly(p, Val{lambertwbp_series_length(x)}()) + +_lambertwbp(x::Number, ::Val{0}) = + lambertwbp_series(sqrt(2 * MathConstants.e * x), x) + +_lambertwbp(x::Number, ::Val{-1}) = + lambertwbp_series(-sqrt(2 * MathConstants.e * x), x) + +_lambertwbp(_::Number, k::Val) = + throw(ArgumentError("lambertw() expansion about branch point for k=$k not implemented (only implemented for 0 and -1).")) + +""" + lambertwbp(z, k=0) + +Compute accurate value of `1 + W(-1/e + z)`, for `abs(z)` in `[0, 1/e]` for `k` either `0` or `-1`. +This function is faster and more accurate near the branch point `-1/e` between `k=0` and `k=1`. +The result is accurate to Float64 precision for abs(z) < 0.32. +If `k=-1` and `imag(z) < 0`, the value on the branch `k=1` is returned. + +# Example +```jldoctest; setup=:(using SpecialFunctions) +julia> lambertw(-1/ℯ + 1e-18, -1) +-1.0 + +julia> lambertwbp(1e-18, -1) +-2.331643983409312e-9 + +julia> convert(Float64, (lambertw(-big(1)/ℯ + big(10)^(-18), -1) + 1)) # Same result, but 1000 times slower +-2.331643983409312e-9 +``` + +!!! note + `lambertwbp` uses a series expansion about the branch point `z=-1/ℯ` to avoid loss of precision. + The loss of precision in `lambertw` is analogous to the loss of precision + in computing the `sqrt(1-x)` for `x` close to `1`. +""" +lambertwbp(x::Number, k::Integer=0) = _lambertwbp(x, Val(Int(k))) diff --git a/test/lambertw.jl b/test/lambertw.jl new file mode 100644 index 00000000..354da1ed --- /dev/null +++ b/test/lambertw.jl @@ -0,0 +1,153 @@ +### domain errors +using SpecialFunctions: inve # FIXME temporary until the fate of inve is decided +#using IrrationalConstants + +@test_throws DomainError lambertw(-2.0, 0) +@test_throws DomainError lambertw(-2.0, -1) +@test_throws DomainError lambertw(-2.0, 1) +@test isnan(@inferred(lambertw(NaN))) + +## math constant e +@test_throws DomainError lambertw(MathConstants.e, 1) +@test_throws DomainError lambertw(MathConstants.e, -1) +@test_throws DomainError lambertw(.3, 2) + +## integer arguments return floating point types +@test @inferred(lambertw(0)) isa AbstractFloat +@test @inferred(lambertw(0)) == 0 + +### math constant, MathConstants.e e + +# could return math const e, but this would break type stability +@test @inferred(lambertw(1)) isa AbstractFloat +@test @inferred(lambertw(1)) == float(LambertW.Omega) +@test @inferred(lambertw(big(1))) == big(LambertW.Omega) +@test @inferred(lambertw(MathConstants.e, 0)) == 1 + +## value at branch point where real branches meet +@test lambertw(-inve, 0) == lambertw(-inve, -1) == -1 +@test typeof(lambertw(-inve, 0)) == typeof(lambertw(-inve, -1)) <: AbstractFloat + +## convert irrationals to float + +@test @inferred(lambertw(pi)) ≈ 1.0736581947961492 +@test @inferred(lambertw(pi, 0)) ≈ 1.0736581947961492 + +### infinite args or return values + +@test @inferred(lambertw(0, -1)) == @inferred(lambertw(0.0, -1)) == -Inf +@test @inferred(lambertw(Inf, 0)) == Inf +@test @inferred(lambertw(complex(Inf, 1), 0)) == complex(Inf, 1) +@test @inferred(lambertw(complex(Inf, 0), 1)) == complex(Inf, 2pi) +@test @inferred(lambertw(complex(-Inf, 0), 1)) == complex(Inf, 3pi) +@test @inferred(lambertw(complex(0.0, 0.0), -1)) == complex(-Inf, 0.0) + +## default branch is k = 0 +@test @inferred(lambertw(1.0)) == @inferred(lambertw(1.0, 0)) + +## BigInt args return BigFloats +@test @inferred(lambertw(BigInt(0))) isa BigFloat +@test @inferred(lambertw(BigInt(3))) isa BigFloat + +## Any Integer type allowed for second argument +@test lambertw(-0.2, -1) == lambertw(-0.2, BigInt(-1)) + +## BigInt for second arg does not promote the type +@test typeof(lambertw(-0.2, -1)) == typeof(lambertw(-0.2, BigInt(-1))) + +for (z, k, res) in [(0, 0 , 0), (complex(0, 0), 0 , 0), + (complex(0.0, 0), 0 , 0), (complex(1.0, 0), 0, 0.567143290409783873)] + if Int != Int32 + @test lambertw(z, k) ≈ res + @test lambertw(z) ≈ res + else + @test lambertw(z, k) ≈ res rtol=1e-14 + @test lambertw(z) ≈ res rtol=1e-14 + end +end + +@testset "complex z=$z, k=$k" for (z, k) in + ((complex(1, 1), 2), (complex(1, 1), 0), (complex(.6, .6), 0), + (complex(.6, -.6), 0)) + w = lambertw(z, k) + @test w*exp(w) ≈ z atol=1e-15 +end + +@test lambertw(complex(-3.0, -4.0), 0) ≈ Complex(1.075073066569255, -1.3251023817343588) atol=1e-14 +@test lambertw(complex(-3.0, -4.0), 1) ≈ Complex(0.5887666813694675, 2.7118802109452247) atol=1e-14 +@test lambertw(complex(.3, .3)) ≈ Complex(0.26763519642648767, 0.1837481231767825) + +# test maxiter keyword and convergence warning +@test_logs (:warn, "lambertw(-0.2) did not converge in 3 iterations.") @inferred(lambertw(-0.2, -1, maxiter=3)) +@test lambertw(-0.2, -1, maxiter=5) == lambertw(-0.2, -1) +@test_logs (:warn, "lambertw(0.3 + 0.3im) did not converge in 3 iterations.") @inferred(lambertw(complex(.3, .3), maxiter=3)) +@test lambertw(complex(.3, .3), maxiter=5) == lambertw(complex(.3, .3)) + +# bug fix +# The routine will start at -1/e + eps * im, rather than -1/e + 0im, +# otherwise root finding will fail +@test lambertw(-inve + 0im, -1) ≈ -1 atol=1e-7 + +# lambertw for BigFloat is more precise than Float64. Note +# that 70 digits in test is about 35 digits in W +@testset "lambertw() for BigFloat z=$z" for z in + [BigFloat(1), BigFloat(2), complex(BigFloat(1), BigFloat(1))] + W = lambertw(z) + @test z ≈ W*exp(W) atol=BigFloat(10)^(-70) +end + +@testset "LambertW.Omega" begin + @test isapprox(LambertW.Ω * exp(LambertW.Ω), 1) + @test LambertW.Omega === LambertW.Ω + + # lower than default precision + setprecision(BigFloat, 196) do + o = big(LambertW.Ω) + @test precision(o) == 196 + @test isapprox(o * exp(o), 1, atol=eps(BigFloat)) + + oalias = big(LambertW.Omega) + @test o == oalias + end + + # higher than default precision + setprecision(BigFloat, 2048) do + o = big(LambertW.Ω) + @test precision(o) == 2048 + @test isapprox(o * exp(o), 1, atol=eps(BigFloat)) + + oalias = big(LambertW.Omega) + @test o == oalias + end +end + +### expansion about branch point +@testset "lambertwbp()" begin + # not a domain error, but not implemented + @test_throws ArgumentError lambertwbp(1, 1) + @test_throws ArgumentError lambertwbp(inve + 1e-5, 2) + @test_throws DomainError lambertwbp(inve + 1e-5, 0) + @test_throws DomainError lambertwbp(inve + 1e-5, -1) + + # Expansions about branch point converges almost to machine precision + # except near the radius of convergence. + # Complex args are not tested here. + + @testset "double-precision expansion near branch point using BigFloats" begin + setprecision(2048) do + z = BigFloat(10)^(-12) + for _ in 1:300 + @test lambertwbp(Float64(z)) ≈ 1 + lambertw(z - big(inve)) atol=5e-16 + @test lambertwbp(Float64(z), -1) ≈ 1 + lambertw(z - big(inve), -1) atol=1e-15 + + z *= 1.1 + if z > 0.23 break end + end + end + end + + # test the expansion about branch point for k=-1, + # by comparing to exact BigFloat calculation. + @test @inferred(lambertwbp(1e-20, -1)) ≈ 1 + lambertw(-big(inve) + BigFloat(10)^(-20), -1) atol=1e-16 + @test @inferred(lambertwbp(Complex(.01, .01), -1)) ≈ Complex(-0.27550382080412062443536, -0.12778889284946406573511) atol=1e-16 +end diff --git a/test/runtests.jl b/test/runtests.jl index 69dcafe7..f9d27d31 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -32,6 +32,7 @@ tests = [ "gamma", "logabsgamma", "sincosint", + "lambertw", "other_tests", "chainrules" ]