Add pure Julia ldexp function #19491

musm · 2016-12-03T06:35:39Z

Get rid of openlibm ldexp function for a pure Julia implementation, which is also faster.

benchmark (median) using BenchmarkTools (master) and julia 0.5 for Float64

regular branch is
~ 1.5x speedup
subnormal number and subnormal output
~ 3.6x speedup
normal number and subnormal output
~ 2.6x speedup
subnormal number and normal output
~ 13x speedup

Test values for all branches (T = Float64, Float32)

    @test Amal.ldexp(T(0.0), 0)                       === T(0.0)
    @test Amal.ldexp(T(-0.0), 0)                      === T(-0.0)
    @test Amal.ldexp(T(Inf), 1)                       === T(Inf)
    @test Amal.ldexp(T(-Inf), 1)                      === T(-Inf)
    @test Amal.ldexp(T(NaN), 10)                      === T(NaN)
    @test Amal.ldexp(T(0.8), 4)                       === T(12.8)
    @test Amal.ldexp(T(-0.854375), 5)                 === T(-27.34)
    @test Amal.ldexp(T(1.0), 0)                       === T(1.0)
    @test  Amal.ldexp( realmin(T)/10, typemin(Int64)) === T(0.0)
    @test  Amal.ldexp(-realmin(T)/10, typemin(Int64)) === T(-0.0)
    @test Amal.ldexp(T(1.0), typemax(Int64))          === T(Inf)
    @test Amal.ldexp(T(1.0), typemin(Int64))          === T(0.0)

	# test against reference 
    @test  Amal.ldexp(realmin(T), -1)      == Base.ldexp(realmin(T), -1)
    @test  Amal.ldexp(realmin(T), -2)      == Base.ldexp(realmin(T), -2)
    @test  Amal.ldexp(realmin(T)/2, 0)     == Base.ldexp(realmin(T)/2, 0)
    @test  Amal.ldexp(realmin(T)/3, 0)     == Base.ldexp(realmin(T)/3, 0)
    @test  Amal.ldexp(realmin(T)/3, -1)    == Base.ldexp(realmin(T)/3, -1)
    @test  Amal.ldexp(realmin(T)/3, 11)    == Base.ldexp(realmin(T)/3, 11)
    @test  Amal.ldexp(realmin(T)/11, -10)  == Base.ldexp(realmin(T)/11, -10)
    @test  Amal.ldexp(-realmin(T)/11, -10) == Base.ldexp(-realmin(T)/11, -10)

simonbyrne · 2016-12-03T13:16:46Z

base/math.jl

-ldexp(x::Float64,e::Integer) = ccall((:scalbn,libm), Float64, (Float64,Int32), x, Int32(e))
-ldexp(x::Float32,e::Integer) = ccall((:scalbnf,libm), Float32, (Float32,Int32), x, Int32(e))
+function ldexp{T<:AbstractFloat}(x::T, q::Integer)
+ n = Int(q)


We might as well correctly handle values outside of the range of Int.

The reason I have this is if q::Int128 then things are type unstable with k + q . It's not a problem for Int32 or Int64

simonbyrne · 2016-12-03T13:19:14Z

base/math.jl

+ return flipsign(T(0.0), x)
+ end
+ end
+ k += n


What happens if this overflows (e.g. if n = typemax(Int)?)

k will overflow and be negative and line 378 catches that case

We could handle Int128 and BigInt by moving line 353 here, and

if q > 50000 return flipsign(T(Inf), x) elseif q < -50000 return flipsign(T(0), x) end n = q % Int

This would add another branch, is it worth it, because you still have to check that k is not overflowing and underflowing in addition to this, which currently is cleverly handled. This adds, based on benchmarking the various branches, approximately a 1ns penalty.

Currently we do a lot worse and cast to Int32.

I don't necessarily have a problem with this, but question if this is the best way to do it. The quest for the most generic function is satisfying, but it typically entails performance penalties.

Internally for my codes that need this function, I don't use this since it's a bit too slow in practice but instead use the following:
https://github.com/JuliaMath/Amal.jl/blob/master/src/ldexp.jl
it has some important cases where it doesn't pass all cases, but that's not important since those cases are handeled externally.

Thoughts?

In that case, another option would be

if q > typemax(Int) return flipsign(T(Inf), x) elseif q < typemin(Int) return flipsign(T(0), x) end n = q % Int

which should be able to optimise away the check if q::Int.

musm · 2016-12-03T23:41:47Z

base/math.jl

+ ys = xs << m
+ xu = ys | (xu & sign_mask(T))
+ k = 1 - signed(m)
+ if n < -50000 # underflow


This only needs to checked in the subnormal number case, otherwise k > 0 and it's not possible to underflow, thereby eliminating an unecessary branch.

musm · 2016-12-04T19:50:40Z

base/math.jl

+inttype(::Type{Float64}) = Int64
+inttype(::Type{Float32}) = Int32
+inttype(::Type{Float16}) = Int16
+exponent_max{T<:AbstractFloat}(::Type{T}) = Int(exponent_mask(T) >> significand_bits(T))


@simonbyrne what do you think about these?

We might be able to do without inttype (see comment below), but if not perhaps a more descriptive name (maybe intofwidth or something like that?).

an exponent_max function seems reasonable, but as it's really the biased exponent, maybe exponent_raw_max?

Also, it would probably be more appropriate to have these in float.jl rather than here.

Yeah we don't need inttype using your comment below.

Another case where I think this function would be useful is for something like this
intasfloat{T<:IEEEFloat}(::Type{T}, m::Integer) = reinterpret(T, (m % inttype(T)) << significand_bits(T))

but personally I think it's best to not include it, until there's actually another function that would benefit.

for exponent_max I've just inlined it. It seems rather special case to add it as a another function in float.jl. However, if you have objections I don't mind to put it in float.jl or similar

simonbyrne · 2016-12-04T20:34:29Z

base/math.jl

+ return flipsign(T(Inf), x)
+ end
+ if k > 0 # normal case
+ xu = (xu & ~exponent_mask(T)) | unsigned((k % inttype(T)) << significand_bits(T))


why not k % typeof(xu)?

musm · 2016-12-04T22:37:08Z

base/math.jl

+ ys = xs << unsigned(m)
+ xu = ys | (xu & sign_mask(T))
+ k = 1 - m
+ if q < -50000 # underflow


This number is more than safe for up to Float128

musm · 2016-12-05T18:23:50Z

Fixed some typos in the comment and cleaned things up.
@simonbyrne I believe I have addressed all your comments.

simonbyrne · 2016-12-07T21:10:22Z

@nanosoldier runbenchmarks("scalar", vs = ":master")

nanosoldier · 2016-12-07T23:25:57Z

Your benchmark job has completed - possible performance regressions were detected. A full report can be found here. cc @jrevels

musm · 2016-12-08T02:27:06Z

exponent (didn't touch that function) and frexp (same change made here #19512 and no difference on nanosoldier) are totally spurios
cc @jrevels

simonbyrne · 2016-12-08T08:45:56Z

Would you mind adding some tests for Int128 or BigInt (including a case where they cause under/overflow)?

Get rid of openlibm ldexp function for a pure Julia implementation, which is also faster.

musm · 2016-12-08T19:55:51Z

hi @simonbyrne I added tests, thanks.

musm · 2016-12-08T21:55:19Z

I made a very trivial change to using == instead of === for comparison of finite nonzero values. Can we merge this or do we need to wait on tests again?

simonbyrne · 2016-12-08T21:56:50Z

I've already mistakenly merged 1 "trivial" change that broke things in the last 24 hours, so maybe best to wait.

tkelman · 2016-12-08T22:24:13Z

test/math.jl

@@ -72,22 +72,22 @@ end
 @test y == 0

 # more ldexp tests
- @test ldexp(T(0.8), 4) === T(12.8)
- @test ldexp(T(-0.854375), 5) === T(-27.34)
+ @test ldexp(T(0.8), 4) == T(12.8)


Why test less specifically? === verifies that the types are the same. Would it fail here?

no it doesn't fail

then === is better, as == could miss return-type regressions

the rest of the file is inconsistent anyways

ok w/e I dropped the commit, we can now merge

tkelman · 2017-02-02T17:25:22Z

base/math.jl

@@ -349,8 +349,49 @@ julia> ldexp(5., 2)
 20.0
 ```
 """
-ldexp(x::Float64,e::Integer) = ccall((:scalbn,libm), Float64, (Float64,Int32), x, Int32(e))
-ldexp(x::Float32,e::Integer) = ccall((:scalbnf,libm), Float32, (Float32,Int32), x, Int32(e))
+function ldexp{T<:AbstractFloat}(x::T, e::Integer)


this isn't really valid for every single AbstractFloat subtype, is it?

Probably not: realistically it's probably only valid for Float16,Float32 and Float64

I guess it assumes:

Is a bitstype.

Is base-2

Uses the same layout as the IEEE format (sign | biased exponent | significand )

Uses gradual underflow for subnormals

So it could be valid for, e.g. an x87 Float80 or an IEEE Float128, but not BigFloat, IEEE decimal (as in DecFP.jl), DoubleDoubles, IBM floats or formats which do funny things like two's complement exponents.

musm force-pushed the ldexp branch 2 times, most recently from 22c56ef to d680b8f Compare December 3, 2016 06:38

simonbyrne requested changes Dec 3, 2016

View reviewed changes

musm force-pushed the ldexp branch 3 times, most recently from 9b88338 to ce4a913 Compare December 3, 2016 22:48

musm commented Dec 3, 2016

View reviewed changes

musm force-pushed the ldexp branch 3 times, most recently from 3782359 to 928c973 Compare December 4, 2016 19:19

musm commented Dec 4, 2016

View reviewed changes

simonbyrne reviewed Dec 4, 2016

View reviewed changes

musm commented Dec 4, 2016

View reviewed changes

musm force-pushed the ldexp branch from 3fe2640 to 7ceb585 Compare December 5, 2016 18:22

simonbyrne added the domain:maths Mathematical functions label Dec 7, 2016

musm added 3 commits December 8, 2016 14:21

Add pure Julia ldexp function

1507070

Get rid of openlibm ldexp function for a pure Julia implementation, which is also faster.

Improve Integer handling

7020f27

Improve readability

08c15b8

musm force-pushed the ldexp branch from e3d1493 to 9bef2cf Compare December 8, 2016 19:52

Add ldexp tests

3389704

musm force-pushed the ldexp branch from 9bef2cf to 3389704 Compare December 8, 2016 19:56

StefanKarpinski approved these changes Dec 8, 2016

View reviewed changes

simonbyrne approved these changes Dec 8, 2016

View reviewed changes

tkelman reviewed Dec 8, 2016

View reviewed changes

musm force-pushed the ldexp branch from 9b6dc9f to 3389704 Compare December 8, 2016 22:31

simonbyrne merged commit 374c3d6 into JuliaLang:master Dec 8, 2016

musm deleted the ldexp branch December 8, 2016 22:58

martinholters mentioned this pull request Dec 13, 2016

RFC: Make any and all always short-circuit #19543

Merged

tkelman reviewed Feb 2, 2017

View reviewed changes

musm mentioned this pull request Feb 3, 2017

Fix allowed argument types for float manipulation function #20436

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add pure Julia ldexp function #19491

Add pure Julia ldexp function #19491

musm commented Dec 3, 2016 •

edited

Loading

simonbyrne Dec 3, 2016

musm Dec 3, 2016

simonbyrne Dec 3, 2016

musm Dec 3, 2016

simonbyrne Dec 3, 2016 •

edited

Loading

musm Dec 3, 2016

simonbyrne Dec 4, 2016

musm Dec 4, 2016

musm Dec 3, 2016 •

edited

Loading

musm Dec 4, 2016

simonbyrne Dec 4, 2016

musm Dec 4, 2016

simonbyrne Dec 4, 2016

musm Dec 4, 2016

musm Dec 4, 2016 •

edited

Loading

musm commented Dec 5, 2016

simonbyrne commented Dec 7, 2016

nanosoldier commented Dec 7, 2016

musm commented Dec 8, 2016 •

edited

Loading

simonbyrne commented Dec 8, 2016

musm commented Dec 8, 2016

musm commented Dec 8, 2016

simonbyrne commented Dec 8, 2016

tkelman Dec 8, 2016

musm Dec 8, 2016

tkelman Dec 8, 2016

musm Dec 8, 2016

musm Dec 8, 2016

tkelman Feb 2, 2017

simonbyrne Feb 2, 2017

simonbyrne Feb 2, 2017 •

edited

Loading

Add pure Julia ldexp function #19491

Add pure Julia ldexp function #19491

Conversation

musm commented Dec 3, 2016 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

simonbyrne Dec 3, 2016 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

musm Dec 3, 2016 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

musm Dec 4, 2016 • edited Loading

Choose a reason for hiding this comment

musm commented Dec 5, 2016

simonbyrne commented Dec 7, 2016

nanosoldier commented Dec 7, 2016

musm commented Dec 8, 2016 • edited Loading

simonbyrne commented Dec 8, 2016

musm commented Dec 8, 2016

musm commented Dec 8, 2016

simonbyrne commented Dec 8, 2016

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

simonbyrne Feb 2, 2017 • edited Loading

Choose a reason for hiding this comment

musm commented Dec 3, 2016 •

edited

Loading

simonbyrne Dec 3, 2016 •

edited

Loading

musm Dec 3, 2016 •

edited

Loading

musm Dec 4, 2016 •

edited

Loading

musm commented Dec 8, 2016 •

edited

Loading

simonbyrne Feb 2, 2017 •

edited

Loading