• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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