Add pure Julia ldexp function (#19491)

* Add pure Julia ldexp function

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

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)
+    xu = reinterpret(Unsigned, x)
+    xs = xu & ~sign_mask(T)
+    xs >= exponent_mask(T) && return x # NaN or Inf
+    k = Int(xs >> significand_bits(T))
+    if k == 0 # x is subnormal
+        xs == 0 && return x # +-0
+        m = leading_zeros(xs) - exponent_bits(T)
+        ys = xs << unsigned(m)
+        xu = ys | (xu & sign_mask(T))
+        k = 1 - m
+        # underflow, otherwise may have integer underflow in the following n + k
+        e < -50000 && return flipsign(T(0.0), x)
+    end
+    # For cases where e of an Integer larger than Int make sure we properly
+    # overlfow/underflow; this is optimized away otherwise.
+    if e > typemax(Int)
+        return flipsign(T(Inf), x)
+    elseif e < typemin(Int)
+        return flipsign(T(0.0), x)
+    end
+    n = e % Int
+    k += n
+    # overflow, if k is larger than maximum posible exponent
+    if k >= Int(exponent_mask(T) >> significand_bits(T))
+        return flipsign(T(Inf), x)
+    end
+    if k > 0 # normal case
+        xu = (xu & ~exponent_mask(T)) | (k % typeof(xu) << significand_bits(T))
+        return reinterpret(T, xu)
+    else # subnormal case
+        if k <= -significand_bits(T) # underflow
+            # overflow, for the case of integer overflow in n + k
+            e > 50000 && return flipsign(T(Inf), x)
+            return flipsign(T(0.0), x)
+        end
+        k += significand_bits(T)
+        z = T(2.0)^-significand_bits(T)
+        xu = (xu & ~exponent_mask(T)) | (k % typeof(xu) << significand_bits(T))
+        return z*reinterpret(T, xu)
+    end
+end
+ldexp(x::Float16, q::Integer) = Float16(ldexp(Float32(x), q))
 
 """
     exponent(x) -> Int
@@ -604,7 +645,6 @@ for func in (:atan2,:hypot)
     end
 end
 
-ldexp(a::Float16, b::Integer) = Float16(ldexp(Float32(a), b))
 cbrt(a::Float16) = Float16(cbrt(Float32(a)))
 
 # More special functions

test/math.jl

@@ -70,12 +70,43 @@ end
         x,y = frexp(convert(T,NaN))
         @test isnan(x)
         @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(1.0), typemax(Int)) === T(Inf)
+        @test ldexp(T(1.0), typemin(Int)) === T(0.0)
+        @test ldexp(prevfloat(realmin(T)), typemax(Int)) === T(Inf)
+        @test ldexp(prevfloat(realmin(T)), typemin(Int)) === T(0.0)
+
+        @test ldexp(T(0.8), Int128(4)) === T(12.8)
+        @test ldexp(T(-0.854375), Int128(5)) === T(-27.34)
+        @test ldexp(T(1.0), typemax(Int128)) === T(Inf)
+        @test ldexp(T(1.0), typemin(Int128)) === T(0.0)
+        @test ldexp(prevfloat(realmin(T)), typemax(Int128)) === T(Inf)
+        @test ldexp(prevfloat(realmin(T)), typemin(Int128)) === T(0.0)
+
+        @test ldexp(T(0.8), BigInt(4)) === T(12.8)
+        @test ldexp(T(-0.854375), BigInt(5)) === T(-27.34)
+        @test ldexp(T(1.0), BigInt(typemax(Int128))) === T(Inf)
+        @test ldexp(T(1.0), BigInt(typemin(Int128))) === T(0.0)
+        @test ldexp(prevfloat(realmin(T)), BigInt(typemax(Int128))) === T(Inf)
+        @test ldexp(prevfloat(realmin(T)), BigInt(typemin(Int128))) === T(0.0)
+
+        # Test also against BigFloat reference. Needs to be exactly rounded.
+        @test ldexp(realmin(T), -1) == T(ldexp(big(realmin(T)), -1))
+        @test ldexp(realmin(T), -2) == T(ldexp(big(realmin(T)), -2))
+        @test ldexp(realmin(T)/2, 0) == T(ldexp(big(realmin(T)/2), 0))
+        @test ldexp(realmin(T)/3, 0) == T(ldexp(big(realmin(T)/3), 0))
+        @test ldexp(realmin(T)/3, -1) == T(ldexp(big(realmin(T)/3), -1))
+        @test ldexp(realmin(T)/3, 11) == T(ldexp(big(realmin(T)/3), 11))
+        @test ldexp(realmin(T)/11, -10) == T(ldexp(big(realmin(T)/11), -10))
+        @test ldexp(-realmin(T)/11, -10) == T(ldexp(big(-realmin(T)/11), -10))
     end
 end
 
 # We compare to BigFloat instead of hard-coding
-# values, assuming that BigFloat has an independent and independently
-# tested implementation.
+# values, assuming that BigFloat has an independently tested implementation.
 @testset "basic math functions" begin
     @testset "$T" for T in (Float32, Float64)
         x = T(1//3)