gistfile1.txt

import Base.Math: significand_bits,exponent_bias

@inline exp_kernel(x::Float64) = evalpoly(x, (1.0, 1.0, .5, .1666666666666666666))

const J_TABLE = Float64[2^((j-1)/1024) for j in 1:1024]

function myexp(x::T) where T<:Float64
	# Negation so NaN gets caught
    if !(abs(x) < 708.3964185322641)
        x <= -708.3964185322641 && return 0.0
        x >=  709.782712893384 && return Inf
        isnan(x) && return x
    end
    N = unsafe_trunc(Int, x*1477.3197218702985 + .5)   # 1024/ln2
    r = x - N*0.0006769015435155716                    # ln2/1024
    k = N >> 10

    two_to_k = reinterpret(T, (k+exponent_bias(T)) << significand_bits(T))
    return two_to_k * @inbounds J_TABLE[N&1023 + 1] * exp_kernel(r)
end