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  * Copyright (C) 2018-2019 Intel Corporation
   8  *
   9  * Permission is hereby granted, free of charge, to any person obtaining a
  10  * copy of this software and associated documentation files (the "Software"),
  11  * to deal in the Software without restriction, including without limitation
  12  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
  13  * and/or sell copies of the Software, and to permit persons to whom the
  14  * Software is furnished to do so, subject to the following conditions:
  15  *
  16  * The above copyright notice and this permission notice shall be included
  17  * in all copies or substantial portions of the Software.
  18  *
  19  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  20  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  21  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
  22  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
  23  * OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
  24  * ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
  25  * OTHER DEALINGS IN THE SOFTWARE.
  26  */
  27
  28 #include <math.h>
  29 #include <assert.h>
  30 #include "half_float.h"
  31 #include "rounding.h"
  32 #include "softfloat.h"
  33 #include "macros.h"
  34 #include "u_math.h"
  35
  36 typedef union { float f; int32_t i; uint32_t u; } fi_type;
  37
  38 /**
  39  * Convert a 4-byte float to a 2-byte half float.
  40  *
  41  * Not all float32 values can be represented exactly as a float16 value. We
  42  * round such intermediate float32 values to the nearest float16. When the
  43  * float32 lies exactly between to float16 values, we round to the one with
  44  * an even mantissa.
  45  *
  46  * This rounding behavior has several benefits:
  47  *   - It has no sign bias.
  48  *
  49  *   - It reproduces the behavior of real hardware: opcode F32TO16 in Intel's
  50  *     GPU ISA.
  51  *
  52  *   - By reproducing the behavior of the GPU (at least on Intel hardware),
  53  *     compile-time evaluation of constant packHalf2x16 GLSL expressions will
  54  *     result in the same value as if the expression were executed on the GPU.
  55  */
  56 uint16_t
  57 _mesa_float_to_half_slow(float val)
  58 {
  59    const fi_type fi = {val};
  60    const int flt_m = fi.i & 0x7fffff;
  61    const int flt_e = (fi.i >> 23) & 0xff;
  62    const int flt_s = (fi.i >> 31) & 0x1;
  63    int s, e, m = 0;
  64    uint16_t result;
  65
  66    /* sign bit */
  67    s = flt_s;
  68
  69    /* handle special cases */
  70    if ((flt_e == 0) && (flt_m == 0)) {
  71       /* zero */
  72       /* m = 0; - already set */
  73       e = 0;
  74    }
  75    else if ((flt_e == 0) && (flt_m != 0)) {
  76       /* denorm -- denorm float maps to 0 half */
  77       /* m = 0; - already set */
  78       e = 0;
  79    }
  80    else if ((flt_e == 0xff) && (flt_m == 0)) {
  81       /* infinity */
  82       /* m = 0; - already set */
  83       e = 31;
  84    }
  85    else if ((flt_e == 0xff) && (flt_m != 0)) {
  86       /* Retain the top bits of a NaN to make sure that the quiet/signaling
  87        * status stays the same.
  88        */
  89       m = flt_m >> 13;
  90       if (!m)
  91          m = 1;
  92       e = 31;
  93    }
  94    else {
  95       /* regular number */
  96       const int new_exp = flt_e - 127;
  97       if (new_exp < -14) {
  98          /* The float32 lies in the range (0.0, min_normal16) and is rounded
  99           * to a nearby float16 value. The result will be either zero, subnormal,
 100           * or normal.
 101           */
 102          e = 0;
 103          m = _mesa_lroundevenf((1 << 24) * fabsf(fi.f));
 104       }
 105       else if (new_exp > 15) {
 106          /* map this value to infinity */
 107          /* m = 0; - already set */
 108          e = 31;
 109       }
 110       else {
 111          /* The float32 lies in the range
 112           *   [min_normal16, max_normal16 + max_step16)
 113           * and is rounded to a nearby float16 value. The result will be
 114           * either normal or infinite.
 115           */
 116          e = new_exp + 15;
 117          m = _mesa_lroundevenf(flt_m / (float) (1 << 13));
 118       }
 119    }
 120
 121    assert(0 <= m && m <= 1024);
 122    if (m == 1024) {
 123       /* The float32 was rounded upwards into the range of the next exponent,
 124        * so bump the exponent. This correctly handles the case where f32
 125        * should be rounded up to float16 infinity.
 126        */
 127       ++e;
 128       m = 0;
 129    }
 130
 131    result = (s << 15) | (e << 10) | m;
 132    return result;
 133 }
 134
 135 uint16_t
 136 _mesa_float_to_float16_rtz_slow(float val)
 137 {
 138     return _mesa_float_to_half_rtz_slow(val);
 139 }
 140
 141 /**
 142  * Convert a 2-byte half float to a 4-byte float.
 143  * Based on code from:
 144  * http://www.opengl.org/discussion_boards/ubb/Forum3/HTML/008786.html
 145  */
 146 float
 147 _mesa_half_to_float_slow(uint16_t val)
 148 {
 149    union fi infnan;
 150    union fi magic;
 151    union fi f32;
 152
 153    infnan.ui = 0x8f << 23;
 154    infnan.f = 65536.0f;
 155    magic.ui  = 0xef << 23;
 156
 157    /* Exponent / Mantissa */
 158    f32.ui = (val & 0x7fff) << 13;
 159
 160    /* Adjust */
 161    f32.f *= magic.f;
 162    /* XXX: The magic mul relies on denorms being available */
 163
 164    /* Inf / NaN */
 165    if (f32.f >= infnan.f)
 166       f32.ui |= 0xff << 23;
 167
 168    /* Sign */
 169    f32.ui |= (uint32_t)(val & 0x8000) << 16;
 170
 171    return f32.f;
 172 }
 173
 174 /**
 175   * Takes a uint16_t, divides by 65536, converts the infinite-precision
 176   * result to fp16 with round-to-zero. Used by the ASTC decoder.
 177   */
 178 uint16_t _mesa_uint16_div_64k_to_half(uint16_t v)
 179 {
 180    /* Zero or subnormal. Set the mantissa to (v << 8) and return. */
 181    if (v < 4)
 182       return v << 8;
 183
 184    /* Count the leading 0s in the uint16_t */
 185 #ifdef HAVE___BUILTIN_CLZ
 186    int n = __builtin_clz(v) - 16;
 187 #else
 188    int n = 16;
 189    for (int i = 15; i >= 0; i--) {
 190       if (v & (1 << i)) {
 191          n = 15 - i;
 192          break;
 193       }
 194    }
 195 #endif
 196
 197    /* Shift the mantissa up so bit 16 is the hidden 1 bit,
 198     * mask it off, then shift back down to 10 bits
 199     */
 200    int m = ( ((uint32_t)v << (n + 1)) & 0xffff ) >> 6;
 201
 202    /*  (0{n} 1 X{15-n}) * 2^-16
 203     * = 1.X * 2^(15-n-16)
 204     * = 1.X * 2^(14-n - 15)
 205     * which is the FP16 form with e = 14 - n
 206     */
 207    int e = 14 - n;
 208
 209    assert(e >= 1 && e <= 30);
 210    assert(m >= 0 && m < 0x400);
 211
 212    return (e << 10) | m;
 213 }