src/util/half_float.c

   1 /*
   2  * Mesa 3-D graphics library
   3  *
   4  * Copyright (C) 1999-2007  Brian Paul   All Rights Reserved.
   5  * Copyright 2015 Philip Taylor <philip@zaynar.co.uk>
   6  * Copyright 2018 Advanced Micro Devices, Inc.
   7  *
   8  * Permission is hereby granted, free of charge, to any person obtaining a
   9  * copy of this software and associated documentation files (the "Software"),
  10  * to deal in the Software without restriction, including without limitation
  11  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
  12  * and/or sell copies of the Software, and to permit persons to whom the
  13  * Software is furnished to do so, subject to the following conditions:
  14  *
  15  * The above copyright notice and this permission notice shall be included
  16  * in all copies or substantial portions of the Software.
  17  *
  18  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  19  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  20  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
  21  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
  22  * OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
  23  * ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
  24  * OTHER DEALINGS IN THE SOFTWARE.
  25  */
  26
  27 #include <math.h>
  28 #include <assert.h>
  29 #include "half_float.h"
  30 #include "rounding.h"
  31 #include "macros.h"
  32
  33 typedef union { float f; int32_t i; uint32_t u; } fi_type;
  34
  35 /**
  36  * Convert a 4-byte float to a 2-byte half float.
  37  *
  38  * Not all float32 values can be represented exactly as a float16 value. We
  39  * round such intermediate float32 values to the nearest float16. When the
  40  * float32 lies exactly between to float16 values, we round to the one with
  41  * an even mantissa.
  42  *
  43  * This rounding behavior has several benefits:
  44  *   - It has no sign bias.
  45  *
  46  *   - It reproduces the behavior of real hardware: opcode F32TO16 in Intel's
  47  *     GPU ISA.
  48  *
  49  *   - By reproducing the behavior of the GPU (at least on Intel hardware),
  50  *     compile-time evaluation of constant packHalf2x16 GLSL expressions will
  51  *     result in the same value as if the expression were executed on the GPU.
  52  */
  53 uint16_t
  54 _mesa_float_to_half(float val)
  55 {
  56    const fi_type fi = {val};
  57    const int flt_m = fi.i & 0x7fffff;
  58    const int flt_e = (fi.i >> 23) & 0xff;
  59    const int flt_s = (fi.i >> 31) & 0x1;
  60    int s, e, m = 0;
  61    uint16_t result;
  62
  63    /* sign bit */
  64    s = flt_s;
  65
  66    /* handle special cases */
  67    if ((flt_e == 0) && (flt_m == 0)) {
  68       /* zero */
  69       /* m = 0; - already set */
  70       e = 0;
  71    }
  72    else if ((flt_e == 0) && (flt_m != 0)) {
  73       /* denorm -- denorm float maps to 0 half */
  74       /* m = 0; - already set */
  75       e = 0;
  76    }
  77    else if ((flt_e == 0xff) && (flt_m == 0)) {
  78       /* infinity */
  79       /* m = 0; - already set */
  80       e = 31;
  81    }
  82    else if ((flt_e == 0xff) && (flt_m != 0)) {
  83       /* NaN */
  84       m = 1;
  85       e = 31;
  86    }
  87    else {
  88       /* regular number */
  89       const int new_exp = flt_e - 127;
  90       if (new_exp < -14) {
  91          /* The float32 lies in the range (0.0, min_normal16) and is rounded
  92           * to a nearby float16 value. The result will be either zero, subnormal,
  93           * or normal.
  94           */
  95          e = 0;
  96          m = _mesa_lroundevenf((1 << 24) * fabsf(fi.f));
  97       }
  98       else if (new_exp > 15) {
  99          /* map this value to infinity */
 100          /* m = 0; - already set */
 101          e = 31;
 102       }
 103       else {
 104          /* The float32 lies in the range
 105           *   [min_normal16, max_normal16 + max_step16)
 106           * and is rounded to a nearby float16 value. The result will be
 107           * either normal or infinite.
 108           */
 109          e = new_exp + 15;
 110          m = _mesa_lroundevenf(flt_m / (float) (1 << 13));
 111       }
 112    }
 113
 114    assert(0 <= m && m <= 1024);
 115    if (m == 1024) {
 116       /* The float32 was rounded upwards into the range of the next exponent,
 117        * so bump the exponent. This correctly handles the case where f32
 118        * should be rounded up to float16 infinity.
 119        */
 120       ++e;
 121       m = 0;
 122    }
 123
 124    result = (s << 15) | (e << 10) | m;
 125    return result;
 126 }
 127
 128
 129 /**
 130  * Convert a 2-byte half float to a 4-byte float.
 131  * Based on code from:
 132  * http://www.opengl.org/discussion_boards/ubb/Forum3/HTML/008786.html
 133  */
 134 float
 135 _mesa_half_to_float(uint16_t val)
 136 {
 137    /* XXX could also use a 64K-entry lookup table */
 138    const int m = val & 0x3ff;
 139    const int e = (val >> 10) & 0x1f;
 140    const int s = (val >> 15) & 0x1;
 141    int flt_m, flt_e, flt_s;
 142    fi_type fi;
 143    float result;
 144
 145    /* sign bit */
 146    flt_s = s;
 147
 148    /* handle special cases */
 149    if ((e == 0) && (m == 0)) {
 150       /* zero */
 151       flt_m = 0;
 152       flt_e = 0;
 153    }
 154    else if ((e == 0) && (m != 0)) {
 155       /* denorm -- denorm half will fit in non-denorm single */
 156       const float half_denorm = 1.0f / 16384.0f; /* 2^-14 */
 157       float mantissa = ((float) (m)) / 1024.0f;
 158       float sign = s ? -1.0f : 1.0f;
 159       return sign * mantissa * half_denorm;
 160    }
 161    else if ((e == 31) && (m == 0)) {
 162       /* infinity */
 163       flt_e = 0xff;
 164       flt_m = 0;
 165    }
 166    else if ((e == 31) && (m != 0)) {
 167       /* NaN */
 168       flt_e = 0xff;
 169       flt_m = 1;
 170    }
 171    else {
 172       /* regular */
 173       flt_e = e + 112;
 174       flt_m = m << 13;
 175    }
 176
 177    fi.i = (flt_s << 31) | (flt_e << 23) | flt_m;
 178    result = fi.f;
 179    return result;
 180 }
 181
 182 /**
 183   * Convert 0.0 to 0x00, 1.0 to 0xff.
 184   * Values outside the range [0.0, 1.0] will give undefined results.
 185   */
 186 uint8_t _mesa_half_to_unorm8(uint16_t val)
 187 {
 188    const int m = val & 0x3ff;
 189    const int e = (val >> 10) & 0x1f;
 190    MAYBE_UNUSED const int s = (val >> 15) & 0x1;
 191
 192    /* v = round_to_nearest(1.mmmmmmmmmm * 2^(e-15) * 255)
 193     *   = round_to_nearest((1.mmmmmmmmmm * 255) * 2^(e-15))
 194     *   = round_to_nearest((1mmmmmmmmmm * 255) * 2^(e-25))
 195     *   = round_to_zero((1mmmmmmmmmm * 255) * 2^(e-25) + 0.5)
 196     *   = round_to_zero(((1mmmmmmmmmm * 255) * 2^(e-24) + 1) / 2)
 197     *
 198     * This happens to give the correct answer for zero/subnormals too
 199     */
 200    assert(s == 0 && val <= FP16_ONE); /* check 0 <= this <= 1 */
 201    /* (implies e <= 15, which means the bit-shifts below are safe) */
 202
 203    uint32_t v = ((1 << 10) | m) * 255;
 204    v = ((v >> (24 - e)) + 1) >> 1;
 205    return v;
 206 }
 207
 208 /**
 209   * Takes a uint16_t, divides by 65536, converts the infinite-precision
 210   * result to fp16 with round-to-zero. Used by the ASTC decoder.
 211   */
 212 uint16_t _mesa_uint16_div_64k_to_half(uint16_t v)
 213 {
 214    /* Zero or subnormal. Set the mantissa to (v << 8) and return. */
 215    if (v < 4)
 216       return v << 8;
 217
 218    /* Count the leading 0s in the uint16_t */
 219 #ifdef HAVE___BUILTIN_CLZ
 220    int n = __builtin_clz(v) - 16;
 221 #else
 222    int n = 16;
 223    for (int i = 15; i >= 0; i--) {
 224       if (v & (1 << i)) {
 225          n = 15 - i;
 226          break;
 227       }
 228    }
 229 #endif
 230
 231    /* Shift the mantissa up so bit 16 is the hidden 1 bit,
 232     * mask it off, then shift back down to 10 bits
 233     */
 234    int m = ( ((uint32_t)v << (n + 1)) & 0xffff ) >> 6;
 235
 236    /*  (0{n} 1 X{15-n}) * 2^-16
 237     * = 1.X * 2^(15-n-16)
 238     * = 1.X * 2^(14-n - 15)
 239     * which is the FP16 form with e = 14 - n
 240     */
 241    int e = 14 - n;
 242
 243    assert(e >= 1 && e <= 30);
 244    assert(m >= 0 && m < 0x400);
 245
 246    return (e << 10) | m;
 247 }