Add pure Julia ldexp function

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

base/math.jl

@@ -339,6 +339,9 @@ minmax{T<:AbstractFloat}(x::T, y::T) =
            ifelse((y > x) | (signbit(x) > signbit(y)), (x,y), (y,x)))
 
 
+inttype(::Type{Float64}) = Int64
+inttype(::Type{Float32}) = Int32
+inttype(::Type{Float16}) = Int16
 """
     ldexp(x, n)
 
@@ -349,8 +352,45 @@ 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, q::Integer)
+    n = Int(q)
+    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 = unsigned(leading_zeros(xs) - exponent_bits(T))
+        ys = xs << m
+        xu = ys | (xu & sign_mask(T))
+        k = 1 - signed(m)
+        if n < -50000 # underflow
+            # otherwise may have have an integer underflow in the following k += n computation
+            return flipsign(T(0.0), x)
+        end
+    end
+    k += n
+    if k >= exponent_max(T) # overflow
+        return flipsign(T(Inf), x)
+    end
+    if k > 0 # normal case
+        xu = (xu & ~exponent_mask(T)) | unsigned((k % inttype(T)) << significand_bits(T))
+        return reinterpret(T, xu)
+    else # subnormal case
+        if k <= -significand_bits(T) # underflow
+            if n > 50000 # case of integer overflow in n + k
+                return flipsign(T(Inf), x)
+            else
+                return flipsign(T(0.0), x)
+            end
+        end
+        k += significand_bits(T)
+        xu = (xu & ~exponent_mask(T)) | unsigned((k % inttype(T)) << significand_bits(T))
+        z = T(1.0)/(1 << significand_bits(T))
+        return z*reinterpret(T, xu)
+    end
+end
+ldexp(x::Float16, q::Integer) = Float16(ldexp(Float32(x), q))
 
 """
     exponent(x) -> Int
@@ -610,7 +650,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