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
  32 typedef union { float f; int32_t i; uint32_t u; } fi_type;
  33
  34 /**
  35  * Convert a 4-byte float to a 2-byte half float.
  36  *
  37  * Not all float32 values can be represented exactly as a float16 value. We
  38  * round such intermediate float32 values to the nearest float16. When the
  39  * float32 lies exactly between to float16 values, we round to the one with
  40  * an even mantissa.
  41  *
  42  * This rounding behavior has several benefits:
  43  *   - It has no sign bias.
  44  *
  45  *   - It reproduces the behavior of real hardware: opcode F32TO16 in Intel's
  46  *     GPU ISA.
  47  *
  48  *   - By reproducing the behavior of the GPU (at least on Intel hardware),
  49  *     compile-time evaluation of constant packHalf2x16 GLSL expressions will
  50  *     result in the same value as if the expression were executed on the GPU.
  51  */
  52 uint16_t
  53 _mesa_float_to_half(float val)
  54 {
  55    const fi_type fi = {val};
  56    const int flt_m = fi.i & 0x7fffff;
  57    const int flt_e = (fi.i >> 23) & 0xff;
  58    const int flt_s = (fi.i >> 31) & 0x1;
  59    int s, e, m = 0;
  60    uint16_t result;
  61
  62    /* sign bit */
  63    s = flt_s;
  64
  65    /* handle special cases */
  66    if ((flt_e == 0) && (flt_m == 0)) {
  67       /* zero */
  68       /* m = 0; - already set */
  69       e = 0;
  70    }
  71    else if ((flt_e == 0) && (flt_m != 0)) {
  72       /* denorm -- denorm float maps to 0 half */
  73       /* m = 0; - already set */
  74       e = 0;
  75    }
  76    else if ((flt_e == 0xff) && (flt_m == 0)) {
  77       /* infinity */
  78       /* m = 0; - already set */
  79       e = 31;
  80    }
  81    else if ((flt_e == 0xff) && (flt_m != 0)) {
  82       /* NaN */
  83       m = 1;
  84       e = 31;
  85    }
  86    else {
  87       /* regular number */
  88       const int new_exp = flt_e - 127;
  89       if (new_exp < -14) {
  90          /* The float32 lies in the range (0.0, min_normal16) and is rounded
  91           * to a nearby float16 value. The result will be either zero, subnormal,
  92           * or normal.
  93           */
  94          e = 0;
  95          m = _mesa_lroundevenf((1 << 24) * fabsf(fi.f));
  96       }
  97       else if (new_exp > 15) {
  98          /* map this value to infinity */
  99          /* m = 0; - already set */
 100          e = 31;
 101       }
 102       else {
 103          /* The float32 lies in the range
 104           *   [min_normal16, max_normal16 + max_step16)
 105           * and is rounded to a nearby float16 value. The result will be
 106           * either normal or infinite.
 107           */
 108          e = new_exp + 15;
 109          m = _mesa_lroundevenf(flt_m / (float) (1 << 13));
 110       }
 111    }
 112
 113    assert(0 <= m && m <= 1024);
 114    if (m == 1024) {
 115       /* The float32 was rounded upwards into the range of the next exponent,
 116        * so bump the exponent. This correctly handles the case where f32
 117        * should be rounded up to float16 infinity.
 118        */
 119       ++e;
 120       m = 0;
 121    }
 122
 123    result = (s << 15) | (e << 10) | m;
 124    return result;
 125 }
 126
 127
 128 /**
 129  * Convert a 2-byte half float to a 4-byte float.
 130  * Based on code from:
 131  * http://www.opengl.org/discussion_boards/ubb/Forum3/HTML/008786.html
 132  */
 133 float
 134 _mesa_half_to_float(uint16_t val)
 135 {
 136    /* XXX could also use a 64K-entry lookup table */
 137    const int m = val & 0x3ff;
 138    const int e = (val >> 10) & 0x1f;
 139    const int s = (val >> 15) & 0x1;
 140    int flt_m, flt_e, flt_s;
 141    fi_type fi;
 142    float result;
 143
 144    /* sign bit */
 145    flt_s = s;
 146
 147    /* handle special cases */
 148    if ((e == 0) && (m == 0)) {
 149       /* zero */
 150       flt_m = 0;
 151       flt_e = 0;
 152    }
 153    else if ((e == 0) && (m != 0)) {
 154       /* denorm -- denorm half will fit in non-denorm single */
 155       const float half_denorm = 1.0f / 16384.0f; /* 2^-14 */
 156       float mantissa = ((float) (m)) / 1024.0f;
 157       float sign = s ? -1.0f : 1.0f;
 158       return sign * mantissa * half_denorm;
 159    }
 160    else if ((e == 31) && (m == 0)) {
 161       /* infinity */
 162       flt_e = 0xff;
 163       flt_m = 0;
 164    }
 165    else if ((e == 31) && (m != 0)) {
 166       /* NaN */
 167       flt_e = 0xff;
 168       flt_m = 1;
 169    }
 170    else {
 171       /* regular */
 172       flt_e = e + 112;
 173       flt_m = m << 13;
 174    }
 175
 176    fi.i = (flt_s << 31) | (flt_e << 23) | flt_m;
 177    result = fi.f;
 178    return result;
 179 }
 180
 181 /**
 182   * Convert 0.0 to 0x00, 1.0 to 0xff.
 183   * Values outside the range [0.0, 1.0] will give undefined results.
 184   */
 185 uint8_t _mesa_half_to_unorm8(uint16_t val)
 186 {
 187    const int m = val & 0x3ff;
 188    const int e = (val >> 10) & 0x1f;
 189    const int s = (val >> 15) & 0x1;
 190
 191    /* v = round_to_nearest(1.mmmmmmmmmm * 2^(e-15) * 255)
 192     *   = round_to_nearest((1.mmmmmmmmmm * 255) * 2^(e-15))
 193     *   = round_to_nearest((1mmmmmmmmmm * 255) * 2^(e-25))
 194     *   = round_to_zero((1mmmmmmmmmm * 255) * 2^(e-25) + 0.5)
 195     *   = round_to_zero(((1mmmmmmmmmm * 255) * 2^(e-24) + 1) / 2)
 196     *
 197     * This happens to give the correct answer for zero/subnormals too
 198     */
 199    assert(s == 0 && val <= FP16_ONE); /* check 0 <= this <= 1 */
 200    /* (implies e <= 15, which means the bit-shifts below are safe) */
 201
 202    uint32_t v = ((1 << 10) | m) * 255;
 203    v = ((v >> (24 - e)) + 1) >> 1;
 204    return v;
 205 }
 206
 207 /**
 208   * Takes a uint16_t, divides by 65536, converts the infinite-precision
 209   * result to fp16 with round-to-zero. Used by the ASTC decoder.
 210   */
 211 uint16_t _mesa_uint16_div_64k_to_half(uint16_t v)
 212 {
 213    /* Zero or subnormal. Set the mantissa to (v << 8) and return. */
 214    if (v < 4)
 215       return v << 8;
 216
 217    /* Count the leading 0s in the uint16_t */
 218 #ifdef HAVE___BUILTIN_CLZ
 219    int n = __builtin_clz(v) - 16;
 220 #else
 221    int n = 16;
 222    for (int i = 15; i >= 0; i--) {
 223       if (v & (1 << i)) {
 224          n = 15 - i;
 225          break;
 226       }
 227    }
 228 #endif
 229
 230    /* Shift the mantissa up so bit 16 is the hidden 1 bit,
 231     * mask it off, then shift back down to 10 bits
 232     */
 233    int m = ( ((uint32_t)v << (n + 1)) & 0xffff ) >> 6;
 234
 235    /*  (0{n} 1 X{15-n}) * 2^-16
 236     * = 1.X * 2^(15-n-16)
 237     * = 1.X * 2^(14-n - 15)
 238     * which is the FP16 form with e = 14 - n
 239     */
 240    int e = 14 - n;
 241
 242    assert(e >= 1 && e <= 30);
 243    assert(m >= 0 && m < 0x400);
 244
 245    return (e << 10) | m;
 246 }