1 /*******************************************************************************
2 * Copyright (c) 2019-2020 The Khronos Group Inc.
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 ******************************************************************************/
16
17 /**
18 * This is a header-only utility library that provides OpenCL host code with
19 * routines for converting to/from cl_half values.
20 *
21 * Example usage:
22 *
23 * #include <CL/cl_half.h>
24 * ...
25 * cl_half h = cl_half_from_float(0.5f, CL_HALF_RTE);
26 * cl_float f = cl_half_to_float(h);
27 */
28
29 #ifndef OPENCL_CL_HALF_H
30 #define OPENCL_CL_HALF_H
31
32 #include <CL/cl_platform.h>
33
34 #include <stdint.h>
35
36 #ifdef __cplusplus
37 extern "C" {
38 #endif
39
40
41 /**
42 * Rounding mode used when converting to cl_half.
43 */
44 typedef enum
45 {
46 CL_HALF_RTE, // round to nearest even
47 CL_HALF_RTZ, // round towards zero
48 CL_HALF_RTP, // round towards positive infinity
49 CL_HALF_RTN, // round towards negative infinity
50 } cl_half_rounding_mode;
51
52
53 /* Private utility macros. */
54 #define CL_HALF_EXP_MASK 0x7C00
55 #define CL_HALF_MAX_FINITE_MAG 0x7BFF
56
57
58 /*
59 * Utility to deal with values that overflow when converting to half precision.
60 */
cl_half_handle_overflow(cl_half_rounding_mode rounding_mode,uint16_t sign)61 static inline cl_half cl_half_handle_overflow(cl_half_rounding_mode rounding_mode,
62 uint16_t sign)
63 {
64 if (rounding_mode == CL_HALF_RTZ)
65 {
66 // Round overflow towards zero -> largest finite number (preserving sign)
67 return (sign << 15) | CL_HALF_MAX_FINITE_MAG;
68 }
69 else if (rounding_mode == CL_HALF_RTP && sign)
70 {
71 // Round negative overflow towards positive infinity -> most negative finite number
72 return (1 << 15) | CL_HALF_MAX_FINITE_MAG;
73 }
74 else if (rounding_mode == CL_HALF_RTN && !sign)
75 {
76 // Round positive overflow towards negative infinity -> largest finite number
77 return CL_HALF_MAX_FINITE_MAG;
78 }
79
80 // Overflow to infinity
81 return (sign << 15) | CL_HALF_EXP_MASK;
82 }
83
84 /*
85 * Utility to deal with values that underflow when converting to half precision.
86 */
cl_half_handle_underflow(cl_half_rounding_mode rounding_mode,uint16_t sign)87 static inline cl_half cl_half_handle_underflow(cl_half_rounding_mode rounding_mode,
88 uint16_t sign)
89 {
90 if (rounding_mode == CL_HALF_RTP && !sign)
91 {
92 // Round underflow towards positive infinity -> smallest positive value
93 return (sign << 15) | 1;
94 }
95 else if (rounding_mode == CL_HALF_RTN && sign)
96 {
97 // Round underflow towards negative infinity -> largest negative value
98 return (sign << 15) | 1;
99 }
100
101 // Flush to zero
102 return (sign << 15);
103 }
104
105
106 /**
107 * Convert a cl_float to a cl_half.
108 */
cl_half_from_float(cl_float f,cl_half_rounding_mode rounding_mode)109 static inline cl_half cl_half_from_float(cl_float f, cl_half_rounding_mode rounding_mode)
110 {
111 // Type-punning to get direct access to underlying bits
112 union
113 {
114 cl_float f;
115 uint32_t i;
116 } f32;
117 f32.f = f;
118
119 // Extract sign bit
120 uint16_t sign = f32.i >> 31;
121
122 // Extract FP32 exponent and mantissa
123 uint32_t f_exp = (f32.i >> (CL_FLT_MANT_DIG - 1)) & 0xFF;
124 uint32_t f_mant = f32.i & ((1 << (CL_FLT_MANT_DIG - 1)) - 1);
125
126 // Remove FP32 exponent bias
127 int32_t exp = f_exp - CL_FLT_MAX_EXP + 1;
128
129 // Add FP16 exponent bias
130 uint16_t h_exp = (uint16_t)(exp + CL_HALF_MAX_EXP - 1);
131
132 // Position of the bit that will become the FP16 mantissa LSB
133 uint32_t lsb_pos = CL_FLT_MANT_DIG - CL_HALF_MANT_DIG;
134
135 // Check for NaN / infinity
136 if (f_exp == 0xFF)
137 {
138 if (f_mant)
139 {
140 // NaN -> propagate mantissa and silence it
141 uint16_t h_mant = (uint16_t)(f_mant >> lsb_pos);
142 h_mant |= 0x200;
143 return (sign << 15) | CL_HALF_EXP_MASK | h_mant;
144 }
145 else
146 {
147 // Infinity -> zero mantissa
148 return (sign << 15) | CL_HALF_EXP_MASK;
149 }
150 }
151
152 // Check for zero
153 if (!f_exp && !f_mant)
154 {
155 return (sign << 15);
156 }
157
158 // Check for overflow
159 if (exp >= CL_HALF_MAX_EXP)
160 {
161 return cl_half_handle_overflow(rounding_mode, sign);
162 }
163
164 // Check for underflow
165 if (exp < (CL_HALF_MIN_EXP - CL_HALF_MANT_DIG - 1))
166 {
167 return cl_half_handle_underflow(rounding_mode, sign);
168 }
169
170 // Check for value that will become denormal
171 if (exp < -14)
172 {
173 // Denormal -> include the implicit 1 from the FP32 mantissa
174 h_exp = 0;
175 f_mant |= 1 << (CL_FLT_MANT_DIG - 1);
176
177 // Mantissa shift amount depends on exponent
178 lsb_pos = -exp + (CL_FLT_MANT_DIG - 25);
179 }
180
181 // Generate FP16 mantissa by shifting FP32 mantissa
182 uint16_t h_mant = (uint16_t)(f_mant >> lsb_pos);
183
184 // Check whether we need to round
185 uint32_t halfway = 1 << (lsb_pos - 1);
186 uint32_t mask = (halfway << 1) - 1;
187 switch (rounding_mode)
188 {
189 case CL_HALF_RTE:
190 if ((f_mant & mask) > halfway)
191 {
192 // More than halfway -> round up
193 h_mant += 1;
194 }
195 else if ((f_mant & mask) == halfway)
196 {
197 // Exactly halfway -> round to nearest even
198 if (h_mant & 0x1)
199 h_mant += 1;
200 }
201 break;
202 case CL_HALF_RTZ:
203 // Mantissa has already been truncated -> do nothing
204 break;
205 case CL_HALF_RTP:
206 if ((f_mant & mask) && !sign)
207 {
208 // Round positive numbers up
209 h_mant += 1;
210 }
211 break;
212 case CL_HALF_RTN:
213 if ((f_mant & mask) && sign)
214 {
215 // Round negative numbers down
216 h_mant += 1;
217 }
218 break;
219 }
220
221 // Check for mantissa overflow
222 if (h_mant & 0x400)
223 {
224 h_exp += 1;
225 h_mant = 0;
226 }
227
228 return (sign << 15) | (h_exp << 10) | h_mant;
229 }
230
231
232 /**
233 * Convert a cl_double to a cl_half.
234 */
cl_half_from_double(cl_double d,cl_half_rounding_mode rounding_mode)235 static inline cl_half cl_half_from_double(cl_double d, cl_half_rounding_mode rounding_mode)
236 {
237 // Type-punning to get direct access to underlying bits
238 union
239 {
240 cl_double d;
241 uint64_t i;
242 } f64;
243 f64.d = d;
244
245 // Extract sign bit
246 uint16_t sign = f64.i >> 63;
247
248 // Extract FP64 exponent and mantissa
249 uint64_t d_exp = (f64.i >> (CL_DBL_MANT_DIG - 1)) & 0x7FF;
250 uint64_t d_mant = f64.i & (((uint64_t)1 << (CL_DBL_MANT_DIG - 1)) - 1);
251
252 // Remove FP64 exponent bias
253 int64_t exp = d_exp - CL_DBL_MAX_EXP + 1;
254
255 // Add FP16 exponent bias
256 uint16_t h_exp = (uint16_t)(exp + CL_HALF_MAX_EXP - 1);
257
258 // Position of the bit that will become the FP16 mantissa LSB
259 uint32_t lsb_pos = CL_DBL_MANT_DIG - CL_HALF_MANT_DIG;
260
261 // Check for NaN / infinity
262 if (d_exp == 0x7FF)
263 {
264 if (d_mant)
265 {
266 // NaN -> propagate mantissa and silence it
267 uint16_t h_mant = (uint16_t)(d_mant >> lsb_pos);
268 h_mant |= 0x200;
269 return (sign << 15) | CL_HALF_EXP_MASK | h_mant;
270 }
271 else
272 {
273 // Infinity -> zero mantissa
274 return (sign << 15) | CL_HALF_EXP_MASK;
275 }
276 }
277
278 // Check for zero
279 if (!d_exp && !d_mant)
280 {
281 return (sign << 15);
282 }
283
284 // Check for overflow
285 if (exp >= CL_HALF_MAX_EXP)
286 {
287 return cl_half_handle_overflow(rounding_mode, sign);
288 }
289
290 // Check for underflow
291 if (exp < (CL_HALF_MIN_EXP - CL_HALF_MANT_DIG - 1))
292 {
293 return cl_half_handle_underflow(rounding_mode, sign);
294 }
295
296 // Check for value that will become denormal
297 if (exp < -14)
298 {
299 // Include the implicit 1 from the FP64 mantissa
300 h_exp = 0;
301 d_mant |= (uint64_t)1 << (CL_DBL_MANT_DIG - 1);
302
303 // Mantissa shift amount depends on exponent
304 lsb_pos = (uint32_t)(-exp + (CL_DBL_MANT_DIG - 25));
305 }
306
307 // Generate FP16 mantissa by shifting FP64 mantissa
308 uint16_t h_mant = (uint16_t)(d_mant >> lsb_pos);
309
310 // Check whether we need to round
311 uint64_t halfway = (uint64_t)1 << (lsb_pos - 1);
312 uint64_t mask = (halfway << 1) - 1;
313 switch (rounding_mode)
314 {
315 case CL_HALF_RTE:
316 if ((d_mant & mask) > halfway)
317 {
318 // More than halfway -> round up
319 h_mant += 1;
320 }
321 else if ((d_mant & mask) == halfway)
322 {
323 // Exactly halfway -> round to nearest even
324 if (h_mant & 0x1)
325 h_mant += 1;
326 }
327 break;
328 case CL_HALF_RTZ:
329 // Mantissa has already been truncated -> do nothing
330 break;
331 case CL_HALF_RTP:
332 if ((d_mant & mask) && !sign)
333 {
334 // Round positive numbers up
335 h_mant += 1;
336 }
337 break;
338 case CL_HALF_RTN:
339 if ((d_mant & mask) && sign)
340 {
341 // Round negative numbers down
342 h_mant += 1;
343 }
344 break;
345 }
346
347 // Check for mantissa overflow
348 if (h_mant & 0x400)
349 {
350 h_exp += 1;
351 h_mant = 0;
352 }
353
354 return (sign << 15) | (h_exp << 10) | h_mant;
355 }
356
357
358 /**
359 * Convert a cl_half to a cl_float.
360 */
cl_half_to_float(cl_half h)361 static inline cl_float cl_half_to_float(cl_half h)
362 {
363 // Type-punning to get direct access to underlying bits
364 union
365 {
366 cl_float f;
367 uint32_t i;
368 } f32;
369
370 // Extract sign bit
371 uint16_t sign = h >> 15;
372
373 // Extract FP16 exponent and mantissa
374 uint16_t h_exp = (h >> (CL_HALF_MANT_DIG - 1)) & 0x1F;
375 uint16_t h_mant = h & 0x3FF;
376
377 // Remove FP16 exponent bias
378 int32_t exp = h_exp - CL_HALF_MAX_EXP + 1;
379
380 // Add FP32 exponent bias
381 uint32_t f_exp = exp + CL_FLT_MAX_EXP - 1;
382
383 // Check for NaN / infinity
384 if (h_exp == 0x1F)
385 {
386 if (h_mant)
387 {
388 // NaN -> propagate mantissa and silence it
389 uint32_t f_mant = h_mant << (CL_FLT_MANT_DIG - CL_HALF_MANT_DIG);
390 f_mant |= 0x400000;
391 f32.i = (sign << 31) | 0x7F800000 | f_mant;
392 return f32.f;
393 }
394 else
395 {
396 // Infinity -> zero mantissa
397 f32.i = (sign << 31) | 0x7F800000;
398 return f32.f;
399 }
400 }
401
402 // Check for zero / denormal
403 if (h_exp == 0)
404 {
405 if (h_mant == 0)
406 {
407 // Zero -> zero exponent
408 f_exp = 0;
409 }
410 else
411 {
412 // Denormal -> normalize it
413 // - Shift mantissa to make most-significant 1 implicit
414 // - Adjust exponent accordingly
415 uint32_t shift = 0;
416 while ((h_mant & 0x400) == 0)
417 {
418 h_mant <<= 1;
419 shift++;
420 }
421 h_mant &= 0x3FF;
422 f_exp -= shift - 1;
423 }
424 }
425
426 f32.i = (sign << 31) | (f_exp << 23) | (h_mant << 13);
427 return f32.f;
428 }
429
430
431 #undef CL_HALF_EXP_MASK
432 #undef CL_HALF_MAX_FINITE_MAG
433
434
435 #ifdef __cplusplus
436 }
437 #endif
438
439
440 #endif /* OPENCL_CL_HALF_H */
441