Skip to content

Lambert's Omega constant #12

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
wants to merge 11 commits into from
Closed

Lambert's Omega constant #12

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
b1a500e
add exp(-1) constant
alyst Nov 19, 2021
50056d3
add Lambert's Omega constant
alyst Nov 19, 2021
2552035
add ASCII inve alias
alyst Dec 30, 2021
e22a444
LambertW => lambertw, add lambertw_Omega alias
alyst Dec 30, 2021
ba45b7b
improve lambertw_Omega docstring
alyst Dec 30, 2021
a72e835
lambertw_Omega: better handle BigFloat precision
alyst Dec 30, 2021
b4b2908
try to fix Omega precision for Julia 1.0
alyst Dec 30, 2021
ee9180e
Omega: remove BigFloat lazy-init to fix Julia 1.0
alyst Dec 30, 2021
d5aae00
Omega: ref. for computation method
alyst Dec 30, 2021
3fc6980
Omega: allow 1000 iters, warn if prec not reached
alyst Dec 30, 2021
1ee42a5
rename to LambertW.Omega
alyst Jan 1, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 5 additions & 1 deletion src/IrrationalConstants.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,8 +26,12 @@ export
logten, # log(10)
logπ, # log(π)
log2π, # log(2π)
log4π # log(4π)
log4π, # log(4π)
invℯ, # 1 / ℯ

LambertW_Ω # Ω exp(Ω) = 1
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We should stick with the convention in this package and in base, and use lowercase names and in particular no camel case. I.e., e.g. lambertwω or just ω. Again, also we might want to replace it with or at least add an ASCII alias.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is it "ω" or "Ω" (I've taken the latter from https://en.wikipedia.org/wiki/Omega_constant)? I think upper/lowercase for constants has a difference (π vs Π).
I also prefer Julian convention to avoid underscore in the names, but in this case mixing latin and greek wω doesn't contribute to readability. What about lambertw_Ω?
I also think ASCII alias (lambertw_Omega) would be very useful, what's the conventional way to define it?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since it is called "Omega constant" I think we should just use Ω (or ω if we want to stick with the convention in Base and this package of only using lowercase names). And correspondingly, the ASCII alias would just be Omega or omega.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In the revised version it's called lambertw_Ω and lambertw_Omega alias is added. I think plain Ω would be confusing, since it's less well known and/or may interfere with variable name in the user code.
lambertwΩ IMO is quite unreadable.

Copy link

@jlapeyre jlapeyre Dec 31, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

All of the constants in MathConstants and IrrationalConstants follow a naming convention. The omega constant should follow the same convention. None of your proposals follow this convention. But, in any case, it should stay in LambertW.jl.


include("stats.jl")
include("lambertw_omega.jl")

end # module
32 changes: 32 additions & 0 deletions src/lambertw_omega.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this enough for higher precision? It seems a bit dangerous to upper bound the number of iterations.

Suggested change
for _ in 1:100
while true

or alternatively move the stopping criterion here.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wiki says (I've just added the reference) that this is the quadratic method, the number of correct digits is doubled at each iteration. So 100 should be safe for any reasonable precision.
while true sounds a bit dangerous.
I can increase it to 1000 and add a warning that the precision was not reached.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe just use the stopping criterion, i.e., while abs(next_o - o) > eps?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I want to avoid dead loop if there's some subtle bug in BigFloat implementation or something like that.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What bug could cause a dead loop? Generally, we should assume that BigFloat does not contain any bugs - and if we notice any, they should be fixed upstream.

BTW just came across the following while true loop, I still think it would be a simple and straightforward implementation here as well: https://github.com/JuliaLang/julia/blob/3d11f7db65a3461320542aee3b0f26619c4e65e3/base/irrationals.jl#L54

o_ = (1 + o) / (1 + exp(o))
abs(o - o_) <= myeps && break
o = o_
end
return o
end
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not sure if it's worth caching the value in LambertW_Omega_BigFloat256. I don't think this is done for other constants? Also it seems it would return a value of a different precision than requested in the case precision(BigFloat) < 256.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The other constants either have corresponding mpfr_const_xxx() functions or simple formulas.
But I agree it would be more straightforward to have just

    o = BigFloat("0.5671432904097838729999686622103555497538157871865125081351310792230457930866845666932194") # initial value with 256-bit precision
    precision(BigFloat) <= 256 && return o

if BigFloat("...") is fast enough.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Now the returned BigFloat constant should respect the current precision better, and the appropriate tests were added.
I still have kept the lambertw_Omega_BigFloat256 as I assume it's faster to reuse the initialized BigFloat than parse it each time.

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is much slower than the code in LambertW.jl. The two should be benchmarked explicitly before a PR would be accepted.

But, I agree with the others that it does not make sense to move the omega constant out of LambertW.jl. It should stay in LambertW.jl and the latter should be moved into SpecialFuncions


@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
"""
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The docstring should follow the conventions and recommendations in the Julia documentation: https://docs.julialang.org/en/v1/manual/documentation/

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you please suggest your variant?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have improved the docstring. But if you think it does not conform to the guidelines, please let me know of the specific changes you want to implement or just fix it yourself.

LambertW_Ω

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change

2 changes: 2 additions & 0 deletions src/stats.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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(ℯ))
11 changes: 11 additions & 0 deletions test/runtests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,3 +40,14 @@ end
@test isapprox(log(4pi), log4π)
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