Sheri Zhang | 79cb945 | 2021-09-07 14:51:49 +0100 | [diff] [blame^] | 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 | */ |
| 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 | */ |
| 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 | */ |
| 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 | */ |
| 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 | */ |
| 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 */ |