-
Notifications
You must be signed in to change notification settings - Fork 1
Expand file tree
/
Copy pathtinyfloat.cpp
More file actions
272 lines (221 loc) · 9.16 KB
/
Copy pathtinyfloat.cpp
File metadata and controls
272 lines (221 loc) · 9.16 KB
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
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
#include <bit>
#include <cassert>
#include "tinyfloat.h"
#include "printer.h"
TinyFloat::TinyFloat(bool negative, int16_t exponent, uint32_t mantissa) : negative(negative), exponent(exponent), mantissa(mantissa) {}
TinyFloat::TinyFloat(int i) {
if (!i) {
*this = TinyFloat::zero();
return;
}
negative = i < 0;
exponent = 23;
mantissa = i < 0 ? -i : i;
while (mantissa < (1u<<23)) { // no need to check for exponent > -126
mantissa *= 2;
exponent--;
}
while (mantissa >= (1u<<24)) { // no need to check for exponent < 128
mantissa /= 2;
exponent++;
}
}
TinyFloat::TinyFloat(float f) { // nan/inf are correctly handled
const uint32_t u = std::bit_cast<uint32_t>(f);
uint32_t sign_bit = (u >> 31) % 2;
uint32_t raw_exponent = (u >> 23) % 256;
uint32_t raw_mantissa = u % (1u<<23);
negative = sign_bit;
exponent = raw_exponent - 127;
mantissa = raw_mantissa;
if (exponent==-127) // zero or subnormal
exponent++;
else if (exponent<128) // normal, recover the hidden bit = 1
mantissa = raw_mantissa + (1u<<23);
}
TinyFloat::operator float() const { // nan/inf are correctly handled
uint32_t sign_bit = negative;
uint32_t raw_exponent = exponent+127;
uint32_t raw_mantissa = mantissa % (1u<<23); // clear the hidden bit
if (exponent==-126 && mantissa<(1u<<23))
raw_exponent = 0; // zero or subnormal
return std::bit_cast<float>((sign_bit<<31) + (raw_exponent<<23) + raw_mantissa);
}
std::ostream& operator<<(std::ostream& out, const TinyFloat& f) {
if (f.isnan()) {
out << "nan";
} else {
if (f.negative) out << "-";
if (f.isinf()) out << "inf";
else out << Q128_149(126 + f.exponent, f.mantissa);
}
return out;
}
bool operator==(const TinyFloat& lhs, const TinyFloat& rhs) {
if (lhs.isnan() || rhs.isnan()) return false; // NaNs are unordered
if (lhs.isfinite() && rhs.isfinite() && !lhs.mantissa && !rhs.mantissa) return true; // +0 = -0
return lhs.mantissa == rhs.mantissa && lhs.exponent == rhs.exponent && lhs.negative == rhs.negative;
}
bool operator!=(const TinyFloat& lhs, const TinyFloat& rhs) {
return !(lhs == rhs);
}
bool operator<(const TinyFloat& lhs, const TinyFloat& rhs) {
if (lhs.isnan() || rhs.isnan() || lhs==rhs) return false;
if (lhs.negative != rhs.negative) // positive > negative
return lhs.negative;
return lhs.negative != // same sign and not equal
((lhs.exponent < rhs.exponent) || // => check exponents and then mantissas
(lhs.exponent == rhs.exponent && lhs.mantissa < rhs.mantissa));
}
bool operator>(const TinyFloat& lhs, const TinyFloat& rhs) {
if (lhs.isnan() || rhs.isnan()) return false; // NaNs are unordered
return !(lhs<rhs || lhs==rhs);
}
bool operator<=(const TinyFloat& lhs, const TinyFloat& rhs) {
return lhs<rhs || lhs==rhs;
}
bool operator>=(const TinyFloat& lhs, const TinyFloat& rhs) {
return lhs>rhs || lhs==rhs;
}
TinyFloat operator+(const TinyFloat &lhs, const TinyFloat &rhs) {
TinyFloat a = lhs;
TinyFloat b = rhs;
if (a.isnan() || b.isnan())
return TinyFloat::nan();
if (a.isinf() && b.isinf()) {
if (a.negative == b.negative) return a; // same sign infinity
return TinyFloat::nan(); // inf + -inf = nan
}
if (a.isinf()) return a;
if (b.isinf()) return b;
if (!a.mantissa && !b.mantissa) // handle zeros
return TinyFloat::zero(a.negative && b.negative); // if signs differ, result is +0
if (a.exponent < b.exponent)
std::swap(a, b);
a.mantissa *= 8; // reserve place for GRS bits
b.mantissa *= 8;
while (a.exponent > b.exponent) { // align exponents
b.mantissa = (b.mantissa/2) | (b.mantissa%2); // LSB is sticky
b.exponent++;
}
TinyFloat sum = { a.mantissa >= b.mantissa ? a.negative : b.negative, a.exponent, 0 };
if (a.negative == b.negative)
sum.mantissa = a.mantissa + b.mantissa;
else
if (a.mantissa >= b.mantissa)
sum.mantissa = a.mantissa - b.mantissa;
else
sum.mantissa = b.mantissa - a.mantissa;
while (sum.mantissa < (1u<<(23+3)) && sum.exponent > -126) { // normalize the result
sum.mantissa *= 2;
sum.exponent--;
}
while (sum.mantissa >= (1u<<(24+3))) { // can't be more than one iteration
sum.mantissa = (sum.mantissa/2) | (sum.mantissa%2); // do not forget the sticky bit
sum.exponent++;
}
uint32_t g = (sum.mantissa / 4) % 2; // guard bit
uint32_t r = (sum.mantissa / 2) % 2; // round bit
uint32_t s = sum.mantissa % 2; // sticky bit
sum.mantissa /= 8;
if (g && (r || s || (sum.mantissa % 2))) { // round-to-nearest, even-on-ties
sum.mantissa++;
if (sum.mantissa == (1u<<24)) { // renormalize if necessary
sum.mantissa /= 2;
sum.exponent++;
}
}
if (sum.exponent >= 128) // handle overflow
return TinyFloat::inf(sum.negative);
if (!sum.mantissa) // When the sum of two operands with opposite signs (or the difference of two operands with like signs)
return TinyFloat::zero(); // is exactly zero, the sign of that sum (or difference) shall be +0
return sum;
}
TinyFloat operator-(const TinyFloat &lhs, const TinyFloat &rhs) {
TinyFloat f(!rhs.negative, rhs.exponent, rhs.mantissa);
return lhs + f;
}
TinyFloat operator*(const TinyFloat &lhs, const TinyFloat &rhs) {
TinyFloat a = lhs;
TinyFloat b = rhs;
if (a.isnan() || b.isnan())
return TinyFloat::nan();
if (a.isinf() || b.isinf()) {
if ((a.isfinite() && !a.mantissa) || (b.isfinite() && !b.mantissa)) // inf * 0 = nan
return TinyFloat::nan();
return TinyFloat::inf(a.negative != b.negative);
}
if (!a.mantissa || !b.mantissa)
return TinyFloat::zero(a.negative != b.negative);
int16_t exponent = a.exponent + b.exponent + 1; // +1 comes from the separation of a.mantissa * b.mantissa into two 24-bit variables
bool negative = a.negative != b.negative;
uint32_t a_hi = a.mantissa / (1u<<12); // multiply 2 24-bit mantissas
uint32_t a_lo = a.mantissa % (1u<<12); // into two 24-bit halves mantissa, mantissa_low
uint32_t b_hi = b.mantissa / (1u<<12);
uint32_t b_lo = b.mantissa % (1u<<12);
uint32_t hihi = a_hi * b_hi;
uint32_t hilo = a_hi * b_lo;
uint32_t lohi = a_lo * b_hi;
uint32_t lolo = a_lo * b_lo;
uint32_t mantissa_low = lolo + (hilo % (1u<<12) + lohi % (1u<<12)) * (1u<<12);
uint32_t mantissa = hihi + hilo / (1u<<12) + lohi / (1u<<12) + mantissa_low/(1u<<24);
mantissa_low = mantissa_low % (1u<<24);
while (mantissa < (1u<<23) && exponent > -126) { // normalize the result
mantissa = mantissa * 2 + mantissa_low / (1u<<23);
mantissa_low = (mantissa_low * 2) % (1u<<24);
exponent--;
}
while (exponent < -126) {
mantissa_low = ((mantissa_low + (mantissa % 2) * (1u<<24))/2) | (mantissa_low % 2); // LSB is sticky
mantissa /= 2;
exponent++;
}
if (mantissa_low / (1u<<23) && (mantissa_low % (1u<<23) || mantissa % 2)) { // round-to-nearest, even-on-ties
mantissa++;
if (mantissa == (1u<<24)) { // renormalize if necessary
mantissa /= 2;
exponent++;
}
}
if (exponent >= 128) // handle overflow
return TinyFloat::inf(negative);
return { negative, exponent, mantissa };
}
TinyFloat operator/(const TinyFloat &a, const TinyFloat &b) {
if (a.isnan() || b.isnan() || (a.isinf() && b.isinf()) || (!a.mantissa && !b.mantissa))
return TinyFloat::nan();
bool negative = a.negative != b.negative;
if (a.isinf() || !b.mantissa)
return TinyFloat::inf(negative);
if (!a.mantissa || b.isinf())
return TinyFloat::zero(negative);
assert(a.isfinite() && b.isfinite() && a.mantissa && b.mantissa);
uint32_t mantissa = a.mantissa / b.mantissa;
uint32_t remainder = a.mantissa % b.mantissa;
int16_t exponent = a.exponent - b.exponent + 23;
while (mantissa < (1u<<23) && exponent > -126) { // normalize the result
remainder *= 2;
mantissa = mantissa * 2 + remainder / b.mantissa;
remainder = remainder % b.mantissa;
exponent--;
}
while (exponent < -126) {
remainder = (remainder + (mantissa % 2)*b.mantissa)/2 | (remainder % 2); // LSB is sticky
mantissa /= 2;
exponent++;
}
if (remainder*2 > b.mantissa || (remainder*2 == b.mantissa && mantissa % 2)) { // round-to-nearest, even-on-ties
mantissa++;
if (mantissa == (1u<<24)) { // renormalize if necessary
mantissa /= 2;
exponent++;
}
}
if (exponent >= 128) // handle overflow
return TinyFloat::inf(negative);
return { negative, exponent, mantissa };
}
TinyFloat operator-(const TinyFloat &f) {
if (f.isnan()) return f;
return { !f.negative, f.exponent, f.mantissa };
}