rename omega => LambertW.Omega and update its def

using the code from JuliaMath/IrrationalConstants.jl/issues/12

No __init__() section is required

Author: alyst

docs/src/functions_list.md

@@ -71,5 +71,5 @@ SpecialFunctions.logabsbeta
 SpecialFunctions.logabsbinomial
 SpecialFunctions.lambertw
 SpecialFunctions.lambertwbp
-SpecialFunctions.omega
+SpecialFunctions.LambertW.Ω
 ```

src/SpecialFunctions.jl

@@ -80,14 +80,8 @@ export
     cosint,
     lbinomial,
     lambertw,
-    lambertwbp
-
-const omega_const_bf_ = Ref{BigFloat}()
-function __init__()
-    # allocate storage for this BigFloat constant each time this module is loaded
-    omega_const_bf_[] =
-        parse(BigFloat,"0.5671432904097838729999686622103555497538157871865125081351310792230457930866845666932194")
-end
+    lambertwbp,
+    LambertW
 
 include("bessel.jl")
 include("erf.jl")

src/lambertw.jl

@@ -37,7 +37,7 @@ _lambertw(x::Irrational, k::Integer, maxits::Integer) = _lambertw(float(x), k, m
 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)), omega) # must be a more efficient way
+        x == 1 && return convert(typeof(float(x)), LambertW.Omega) # must be a more efficient way
     end
     return _lambertw(float(x), k, maxits)
 end
@@ -141,54 +141,44 @@ function lambertw_root_finding(z::T, x0::T, maxits::Integer) where T <: Number
     return x
 end
 
-### omega constant
+### Lambert's 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
-
-# compute omega constant via root finding
-# We could compute higher precision. This converges very quickly.
-function omega_const(::Type{BigFloat})
-    precision(BigFloat) <= 256 && return omega_const_bf_[]
+# 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)
-    oc = omega_const_bf_[]
-    for i in 1:100
+    for _ in 1:1000
         nextoc = (1 + oc) / (1 + exp(oc))
-        abs(oc - nextoc) <= myeps && break
+        abs(oc - nextoc) <= myeps && return oc
         oc = nextoc
     end
+    @warn "Omega precision is less than current BigFloat precision ($(precision(BigFloat)))"
     return oc
 end
 
-"""
-    omega
-    ω
+# "private" declaration of Omega constant
+Base.@irrational lambertw_Omega 0.567143290409783872999968662210355 compute_lambertw_Omega()
 
-The constant defined by `ω exp(ω) = 1`.
+module LambertW
 
-# Example
-```jldoctest
-julia> ω
-ω = 0.5671432904097...
-
-julia> omega
-ω = 0.5671432904097...
+"""
+Lambert's Omega (*Ω*) constant.
 
-julia> ω * exp(ω)
-1.0
+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).
 
-julia> big(omega)
-5.67143290409783872999968662210355549753815787186512508135131079223045793086683e-01
-```
+# See also
+  * https://en.wikipedia.org/wiki/Omega_constant
+  * [`lambertw()`][@ref SpecialFunctions.lambertw]
 """
-const ω = Irrational{:ω}()
-@doc (@doc ω) omega = ω
+const Ω = Irrational{:lambertw_Omega}()
+const Omega = Ω # ASCII alias
 
-Base.Float64(::Irrational{:ω}) = _omega_const
-Base.Float32(::Irrational{:ω}) = Float32(_omega_const)
-Base.Float16(::Irrational{:ω}) = Float16(_omega_const)
-Base.BigFloat(::Irrational{:ω}) = omega_const(BigFloat)
+end
 
 ### Expansion about branch point x = -1/e
 

test/lambertw.jl

@@ -19,6 +19,8 @@ using IrrationalConstants
 
 # 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
@@ -87,19 +89,30 @@ end
     @test z ≈ W*exp(W) atol=BigFloat(10)^(-70)
 end
 
-###  ω constant
+@testset "LambertW.Omega" begin
+    @test isapprox(LambertW.Ω * exp(LambertW.Ω), 1)
+    @test LambertW.Omega === LambertW.Ω
 
-## get ω from recursion and compare to value from lambertw
-let sp = precision(BigFloat)
-    setprecision(512)
-    @test lambertw(big(1)) == big(SpecialFunctions.omega)
-    setprecision(sp)
-end
+    # 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
 
-@test lambertw(1) == float(SpecialFunctions.omega)
-@test convert(Float16, SpecialFunctions.omega) == convert(Float16, 0.5674)
-@test convert(Float32, SpecialFunctions.omega) == 0.56714326f0
-@test lambertw(BigInt(1)) == big(SpecialFunctions.omega)
+    # 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