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