-
Notifications
You must be signed in to change notification settings - Fork 23
sumlog #48
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.
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
base: master
Are you sure you want to change the base?
sumlog #48
Changes from 1 commit
581b9df
5725aa9
77aa3d9
88d6fb1
9ecf589
9db732f
76533e1
5747205
0f5a927
977723d
4d488cd
afa5d94
e400483
cc1aaac
07809b7
16ee153
1af518b
0eaf8d2
1f478d0
0807f7a
2a0004d
e5809d1
eb1b524
6fe8bb1
0def97d
3ad95b2
a0a9348
fa667ec
989a111
a54a024
dc48433
39ca989
207fce2
55d125e
bef4728
9572e48
3848848
c4c3e89
e0f410e
23b5bf1
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,36 @@ | ||
using IrrationalConstants: logtwo | ||
|
||
cscherrer marked this conversation as resolved.
Show resolved
Hide resolved
|
||
""" | ||
$(SIGNATURES) | ||
|
||
Compute `sum(log.(X))`. | ||
|
||
`sum(log.(X))` can be evaluated much more quickly as `sum(log, X)`. However, | ||
this still requires computing `log` for each element of `X`. | ||
|
||
`sumlog(X)` can be faster still, especially an the length of `X` increases. | ||
cscherrer marked this conversation as resolved.
Show resolved
Hide resolved
|
||
|
||
This works by representing the `j`th element of `X` as `xⱼ = aⱼ * 2 ^ bⱼ`, | ||
cscherrer marked this conversation as resolved.
Show resolved
Hide resolved
|
||
allowing us to write | ||
|
||
∑ⱼ log(xⱼ) = log(∏ⱼ aⱼ) + log(2) * ∑ⱼ bⱼ | ||
|
||
cscherrer marked this conversation as resolved.
Show resolved
Hide resolved
|
||
Since `log(2)` is constant, `sumlog` only requires a single `log` evaluation. | ||
cscherrer marked this conversation as resolved.
Show resolved
Hide resolved
|
||
""" | ||
function sumlog(x::AbstractArray{T}) where {T} | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. It seems this will fail for non-floating point types |
||
sig = one(T) | ||
ex = zero(exponent(one(T))) | ||
bound = floatmax(T) / 2 | ||
for xj in x | ||
sig *= significand(xj) | ||
ex += exponent(xj) | ||
cscherrer marked this conversation as resolved.
Show resolved
Hide resolved
|
||
|
||
# Significands are in the rang [1,2), so multiplication will eventually overflow | ||
cscherrer marked this conversation as resolved.
Show resolved
Hide resolved
|
||
if sig > bound | ||
(a, b) = (significand(sig), exponent(sig)) | ||
sig = a | ||
ex += b | ||
end | ||
end | ||
log(sig) + logtwo * ex | ||
end |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,5 @@ | ||
@testset "sumlog" begin | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Some surprises: julia> sumlog([1,2,-0.1,-0.2])
-3.2188758248682
julia> sumlog([1,2,NaN])
ERROR: DomainError with NaN:
Cannot be NaN or Inf.
Stacktrace:
[1] (::Base.Math.var"#throw1#5")(x::Float64)
@ Base.Math ./math.jl:845
[2] exponent
@ ./math.jl:848 [inlined]
julia> sumlog([-0.0])
ERROR: DomainError with -0.0:
Cannot be ±0.0.
julia> sum(log, -0.0)
-Inf There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Calling Adding There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Also, note that tests right now only test Float64 |
||
for x in [10 .* rand(1000), repeat([nextfloat(1.0)], 1000), repeat([prevfloat(2.0)], 1000)] | ||
cscherrer marked this conversation as resolved.
Show resolved
Hide resolved
|
||
@test (@inferred sumlog(x)) ≈ sum(log, x) | ||
end | ||
end |
Uh oh!
There was an error while loading. Please reload this page.