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
|
package math32
func Exp(x float32) float32 {
if haveArchExp {
return archExp(x)
}
return exp(x)
}
func exp(x float32) float32 {
const (
Ln2Hi = float32(6.9313812256e-01)
Ln2Lo = float32(9.0580006145e-06)
Log2e = float32(1.4426950216e+00)
Overflow = 7.09782712893383973096e+02
Underflow = -7.45133219101941108420e+02
NearZero = 1.0 / (1 << 28) // 2**-28
LogMax = 0x42b2d4fc // The bitmask of log(FLT_MAX), rounded down. This value is the largest input that can be passed to exp() without producing overflow.
LogMin = 0x42aeac50 // The bitmask of |log(REAL_FLT_MIN)|, rounding down
)
// hx := Float32bits(x) & uint32(0x7fffffff)
// special cases
switch {
case IsNaN(x) || IsInf(x, 1):
return x
case IsInf(x, -1):
return 0
case x > Overflow:
return Inf(1)
case x < Underflow:
// case hx > LogMax:
// return Inf(1)
// case x < 0 && hx > LogMin:
return 0
case -NearZero < x && x < NearZero:
return 1 + x
}
// reduce; computed as r = hi - lo for extra precision.
var k int
switch {
case x < 0:
k = int(Log2e*x - 0.5)
case x > 0:
k = int(Log2e*x + 0.5)
}
hi := x - float32(k)*Ln2Hi
lo := float32(k) * Ln2Lo
// compute
return expmulti(hi, lo, k)
}
// Exp2 returns 2**x, the base-2 exponential of x.
//
// Special cases are the same as Exp.
func Exp2(x float32) float32 {
if haveArchExp2 {
return archExp2(x)
}
return exp2(x)
}
func exp2(x float32) float32 {
const (
Ln2Hi = 6.9313812256e-01
Ln2Lo = 9.0580006145e-06
Overflow = 1.0239999999999999e+03
Underflow = -1.0740e+03
)
// special cases
switch {
case IsNaN(x) || IsInf(x, 1):
return x
case IsInf(x, -1):
return 0
case x > Overflow:
return Inf(1)
case x < Underflow:
return 0
}
// argument reduction; x = r×lg(e) + k with |r| ≤ ln(2)/2.
// computed as r = hi - lo for extra precision.
var k int
switch {
case x > 0:
k = int(x + 0.5)
case x < 0:
k = int(x - 0.5)
}
t := x - float32(k)
hi := t * Ln2Hi
lo := -t * Ln2Lo
// compute
return expmulti(hi, lo, k)
}
// exp1 returns e**r × 2**k where r = hi - lo and |r| ≤ ln(2)/2.
func expmulti(hi, lo float32, k int) float32 {
const (
P1 = float32(1.6666667163e-01) /* 0x3e2aaaab */
P2 = float32(-2.7777778450e-03) /* 0xbb360b61 */
P3 = float32(6.6137559770e-05) /* 0x388ab355 */
P4 = float32(-1.6533901999e-06) /* 0xb5ddea0e */
P5 = float32(4.1381369442e-08) /* 0x3331bb4c */
)
r := hi - lo
t := r * r
c := r - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))))
y := 1 - ((lo - (r*c)/(2-c)) - hi)
// TODO(rsc): make sure Ldexp can handle boundary k
return Ldexp(y, k)
}
|
