Skip to content

Commit

Permalink
make several pure-marked math functions actually pure
Browse files Browse the repository at this point in the history
  • Loading branch information
vtjnash committed Feb 23, 2017
1 parent 1c5d733 commit cc7990b
Show file tree
Hide file tree
Showing 3 changed files with 19 additions and 18 deletions.
8 changes: 0 additions & 8 deletions base/float.jl
Original file line number Diff line number Diff line change
Expand Up @@ -812,14 +812,6 @@ fpinttype(::Type{Float64}) = UInt64
fpinttype(::Type{Float32}) = UInt32
fpinttype(::Type{Float16}) = UInt16

@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))
# maximum float exponent
@pure exponent_max{T<:AbstractFloat}(::Type{T}) = Int(exponent_mask(T) >> significand_bits(T)) - exponent_bias(T)
# maximum float exponent without bias
@pure exponent_raw_max{T<:AbstractFloat}(::Type{T}) = Int(exponent_mask(T) >> significand_bits(T))

## TwicePrecision utilities
# The numeric constants are half the number of bits in the mantissa
for (F, T, n) in ((Float16, UInt16, 5), (Float32, UInt32, 12), (Float64, UInt64, 26))
Expand Down
21 changes: 15 additions & 6 deletions base/math.jl
Original file line number Diff line number Diff line change
Expand Up @@ -20,13 +20,23 @@ import Base: log, exp, sin, cos, tan, sinh, cosh, tanh, asin,
max, min, minmax, ^, exp2, muladd, rem,
exp10, expm1, log1p

using Base: sign_mask, exponent_mask, exponent_one, exponent_bias,
exponent_half, exponent_max, exponent_raw_max, fpinttype,
significand_mask, significand_bits, exponent_bits
using Base: sign_mask, exponent_mask, exponent_one,
exponent_half, fpinttype, significand_mask

using Core.Intrinsics: sqrt_llvm

const IEEEFloat = Union{Float16,Float32,Float64}
const IEEEFloat = Union{Float16, Float32, Float64}

for T in (Float16, Float32, Float64)
@eval significand_bits(::Type{$T}) = $(trailing_ones(significand_mask(T)))
@eval exponent_bits(::Type{$T}) = $(sizeof(T)*8 - significand_bits(T) - 1)
@eval exponent_bias(::Type{$T}) = $(Int(exponent_one(T) >> significand_bits(T)))
# maximum float exponent
@eval exponent_max(::Type{$T}) = $(Int(exponent_mask(T) >> significand_bits(T)) - exponent_bias(T))
# maximum float exponent without bias
@eval exponent_raw_max(::Type{$T}) = $(Int(exponent_mask(T) >> significand_bits(T)))
end

# non-type specific math functions

"""
Expand Down Expand Up @@ -229,8 +239,6 @@ for f in (:cbrt, :sinh, :cosh, :tanh, :atan, :asinh, :exp2, :expm1)
($f)(x::Real) = ($f)(float(x))
end
end
# pure julia exp function
include("special/exp.jl")
exp(x::Real) = exp(float(x))

# fallback definitions to prevent infinite loop from $f(x::Real) def above
Expand Down Expand Up @@ -950,6 +958,7 @@ end
cbrt(a::Float16) = Float16(cbrt(Float32(a)))

# More special functions
include("special/exp.jl")
include("special/trig.jl")
include("special/gamma.jl")

Expand Down
8 changes: 4 additions & 4 deletions base/special/exp.jl
Original file line number Diff line number Diff line change
Expand Up @@ -56,8 +56,8 @@ MINEXP(::Type{Float32}) = -103.97207708f0 # log 2^-150
@inline exp_kernel(x::Float32) = @horner(x, 1.6666625440f-1, -2.7667332906f-3)

# for values smaller than this threshold just use a Taylor expansion
exp_small_thres(::Type{Float64}) = 2.0^-28
exp_small_thres(::Type{Float32}) = 2.0f0^-13
@eval exp_small_thres(::Type{Float64}) = $(2.0^-28)
@eval exp_small_thres(::Type{Float32}) = $(2.0f0^-13)

"""
exp(x)
Expand Down Expand Up @@ -114,13 +114,13 @@ function exp{T<:Union{Float32,Float64}}(x::T)
# scale back
if k > -significand_bits(T)
# multiply by 2.0 first to prevent overflow, which helps extends the range
k == exponent_max(T) && return y*T(2.0)*T(2.0)^(exponent_max(T) - 1)
k == exponent_max(T) && return y * T(2.0) * T(2.0)^(exponent_max(T) - 1)
twopk = reinterpret(T, rem(exponent_bias(T) + k, fpinttype(T)) << significand_bits(T))
return y*twopk
else
# add significand_bits(T) + 1 to lift the range outside the subnormals
twopk = reinterpret(T, rem(exponent_bias(T) + significand_bits(T) + 1 + k, fpinttype(T)) << significand_bits(T))
return y*twopk*T(2.0)^(-significand_bits(T) - 1)
return y * twopk * T(2.0)^(-significand_bits(T) - 1)
end
elseif xa < reinterpret(Unsigned, exp_small_thres(T)) # |x| < exp_small_thres
# Taylor approximation for small x
Expand Down

0 comments on commit cc7990b

Please sign in to comment.