1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
|
# This file is a part of Julia. License is MIT: https://julialang.org/license
# 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 ``10^x``.
# Examples
```jldoctest
julia> exp10(2)
100.0
julia> exp10(0.2)
1.5848931924611136
```
"""
exp10(x::Real) = exp10(float(x))
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, uinttype(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, uinttype(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
|
