Add pure julia exp10 function

base/fastmath.jl

@@ -250,7 +250,7 @@ sqrt_fast(x::FloatTypes) = sqrt_llvm_fast(x)
 const libm = Base.libm_name
 
 for f in (:acos, :acosh, :asin, :asinh, :atan, :atanh, :cbrt, :cos,
-          :cosh, :exp2, :exp, :expm1, :lgamma, :log10, :log1p, :log2,
+          :cosh, :exp2, :expm1, :lgamma, :log10, :log1p, :log2,
           :log, :sin, :sinh, :tan, :tanh)
     f_fast = fast_op[f]
     @eval begin
@@ -274,9 +274,6 @@ atan2_fast(x::Float64, y::Float64) =
 # explicit implementations
 
 @fastmath begin
-    exp10_fast(x::T) where {T<:FloatTypes} = exp2(log2(T(10))*x)
-    exp10_fast(x::Integer) = exp10(float(x))
-
     hypot_fast(x::T, y::T) where {T<:FloatTypes} = sqrt(x*x + y*y)
 
     # Note: we use the same comparison for min, max, and minmax, so

base/math.jl

@@ -240,6 +240,7 @@ for f in (:cbrt, :sinh, :cosh, :tanh, :atan, :asinh, :exp2, :expm1)
     end
 end
 exp(x::Real) = exp(float(x))
+exp10(x::Real) = exp10(float(x))
 
 # fallback definitions to prevent infinite loop from $f(x::Real) def above
 
@@ -288,11 +289,6 @@ end
     end
 end
 
-# TODO: GNU libc has exp10 as an extension; should openlibm?
-exp10(x::Float64) = 10.0^x
-exp10(x::Float32) = 10.0f0^x
-exp10(x::Real) = exp10(float(x))
-
 # utility for converting NaN return to DomainError
 # the branch in nan_dom_err prevents its callers from inlining, so be sure to force it
 # until the heuristics can be improved
@@ -959,6 +955,7 @@ cbrt(a::Float16) = Float16(cbrt(Float32(a)))
 
 # More special functions
 include(joinpath("special", "exp.jl"))
+include(joinpath("special", "exp10.jl"))
 include(joinpath("special", "trig.jl"))
 include(joinpath("special", "gamma.jl"))
 

base/special/exp.jl

@@ -32,26 +32,26 @@
 # log(2)
 const LN2 = 6.931471805599453094172321214581765680755001343602552541206800094933936219696955e-01
 # log2(e)
-const LOG2E = 1.442695040888963407359924681001892137426646
+const LOG2_E = 1.442695040888963407359924681001892137426646
 
-# log(2) into upper and lower
+# log(2) into upper and lower bits
 LN2U(::Type{Float64}) = 6.93147180369123816490e-1
 LN2U(::Type{Float32}) = 6.9313812256f-1
 
 LN2L(::Type{Float64}) = 1.90821492927058770002e-10
 LN2L(::Type{Float32}) = 9.0580006145f-6
 
-# max and min arguments for exponential fucntions
-MAXEXP(::Type{Float64}) = 7.09782712893383996732e2 # log 2^1023*(2-2^-52)
-MAXEXP(::Type{Float32}) = 88.72283905206835f0      # log 2^127 *(2-2^-23)
+# max and min arguments
+MAX_EXP(::Type{Float64}) = 7.09782712893383996732e2 # log 2^1023*(2-2^-52)
+MAX_EXP(::Type{Float32}) = 88.72283905206835f0      # log 2^127 *(2-2^-23)
 
 # one less than the min exponent since we can sqeeze a bit more from the exp function
-MINEXP(::Type{Float64}) = -7.451332191019412076235e2 # log 2^-1075
-MINEXP(::Type{Float32}) = -103.97207708f0            # log 2^-150
+MIN_EXP(::Type{Float64}) = -7.451332191019412076235e2 # log 2^-1075
+MIN_EXP(::Type{Float32}) = -103.97207708f0            # log 2^-150
 
 @inline exp_kernel(x::Float64) = @horner(x, 1.66666666666666019037e-1,
-        -2.77777777770155933842e-3, 6.61375632143793436117e-5,
-        -1.65339022054652515390e-6, 4.13813679705723846039e-8)
+    -2.77777777770155933842e-3, 6.61375632143793436117e-5,
+    -1.65339022054652515390e-6, 4.13813679705723846039e-8)
 
 @inline exp_kernel(x::Float32) = @horner(x, 1.6666625440f-1, -2.7667332906f-3)
 
@@ -63,29 +63,33 @@ MINEXP(::Type{Float32}) = -103.97207708f0            # log 2^-150
     exp(x)
 
 Compute the natural base exponential of `x`, in other words ``e^x``.
+
+```jldoctest
+julia> exp(1.0)
+2.718281828459045
+```
 """
 function exp(x::T) where T<:Union{Float32,Float64}
     xa = reinterpret(Unsigned, x) & ~sign_mask(T)
     xsb = signbit(x)
 
     # filter out non-finite arguments
-    if xa > reinterpret(Unsigned, MAXEXP(T))
+    if xa > reinterpret(Unsigned, MAX_EXP(T))
         if xa >= exponent_mask(T)
             xa & significand_mask(T) != 0 && return T(NaN)
             return xsb ? T(0.0) : T(Inf) # exp(+-Inf)
         end
-        x > MAXEXP(T) && return T(Inf)
-        x < MINEXP(T) && return T(0.0)
+        x > MAX_EXP(T) && return T(Inf)
+        x < MIN_EXP(T) && return T(0.0)
     end
-
     # This implementation gives 2.7182818284590455 for exp(1.0) when T ==
     # Float64, which is well within the allowable error; however,
     # 2.718281828459045 is closer to the true value so we prefer that answer,
     # given that 1.0 is such an important argument value.
     if x == T(1.0) && T == Float64
         return 2.718281828459045235360
     end
-
+    # compute approximation
     if xa > reinterpret(Unsigned, T(0.5)*T(LN2)) # |x| > 0.5 log(2)
         # argument reduction
         if xa < reinterpret(Unsigned, T(1.5)*T(LN2)) # |x| < 1.5 log(2)
@@ -99,18 +103,16 @@ function exp(x::T) where T<:Union{Float32,Float64}
                 lo = LN2L(T)
             end
         else
-            n = round(T(LOG2E)*x)
+            n = round(T(LOG2_E)*x)
             k = unsafe_trunc(Int,n)
             hi = muladd(n, -LN2U(T), x)
             lo = n*LN2L(T)
         end
-        r = hi - lo
-
         # compute approximation on reduced argument
+        r = hi - lo
         z = r*r
         p = r - z*exp_kernel(z)
         y = T(1.0) - ((lo - (r*p)/(T(2.0) - p)) - hi)
-
         # scale back
         if k > -significand_bits(T)
             # multiply by 2.0 first to prevent overflow, which helps extends the range
@@ -123,7 +125,7 @@ function exp(x::T) where T<:Union{Float32,Float64}
             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
+        # Taylor approximation for small values: exp(x) ≈ 1.0 + x
         return T(1.0) + x
     else
         # primary range with k = 0, so compute approximation directly

base/special/exp10.jl

@@ -0,0 +1,132 @@
+#  Method
+#  1. Argument reduction: Reduce x to an r so that |r| <= 0.5*log10(2). Given x,
+#     find r and integer k such that
+#
+#                x = k*log10(2) + r,  |r| <= 0.5*log10(2).
+#
+# 2. Approximate exp10(r) by a polynomial on the interval [-0.5*log10(2), 0.5*log10(2)]:
+#
+#           exp10(x) = 1.0 + polynomial(x),
+#
+#    sup norm relative error within the interval of the polynomial approximations:
+#    Float64 : [2.7245504724394698952e-18; 2.7245529895753476720e-18]
+#    Float32 : [9.6026471477842205871e-10; 9.6026560194009888672e-10]
+#
+# 3. Scale back: exp10(x) = 2^k * exp10(r)
+
+# log2(10)
+const LOG2_10 = 3.321928094887362347870319429489390175864831393024580612054756395815934776608624
+# log10(2)
+const LOG10_2 = 3.010299956639811952137388947244930267681898814621085413104274611271081892744238e-01
+# log(10)
+const LN10 = 2.302585092994045684017991454684364207601101488628772976033327900967572609677367
+
+# log10(2) into upper and lower bits
+LOG10_2U(::Type{Float64}) = 3.01025390625000000000e-1
+LOG10_2U(::Type{Float32}) = 3.00781250000000000000f-1
+
+LOG10_2L(::Type{Float64}) = 4.60503898119521373889e-6
+LOG10_2L(::Type{Float32}) = 2.48745663981195213739f-4
+
+# max and min arguments
+MAX_EXP10(::Type{Float64}) = 3.08254715559916743851e2 # log 2^1023*(2-2^-52)
+MAX_EXP10(::Type{Float32}) = 38.531839419103626f0     # log 2^127 *(2-2^-23)
+
+# one less than the min exponent since we can sqeeze a bit more from the exp10 function
+MIN_EXP10(::Type{Float64}) = -3.23607245338779784854769e2 # log10 2^-1075
+MIN_EXP10(::Type{Float32}) = -45.15449934959718f0         # log10 2^-150
+
+@inline exp10_kernel(x::Float64) =
+    @horner(x, 1.0,
+    2.30258509299404590109361379290930926799774169921875,
+    2.6509490552391992146397114993305876851081848144531,
+    2.03467859229323178027470930828712880611419677734375,
+    1.17125514891212478829629617393948137760162353515625,
+    0.53938292928868392106522833273629657924175262451172,
+    0.20699584873167015119932443667494226247072219848633,
+    6.8089348259156870502017966373387025669217109680176e-2,
+    1.9597690535095281527677713029333972372114658355713e-2,
+    5.015553121397981796436571499953060992993414402008e-3,
+    1.15474960721768829356725927226534622604958713054657e-3,
+    1.55440426715227567738830671828509366605430841445923e-4,
+    3.8731032432074128681303432086835414338565897196531e-5,
+    2.3804466459036747669197886523306806338950991630554e-3,
+    9.3881392238209649520573607528461934634833596646786e-5,
+    -2.64330486232183387018679354696359951049089431762695e-2)
+
+@inline exp10_kernel(x::Float32) =
+    @horner(x, 1.0f0,
+    2.302585124969482421875f0,
+    2.650949001312255859375f0,
+    2.0346698760986328125f0,
+    1.17125606536865234375f0,
+    0.5400512218475341796875f0,
+    0.20749187469482421875f0,
+    5.2789829671382904052734375f-2)
+
+@eval exp10_small_thres(::Type{Float64}) = $(2.0^-29)
+@eval exp10_small_thres(::Type{Float32}) = $(2.0f0^-14)
+
+"""
+    exp10(x)
+
+Compute the base `10` exponential of `x`, in other words ``10^x``.
+
+```jldoctest
+julia> exp10(1.0)
+10.0
+```
+"""
+function exp10(x::T) where T<:Union{Float32,Float64}
+    xa = reinterpret(Unsigned, x) & ~sign_mask(T)
+    xsb = signbit(x)
+
+    # filter out non-finite arguments
+    if xa > reinterpret(Unsigned, MAX_EXP10(T))
+        if xa >= exponent_mask(T)
+            xa & significand_mask(T) != 0 && return T(NaN)
+            return xsb ? T(0.0) : T(Inf) # exp10(+-Inf)
+        end
+        x > MAX_EXP10(T) && return T(Inf)
+        x < MIN_EXP10(T) && return T(0.0)
+    end
+    # compute approximation
+    if xa > reinterpret(Unsigned, T(0.5)*T(LOG10_2)) # |x| > 0.5 log10(2).
+        # argument reduction
+        if xa < reinterpret(Unsigned, T(1.5)*T(LOG10_2)) # |x| <= 1.5 log10(2)
+            if xsb
+                k = -1
+                r = LOG10_2U(T) + x
+                r = LOG10_2L(T) + r
+            else
+                k = 1
+                r = x - LOG10_2U(T)
+                r = r - LOG10_2L(T)
+            end
+        else
+            n = round(T(LOG2_10)*x)
+            k = unsafe_trunc(Int,n)
+            r = muladd(n, -LOG10_2U(T), x)
+            r = muladd(n, -LOG10_2L(T), r)
+        end
+        # compute approximation on reduced argument
+        y = exp10_kernel(r)
+        # scale back
+        if k > -significand_bits(T)
+            # multiply by 2.0 first to prevent overflow, extending the range
+            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)
+        end
+    elseif xa < reinterpret(Unsigned, exp10_small_thres(T))  # |x| < exp10_small_thres
+        # Taylor approximation for small values: exp10(x) ≈ 1.0 + log(10)*x
+        return muladd(x, T(LN10), T(1.0))
+    else
+        # primary range with k = 0, so compute approximation directly
+        return exp10_kernel(x)
+    end
+end

test/math.jl

@@ -270,7 +270,32 @@ end
     end
 end
 
-@testset "test abstractarray trig fxns" begin
+@testset "exp10 function" begin
+    @testset "accuracy" begin
+        X = map(Float64, vcat(-10:0.00021:10, -35:0.0023:100, -300:0.001:300))
+        for x in X
+            y, yb = exp10(x), exp10(big(x))
+            @test abs(y-yb) <= 1.2*eps(Float64(yb))
+        end
+        X = map(Float32, vcat(-10:0.00021:10, -35:0.0023:35, -35:0.001:35))
+        for x in X
+            y, yb = exp10(x), exp10(big(x))
+            @test abs(y-yb) <= 1.2*eps(Float32(yb))
+        end
+    end
+    @testset "$T edge cases" for T in (Float64, Float32)
+        @test isnan(exp10(T(NaN)))
+        @test exp10(T(-Inf)) === T(0.0)
+        @test exp10(T(Inf)) === T(Inf)
+        @test exp10(T(0.0)) === T(1.0) # exact
+        @test exp10(T(1.0)) === T(10.0)
+        @test exp10(T(3.0)) === T(1000.0)
+        @test exp10(T(5000.0)) === T(Inf)
+        @test exp10(T(-5000.0)) === T(0.0)
+    end
+end
+
+@testset "test abstractarray trig functions" begin
     TAA = rand(2,2)
     TAA = (TAA + TAA.')/2.
     STAA = Symmetric(TAA)