floats/exp32.go at main · shogo82148/floats · GitHub

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
package floats

// Exp returns e**x, the base-e exponential of a.
//
// Special cases are:
//
//	+Inf.Exp() = +Inf
//	NaN.Exp() = NaN
//
// Very large values overflow to 0 or +Inf.
// Very small values underflow to 1.
func (a Float32) Exp() Float32 {
	const (
		// ln(2) split into high and low parts
		Ln2Hi = Float32(6.9313812256e-01)
		Ln2Lo = Float32(9.0580006145e-06)

		// log2(e)
		Log2e = Float32(1.4426950216e+00)

		// ln(max float32 + 0.5ulp) = ln(2¹²⁷×(2-2⁻²⁴))
		Overflow = Float32(88.7228390818706768)

		// ln(min float32 - 0.5ulp) = ln(2⁻¹⁵⁰)
		Underflow = Float32(-103.972077083991796)

		// The upper limit for underflow
		// when exp(a) ~ 1 + a + a²/2! + ...
		NearZero = Float32(0x1p-14) // 2**-14
	)

	// special cases
	switch {
	case a.IsNaN():
		return NewFloat32NaN()
	case a.IsInf(1):
		return NewFloat32Inf(1)
	case a.IsInf(-1):
		return 0
	case a > Overflow:
		return NewFloat32Inf(1)
	case a < Underflow:
		return 0
	case a.Abs() < NearZero:
		return 1 + a
	}

	// reduce; computed as r = hi - lo for extra precision.
	var k int
	switch {
	case a < 0:
		k = int(Log2e*a - 0.5)
	case a > 0:
		k = int(Log2e*a + 0.5)
	}
	hi := a - Float32(k)*Ln2Hi
	lo := Float32(k) * Ln2Lo

	// compute
	return expmulti32(hi, lo, k)
}

// Exp2 returns 2**x, the base-2 exponential of x.
//
// Special cases are the same as [Exp].
func (a Float32) Exp2() Float32 {
	const (
		// ln(2) split into high and low parts
		Ln2Hi = Float32(6.9313812256e-01)
		Ln2Lo = Float32(9.0580006145e-06)

		// log2(max float32 + 0.5ulp) = log2(2¹²⁷×(2-2⁻²⁴))
		Overflow = Float32(127.99999995700433664361252807573833209430034819546)

		// log2(min float32 - 0.5ulp) = log2(2⁻¹⁵⁰)
		Underflow = Float32(-158)
	)

	// special cases
	switch {
	case a.IsNaN():
		return NewFloat32NaN()
	case a.IsInf(1):
		return NewFloat32Inf(1)
	case a.IsInf(-1):
		return 0
	case a > Overflow:
		return NewFloat32Inf(1)
	case a < 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 a < 0:
		k = int(a - 0.5)
	case a > 0:
		k = int(a + 0.5)
	}
	t := a - Float32(k)
	hi := t * Ln2Hi
	lo := -t * Ln2Lo

	// compute
	return expmulti32(hi, lo, k)
}

func expmulti32(hi, lo Float32, k int) Float32 {
	const (
		P1 = 1.66666666666666657415e-01
		P2 = -2.77777777770155933842e-03
		P3 = 6.61375632143793436117e-05
		P4 = -1.65339022054652515390e-06
		P5 = 4.13813679705723846039e-08
	)

	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)
	return y.Ldexp(k)
}