From b1a500ee7f29d2216b3d748d9a30538b4ac06a4d Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Fri, 19 Nov 2021 15:58:24 +0100 Subject: [PATCH 01/11] add exp(-1) constant useful in many situations, also -exp(-1) is the branching point of Lambert's W_0(t) and W_{-1}(t) --- src/IrrationalConstants.jl | 3 ++- src/stats.jl | 2 ++ test/runtests.jl | 3 +++ 3 files changed, 7 insertions(+), 1 deletion(-) diff --git a/src/IrrationalConstants.jl b/src/IrrationalConstants.jl index e9c2b27..54263e0 100644 --- a/src/IrrationalConstants.jl +++ b/src/IrrationalConstants.jl @@ -26,7 +26,8 @@ export logten, # log(10) logπ, # log(π) log2π, # log(2π) - log4π # log(4π) + log4π, # log(4π) + invℯ # 1 / ℯ include("stats.jl") diff --git a/src/stats.jl b/src/stats.jl index c221961..389dc50 100644 --- a/src/stats.jl +++ b/src/stats.jl @@ -28,3 +28,5 @@ @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(ℯ)) diff --git a/test/runtests.jl b/test/runtests.jl index 0bab199..554ec86 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -40,3 +40,6 @@ end @test isapprox(log(4pi), log4π) end +@testset "1/e" begin + @test isapprox(invℯ, exp(-1)) +end From 50056d3538428f7a19d61ce383fae9412aaa4dfd Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Fri, 19 Nov 2021 16:01:20 +0100 Subject: [PATCH 02/11] add Lambert's Omega constant --- src/IrrationalConstants.jl | 5 ++++- src/lambertw_omega.jl | 32 ++++++++++++++++++++++++++++++++ test/runtests.jl | 8 ++++++++ 3 files changed, 44 insertions(+), 1 deletion(-) create mode 100644 src/lambertw_omega.jl diff --git a/src/IrrationalConstants.jl b/src/IrrationalConstants.jl index 54263e0..e9085a8 100644 --- a/src/IrrationalConstants.jl +++ b/src/IrrationalConstants.jl @@ -27,8 +27,11 @@ export logπ, # log(π) log2π, # log(2π) log4π, # log(4π) - invℯ # 1 / ℯ + invℯ, # 1 / ℯ + + LambertW_Ω # Ω 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..31a52e4 --- /dev/null +++ b/src/lambertw_omega.jl @@ -0,0 +1,32 @@ +# 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 Omega_BigFloat256 + isassigned(LambertW_Omega_BigFloat256) || + (LambertW_Omega_BigFloat256[] = BigFloat("0.5671432904097838729999686622103555497538157871865125081351310792230457930866845666932194")) + o = LambertW_Omega_BigFloat256[] # initial value + precision(BigFloat) <= 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() + +""" +Lambert's Omega (Ω) constant, such that Ω exp(Ω) = 1. + +*W(Ω) = 1*, where *W(t) = t exp(t)* is the *Lambert's W function*. + +# See +https://en.wikipedia.org/wiki/Omega_constant +""" +LambertW_Ω + diff --git a/test/runtests.jl b/test/runtests.jl index 554ec86..1c064e4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -43,3 +43,11 @@ end @testset "1/e" begin @test isapprox(invℯ, exp(-1)) end + +@testset "LambertW_Omega" begin + @test isapprox(LambertW_Ω * exp(LambertW_Ω), 1) + setprecision(BigFloat, 2048) do + o = big(LambertW_Ω) + @test isapprox(o * exp(o), 1, atol=eps(BigFloat)) + end +end From 2552035b7d402b458635935243abb5374142b928 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Thu, 30 Dec 2021 22:39:34 +0100 Subject: [PATCH 03/11] add ASCII inve alias --- src/IrrationalConstants.jl | 2 +- src/stats.jl | 1 + test/runtests.jl | 1 + 3 files changed, 3 insertions(+), 1 deletion(-) diff --git a/src/IrrationalConstants.jl b/src/IrrationalConstants.jl index e9085a8..2f8e8c8 100644 --- a/src/IrrationalConstants.jl +++ b/src/IrrationalConstants.jl @@ -27,7 +27,7 @@ export logπ, # log(π) log2π, # log(2π) log4π, # log(4π) - invℯ, # 1 / ℯ + invℯ, inve, # 1 / ℯ LambertW_Ω # Ω exp(Ω) = 1 diff --git a/src/stats.jl b/src/stats.jl index 389dc50..ce1dacc 100644 --- a/src/stats.jl +++ b/src/stats.jl @@ -30,3 +30,4 @@ @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 1c064e4..f1fb4ce 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -42,6 +42,7 @@ end @testset "1/e" begin @test isapprox(invℯ, exp(-1)) + @test isapprox(inve, exp(-1)) end @testset "LambertW_Omega" begin From e22a444630a268bf55f8fbd10c1ec14cb809133c Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Thu, 30 Dec 2021 22:41:08 +0100 Subject: [PATCH 04/11] LambertW => lambertw, add lambertw_Omega alias --- src/IrrationalConstants.jl | 2 +- src/lambertw_omega.jl | 18 +++++++++--------- test/runtests.jl | 11 ++++++++--- 3 files changed, 18 insertions(+), 13 deletions(-) diff --git a/src/IrrationalConstants.jl b/src/IrrationalConstants.jl index 2f8e8c8..89e316c 100644 --- a/src/IrrationalConstants.jl +++ b/src/IrrationalConstants.jl @@ -29,7 +29,7 @@ export log4π, # log(4π) invℯ, inve, # 1 / ℯ - LambertW_Ω # Ω exp(Ω) = 1 + lambertw_Ω, lambertw_Omega # Ω exp(Ω) = 1 include("stats.jl") include("lambertw_omega.jl") diff --git a/src/lambertw_omega.jl b/src/lambertw_omega.jl index 31a52e4..1fa43ae 100644 --- a/src/lambertw_omega.jl +++ b/src/lambertw_omega.jl @@ -1,12 +1,12 @@ # lazy-initialized LambertW Omega at 256-bit precision -const LambertW_Omega_BigFloat256 = Ref{BigFloat}() +const lambertw_Omega_BigFloat256 = Ref{BigFloat}() # compute BigFloat Omega constant at arbitrary precision -function compute_LambertW_Omega() - # initialize Omega_BigFloat256 - isassigned(LambertW_Omega_BigFloat256) || - (LambertW_Omega_BigFloat256[] = BigFloat("0.5671432904097838729999686622103555497538157871865125081351310792230457930866845666932194")) - o = LambertW_Omega_BigFloat256[] # initial value +function compute_lambertw_Omega() + # initialize lambertw_Omega_BigFloat256 + isassigned(lambertw_Omega_BigFloat256) || + (lambertw_Omega_BigFloat256[] = BigFloat("0.5671432904097838729999686622103555497538157871865125081351310792230457930866845666932194")) + o = lambertw_Omega_BigFloat256[] # initial value precision(BigFloat) <= 256 && return o # iteratively improve the precision of the constant myeps = eps(BigFloat) @@ -18,7 +18,8 @@ function compute_LambertW_Omega() return o end -@irrational LambertW_Ω 0.567143290409783872999968662210355 compute_LambertW_Omega() +@irrational lambertw_Ω 0.567143290409783872999968662210355 compute_lambertw_Omega() +const lambertw_Omega = lambertw_Ω # ASCII alias """ Lambert's Omega (Ω) constant, such that Ω exp(Ω) = 1. @@ -28,5 +29,4 @@ Lambert's Omega (Ω) constant, such that Ω exp(Ω) = 1. # See https://en.wikipedia.org/wiki/Omega_constant """ -LambertW_Ω - +lambertw_Ω \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index f1fb4ce..097d4a0 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -45,10 +45,15 @@ end @test isapprox(inve, exp(-1)) end -@testset "LambertW_Omega" begin - @test isapprox(LambertW_Ω * exp(LambertW_Ω), 1) +@testset "lambertw_Omega" begin + @test isapprox(lambertw_Ω * exp(lambertw_Ω), 1) + @test lambertw_Omega == lambertw_Ω + setprecision(BigFloat, 2048) do - o = big(LambertW_Ω) + o = big(lambertw_Ω) @test isapprox(o * exp(o), 1, atol=eps(BigFloat)) + + oalias = big(lambertw_Omega) + @test o == oalias end end From ba45b7b60c08280383fdb5671c86980c6026e81f Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Thu, 30 Dec 2021 22:41:35 +0100 Subject: [PATCH 05/11] improve lambertw_Omega docstring --- src/lambertw_omega.jl | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/lambertw_omega.jl b/src/lambertw_omega.jl index 1fa43ae..852d648 100644 --- a/src/lambertw_omega.jl +++ b/src/lambertw_omega.jl @@ -22,11 +22,14 @@ end const lambertw_Omega = lambertw_Ω # ASCII alias """ -Lambert's Omega (Ω) constant, such that Ω exp(Ω) = 1. +Lambert's Omega (*Ω*) constant. -*W(Ω) = 1*, where *W(t) = t exp(t)* is the *Lambert's W function*. +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 -https://en.wikipedia.org/wiki/Omega_constant +# See also + * https://en.wikipedia.org/wiki/Omega_constant + * [`lambertw()`][@ref SpecialFunctions.lambertw] """ lambertw_Ω \ No newline at end of file From a72e8350022eacd8dc18c49d52346fbfb0ec1aad Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Thu, 30 Dec 2021 22:42:10 +0100 Subject: [PATCH 06/11] lambertw_Omega: better handle BigFloat precision --- src/lambertw_omega.jl | 6 +++--- test/runtests.jl | 12 ++++++++++++ 2 files changed, 15 insertions(+), 3 deletions(-) diff --git a/src/lambertw_omega.jl b/src/lambertw_omega.jl index 852d648..f78fc13 100644 --- a/src/lambertw_omega.jl +++ b/src/lambertw_omega.jl @@ -5,9 +5,9 @@ const lambertw_Omega_BigFloat256 = Ref{BigFloat}() function compute_lambertw_Omega() # initialize lambertw_Omega_BigFloat256 isassigned(lambertw_Omega_BigFloat256) || - (lambertw_Omega_BigFloat256[] = BigFloat("0.5671432904097838729999686622103555497538157871865125081351310792230457930866845666932194")) - o = lambertw_Omega_BigFloat256[] # initial value - precision(BigFloat) <= 256 && return o + (lambertw_Omega_BigFloat256[] = BigFloat("0.5671432904097838729999686622103555497538157871865125081351310792230457930866845666932194", 256)) + o = BigFloat(lambertw_Omega_BigFloat256[]) # initial value with current precision + precision(o) <= 256 && return o # iteratively improve the precision of the constant myeps = eps(BigFloat) for _ in 1:100 diff --git a/test/runtests.jl b/test/runtests.jl index 097d4a0..2190d9c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -49,8 +49,20 @@ end @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) From b4b290883a876595d40ae116b1728a1dd177bb12 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Thu, 30 Dec 2021 22:57:39 +0100 Subject: [PATCH 07/11] try to fix Omega precision for Julia 1.0 --- src/lambertw_omega.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/lambertw_omega.jl b/src/lambertw_omega.jl index f78fc13..3b79c9f 100644 --- a/src/lambertw_omega.jl +++ b/src/lambertw_omega.jl @@ -6,7 +6,7 @@ function compute_lambertw_Omega() # initialize lambertw_Omega_BigFloat256 isassigned(lambertw_Omega_BigFloat256) || (lambertw_Omega_BigFloat256[] = BigFloat("0.5671432904097838729999686622103555497538157871865125081351310792230457930866845666932194", 256)) - o = BigFloat(lambertw_Omega_BigFloat256[]) # initial value with current precision + 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) From ee9180e14fd88f945081f62ce3d369fecf6612c9 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Thu, 30 Dec 2021 23:11:43 +0100 Subject: [PATCH 08/11] Omega: remove BigFloat lazy-init to fix Julia 1.0 --- src/lambertw_omega.jl | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/src/lambertw_omega.jl b/src/lambertw_omega.jl index 3b79c9f..818d38f 100644 --- a/src/lambertw_omega.jl +++ b/src/lambertw_omega.jl @@ -1,12 +1,6 @@ -# 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 + o = BigFloat("0.5671432904097838729999686622103555497538157871865125081351310792230457930866845666932194") precision(o) <= 256 && return o # iteratively improve the precision of the constant myeps = eps(BigFloat) From d5aae009fa913743526bb2622c815d861ce690d5 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Thu, 30 Dec 2021 23:11:58 +0100 Subject: [PATCH 09/11] Omega: ref. for computation method --- src/lambertw_omega.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/lambertw_omega.jl b/src/lambertw_omega.jl index 818d38f..253bf6c 100644 --- a/src/lambertw_omega.jl +++ b/src/lambertw_omega.jl @@ -2,7 +2,8 @@ function compute_lambertw_Omega() o = BigFloat("0.5671432904097838729999686622103555497538157871865125081351310792230457930866845666932194") precision(o) <= 256 && return o - # iteratively improve the precision of the constant + # iteratively improve the precision + # see https://en.wikipedia.org/wiki/Omega_constant#Computation myeps = eps(BigFloat) for _ in 1:100 o_ = (1 + o) / (1 + exp(o)) From 3fc698095fc2e205d82d74d4049c58a89b8f4896 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Fri, 31 Dec 2021 00:11:23 +0100 Subject: [PATCH 10/11] Omega: allow 1000 iters, warn if prec not reached --- src/lambertw_omega.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/lambertw_omega.jl b/src/lambertw_omega.jl index 253bf6c..a77252f 100644 --- a/src/lambertw_omega.jl +++ b/src/lambertw_omega.jl @@ -5,11 +5,12 @@ function compute_lambertw_Omega() # iteratively improve the precision # see https://en.wikipedia.org/wiki/Omega_constant#Computation myeps = eps(BigFloat) - for _ in 1:100 + for _ in 1:1000 o_ = (1 + o) / (1 + exp(o)) - abs(o - o_) <= myeps && break + abs(o - o_) <= myeps && return o o = o_ end + @warn "lambertw_Omega precision is less than current BigFloat precision ($(precision(BigFloat)))" return o end From 1ee42a5b7b020b46f52085ae5b550e476eb26d0f Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 1 Jan 2022 11:15:39 +0100 Subject: [PATCH 11/11] rename to LambertW.Omega --- src/IrrationalConstants.jl | 2 +- src/lambertw_omega.jl | 12 ++++++++---- test/runtests.jl | 14 +++++++------- 3 files changed, 16 insertions(+), 12 deletions(-) diff --git a/src/IrrationalConstants.jl b/src/IrrationalConstants.jl index 89e316c..99ccd5d 100644 --- a/src/IrrationalConstants.jl +++ b/src/IrrationalConstants.jl @@ -29,7 +29,7 @@ export log4π, # log(4π) invℯ, inve, # 1 / ℯ - lambertw_Ω, lambertw_Omega # Ω exp(Ω) = 1 + LambertW # module that defines Lambert'W Ω: Ω*exp(Ω) = 1 include("stats.jl") include("lambertw_omega.jl") diff --git a/src/lambertw_omega.jl b/src/lambertw_omega.jl index a77252f..71bd088 100644 --- a/src/lambertw_omega.jl +++ b/src/lambertw_omega.jl @@ -10,12 +10,13 @@ function compute_lambertw_Omega() abs(o - o_) <= myeps && return o o = o_ end - @warn "lambertw_Omega precision is less than current BigFloat precision ($(precision(BigFloat)))" + @warn "Omega precision is less than current BigFloat precision ($(precision(BigFloat)))" return o end -@irrational lambertw_Ω 0.567143290409783872999968662210355 compute_lambertw_Omega() -const lambertw_Omega = lambertw_Ω # ASCII alias +@irrational lambertw_Omega 0.567143290409783872999968662210355 compute_lambertw_Omega() + +module LambertW """ Lambert's Omega (*Ω*) constant. @@ -28,4 +29,7 @@ where *W(t) = t exp(t)* is the * https://en.wikipedia.org/wiki/Omega_constant * [`lambertw()`][@ref SpecialFunctions.lambertw] """ -lambertw_Ω \ No newline at end of file +const Ω = parentmodule(@__MODULE__).lambertw_Omega +const Omega = parentmodule(@__MODULE__).lambertw_Omega # ASCII alias + +end diff --git a/test/runtests.jl b/test/runtests.jl index 2190d9c..684f66d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -45,27 +45,27 @@ end @test isapprox(inve, exp(-1)) end -@testset "lambertw_Omega" begin - @test isapprox(lambertw_Ω * exp(lambertw_Ω), 1) - @test lambertw_Omega == lambertw_Ω +@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_Ω) + o = big(LambertW.Ω) @test precision(o) == 196 @test isapprox(o * exp(o), 1, atol=eps(BigFloat)) - oalias = big(lambertw_Omega) + oalias = big(LambertW.Omega) @test o == oalias end # higher than default precision setprecision(BigFloat, 2048) do - o = big(lambertw_Ω) + o = big(LambertW.Ω) @test precision(o) == 2048 @test isapprox(o * exp(o), 1, atol=eps(BigFloat)) - oalias = big(lambertw_Omega) + oalias = big(LambertW.Omega) @test o == oalias end end