Lambert's Omega constant

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_Ω, lambertw_Omega  # Ω exp(Ω) = 1
 
 include("stats.jl")
+include("lambertw_omega.jl")
 
 end # module

src/lambertw_omega.jl

@@ -0,0 +1,35 @@
+# lazy-initialized LambertW Omega at 256-bit precision
+const lambertw_Omega_BigFloat256 = Ref{BigFloat}()
+
+# compute BigFloat Omega constant at arbitrary precision
+function compute_lambertw_Omega()
+    # initialize lambertw_Omega_BigFloat256
+    isassigned(lambertw_Omega_BigFloat256) ||
+        (lambertw_Omega_BigFloat256[] = BigFloat("0.5671432904097838729999686622103555497538157871865125081351310792230457930866845666932194", 256))
+    o = BigFloat(lambertw_Omega_BigFloat256[], precision(BigFloat)) # initial value with current precision
+    precision(o) <= 256 && return o
+    # iteratively improve the precision of the constant
+    myeps = eps(BigFloat)
+    for _ in 1:100
+        o_ = (1 + o) / (1 + exp(o))
+        abs(o - o_) <= myeps && break
+        o = o_
+    end
+    return o
+end
+
+@irrational lambertw_Ω 0.567143290409783872999968662210355 compute_lambertw_Omega()
+const lambertw_Omega = lambertw_Ω # ASCII alias
+
+"""
+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]
+"""
+lambertw_Ω

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

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