Merge pull request #17447 from oxinabox/master

Fast exp2(::Int) (fixs #17412)

base/float.jl

@@ -519,9 +519,9 @@ exponent_one(::Type{Float32}) =     0x3f80_0000
 exponent_half(::Type{Float32}) =    0x3f00_0000
 significand_mask(::Type{Float32}) = 0x007f_ffff
 
-significand_bits{T<:AbstractFloat}(::Type{T}) = trailing_ones(significand_mask(T))
-exponent_bits{T<:AbstractFloat}(::Type{T}) = sizeof(T)*8 - significand_bits(T) - 1
-exponent_bias{T<:AbstractFloat}(::Type{T}) = Int(exponent_one(T) >> significand_bits(T))
+@pure significand_bits{T<:AbstractFloat}(::Type{T}) = trailing_ones(significand_mask(T))
+@pure exponent_bits{T<:AbstractFloat}(::Type{T}) = sizeof(T)*8 - significand_bits(T) - 1
+@pure exponent_bias{T<:AbstractFloat}(::Type{T}) = Int(exponent_one(T) >> significand_bits(T))
 
 ## Array operations on floating point numbers ##
 

base/math.jl

@@ -129,6 +129,21 @@ for f in (:sinh, :cosh, :tanh, :atan, :asinh, :exp, :erf, :erfc, :expm1)
     @eval ($f)(x::AbstractFloat) = error("not implemented for ", typeof(x))
 end
 
+# functions with special cases for integer arguments
+@inline function exp2(x::Base.BitInteger)
+    if x > 1023
+        Inf64
+    elseif x <= -1023
+        # if -1073 < x <= -1023 then Result will be a subnormal number
+        # Hex literal with padding must be use to work on 32bit machine
+        reinterpret(Float64, 0x0000_0000_0000_0001  << ((x + 1074)) % UInt)
+    else
+        # We will cast everything to Int64 to avoid errors in case of Int128
+        # If x is a Int128, and is outside the range of Int64, then it is not -1023<x<=1023
+        reinterpret(Float64, (exponent_bias(Float64) + (x % Int64)) << (significand_bits(Float64)) % UInt)
+    end
+end
+
 # TODO: GNU libc has exp10 as an extension; should openlibm?
 exp10(x::Float64) = 10.0^x
 exp10(x::Float32) = 10.0f0^x

test/bigint.jl

@@ -311,3 +311,22 @@ ndigits(rand(big(-999:999)), big(2)^rand(2:999))
 
 @test big(5)^true == big(5)
 @test big(5)^false == one(BigInt)
+
+
+# operations that when applied to Int64 give Float64, should give BigFloat
+@test typeof(exp(a)) == BigFloat
+@test typeof(exp2(a)) == BigFloat
+@test typeof(exp10(a)) == BigFloat
+@test typeof(expm1(a)) == BigFloat
+@test typeof(erf(a)) == BigFloat
+@test typeof(erfc(a)) == BigFloat
+@test typeof(cosh(a)) == BigFloat
+@test typeof(sinh(a)) == BigFloat
+@test typeof(tanh(a)) == BigFloat
+@test typeof(sech(a)) == BigFloat
+@test typeof(csch(a)) == BigFloat
+@test typeof(coth(a)) == BigFloat
+@test typeof(cbrt(a)) == BigFloat
+@test typeof(tan(a)) == BigFloat
+@test typeof(cos(a)) == BigFloat
+@test typeof(sin(a)) == BigFloat

test/math.jl

@@ -192,6 +192,22 @@ end
 @test exp2(Float16(2.)) ≈ exp2(2.)
 @test log(e) == 1
 
+# check exp2(::Integer) matches exp2(::Float)"
+for ii in -2048:2048
+    expected = exp2(float(ii))
+    @test(exp2(Int16(ii)) == expected)
+    @test(exp2(Int32(ii)) == expected)
+    @test(exp2(Int64(ii)) == expected)
+    @test(exp2(Int128(ii)) == expected)
+    if ii >= 0
+        @test(exp2(UInt16(ii)) == expected)
+        @test(exp2(UInt32(ii)) == expected)
+        @test(exp2(UInt64(ii)) == expected)
+        @test(exp2(UInt128(ii)) == expected)
+    end
+end
+
+
 for T in (Int, Float64, BigFloat)
     @test deg2rad(T(180)) ≈ 1pi
     @test deg2rad(T[45, 60]) ≈ [pi/T(4), pi/T(3)]