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
140
141
142
143
144
145
146
147
148
149
150
151
152
|
#include <stdint.h>
#include <string.h>
#include <math.h>
uint32_t float_as_uint32 (float a) {
uint32_t r;
memcpy (&r, &a, sizeof r);
return r;
}
float uint32_as_float (uint32_t a) {
float r;
memcpy (&r, &a, sizeof r);
return r;
}
void logf_ext (float a, float *loghi, float *loglo) {
const float LOG2_HI = 6.93147182e-1f;
const float LOG2_LO = -1.90465421e-9f;
const float SQRT_HALF = 0.70710678f;
float m, r, i, s, t, p, qhi, qlo;
int e;
const float POW_TWO_M23 = 1.19209290e-7f;
const float POW_TWO_P23 = 8388608.0f;
const float FP32_MIN_NORM = 1.175494351e-38f;
i = 0.0f;
if (a < FP32_MIN_NORM) {
a = a * POW_TWO_P23;
i = -23.0f;
}
e = (float_as_uint32 (a) - float_as_uint32 (SQRT_HALF)) & 0xff800000;
m = uint32_as_float (float_as_uint32 (a) - e);
i = fmaf ((float)e, POW_TWO_M23, i);
p = m + 1.0f;
m = m - 1.0f;
r = 1.0f / p;
qhi = r * m;
qlo = r * fmaf (qhi, -m, fmaf (qhi, -2.0f, m));
s = qhi * qhi;
r = 0.1293334961f;
r = fmaf (r, s, 0.1419928074f);
r = fmaf (r, s, 0.2000148296f);
r = fmaf (r, s, 0.3333332539f);
t = fmaf (qhi, qlo + qlo, fmaf (qhi, qhi, -s));
p = s * qhi;
t = fmaf (s, qlo, fmaf (t, qhi, fmaf (s, qhi, -p)));
s = fmaf (r, p, fmaf (r, t, qlo));
r = 2 * qhi;
t = fmaf ( LOG2_HI, i, r);
p = fmaf (-LOG2_HI, i, t);
s = fmaf ( LOG2_LO, i, fmaf (2.f, s, r - p));
*loghi = p = t + s;
*loglo = (t - p) + s;
}
float expf_unchecked (float a) {
float f, j, r;
int i;
j = fmaf (1.442695f, a, 12582912.f) - 12582912.f;
f = fmaf (j, -6.93145752e-1f, a);
f = fmaf (j, -1.42860677e-6f, f);
i = (int)j;
r = 1.37805939e-3f;
r = fmaf (r, f, 8.37312452e-3f);
r = fmaf (r, f, 4.16695364e-2f);
r = fmaf (r, f, 1.66664720e-1f);
r = fmaf (r, f, 4.99999851e-1f);
r = fmaf (r, f, 1.00000000e+0f);
r = fmaf (r, f, 1.00000000e+0f);
float s, t;
uint32_t ia = (i > 0) ? 0u : 0x83000000u;
s = uint32_as_float (0x7f000000u + ia);
t = uint32_as_float (((uint32_t)i << 23) - ia);
r = r * s;
r = r * t;
return r;
}
float powf_core (float a, float b) {
const float MAX_IEEE754_FLT = uint32_as_float (0x7f7fffff);
const float EXP_OVFL_BOUND = 88.7228394f; // 0x1.62e430p+6f;
const float EXP_OVFL_UNFL_F = 104.0f;
const float MY_INF_F = uint32_as_float (0x7f800000);
float lhi, llo, thi, tlo, phi, plo, r;
logf_ext (a, &lhi, &llo);
thi = lhi * b;
if (fabsf (thi) > EXP_OVFL_UNFL_F) {
r = (thi < 0.0f) ? 0.0f : MY_INF_F;
} else {
tlo = fmaf (lhi, b, -thi);
tlo = fmaf (llo, b, +tlo);
phi = thi + tlo;
if (phi == EXP_OVFL_BOUND)
phi = uint32_as_float (float_as_uint32 (phi) - 1);
plo = (thi - phi) + tlo;
r = expf_unchecked (phi);
if (fabsf (r) <= MAX_IEEE754_FLT) {
r = fmaf (plo, r, r);
}
}
return r;
}
float powf (float a, float b) {
const float MY_INF_F = uint32_as_float (0x7f800000);
const float MY_NAN_F = uint32_as_float (0xffc00000);
int expo_odd_int;
float r;
expo_odd_int = fmaf (-2.0f, floorf (0.5f * b), b) == 1.0f;
if ((a == 1.0f) || (b == 0.0f)) {
r = 1.0f;
} else if (isnan (a) || isnan (b)) {
r = a + b;
} else if (isinf (b)) {
r = ((fabsf (a) < 1.0f) != (b < 0.0f)) ? 0.0f : MY_INF_F;
if (a == -1.0f) r = 1.0f;
} else if (isinf (a)) {
r = (b < 0.0f) ? 0.0f : MY_INF_F;
if ((a < 0.0f) && expo_odd_int) r = -r;
} else if (a == 0.0f) {
r = (expo_odd_int) ? (a + a) : 0.0f;
if (b < 0.0f) r = copysignf (MY_INF_F, r);
} else if ((a < 0.0f) && (b != floorf (b))) {
r = MY_NAN_F;
} else {
r = powf_core (fabsf (a), b);
if ((a < 0.0f) && expo_odd_int) {
r = -r;
}
}
return r;
}
double pow (double a, double b) {
return powf(a, b); // who cares about precision in a kernel :3
}
|
