Sum of logs

General Usage Performance
5

Ok now we’re talkin’

julia> function sumlog(x::AbstractArray{T}) where {T}
           sig = one(T)
           ex = zero(exponent(one(T)))
           bound = floatmax(T) / 2
           for xj in x
               sig *= significand(xj) 
               ex += exponent(xj) 
               if sig > bound
                   (a, b) = (significand(sig), exponent(sig))
                   sig = a
                   ex += b
               end
           end
           log(sig) + logtwo * ex
       end
sumlog

julia> x = 10 * rand(100);

julia> @btime sumlog($x)
  127.330 ns (0 allocations: 0 bytes)
129.854

julia> @btime sum(log, $x)
  680.462 ns (0 allocations: 0 bytes)
129.854

julia> sumlog(x) - sum(log ∘ big, x)
2.83327e-15

Oh, but 100 iterations isn’t ever enough to overflow, so we should check for larger arrays:

julia> x = 10 .* rand(100000000);

julia> @time sumlog(x)
  0.128695 seconds (1 allocation: 16 bytes)
1.30265e8

julia> @time sum(log, x)
  0.697396 seconds (1 allocation: 16 bytes)
1.30265e8
14 Likes