xref: /third_party/opencl-headers/CL/cl_half.h (revision 80d59932)
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
37extern "C" {
38#endif
39
40
41/**
42 * Rounding mode used when converting to cl_half.
43 */
44typedef 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 */
61static 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 */
87static 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 */
109static 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 */
235static 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 */
361static 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