diff --git a/src/IrrationalConstants.jl b/src/IrrationalConstants.jl index e9c2b27..99ccd5d 100644 --- a/src/IrrationalConstants.jl +++ b/src/IrrationalConstants.jl @@ -26,8 +26,12 @@ export logten, # log(10) logπ, # log(π) log2π, # log(2π) - log4π # log(4π) + log4π, # log(4π) + invℯ, inve, # 1 / ℯ + + LambertW # module that defines Lambert'W Ω: Ω*exp(Ω) = 1 include("stats.jl") +include("lambertw_omega.jl") end # module diff --git a/src/lambertw_omega.jl b/src/lambertw_omega.jl new file mode 100644 index 0000000..71bd088 --- /dev/null +++ b/src/lambertw_omega.jl @@ -0,0 +1,35 @@ +# compute BigFloat Omega constant at arbitrary precision +function compute_lambertw_Omega() + o = BigFloat("0.5671432904097838729999686622103555497538157871865125081351310792230457930866845666932194") + precision(o) <= 256 && return o + # iteratively improve the precision + # see https://en.wikipedia.org/wiki/Omega_constant#Computation + myeps = eps(BigFloat) + for _ in 1:1000 + o_ = (1 + o) / (1 + exp(o)) + abs(o - o_) <= myeps && return o + o = o_ + end + @warn "Omega precision is less than current BigFloat precision ($(precision(BigFloat)))" + return o +end + +@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 Ω = parentmodule(@__MODULE__).lambertw_Omega +const Omega = parentmodule(@__MODULE__).lambertw_Omega # ASCII alias + +end diff --git a/src/stats.jl b/src/stats.jl index c221961..ce1dacc 100644 --- a/src/stats.jl +++ b/src/stats.jl @@ -28,3 +28,6 @@ @irrational logπ 1.1447298858494001741 log(big(π)) @irrational log2π 1.8378770664093454836 log(2 * big(π)) @irrational log4π 2.5310242469692907930 log(4 * big(π)) + +@irrational invℯ 0.367879441171442321595 inv(big(ℯ)) +const inve = invℯ # ASCII alias diff --git a/test/runtests.jl b/test/runtests.jl index 0bab199..684f66d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -40,3 +40,32 @@ end @test isapprox(log(4pi), log4π) end +@testset "1/e" begin + @test isapprox(invℯ, exp(-1)) + @test isapprox(inve, exp(-1)) +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