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
36typedef 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 */
56uint16_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      /* NaN */
87      m = 1;
88      e = 31;
89   }
90   else {
91      /* regular number */
92      const int new_exp = flt_e - 127;
93      if (new_exp < -14) {
94         /* The float32 lies in the range (0.0, min_normal16) and is rounded
95          * to a nearby float16 value. The result will be either zero, subnormal,
96          * or normal.
97          */
98         e = 0;
99         m = _mesa_lroundevenf((1 << 24) * fabsf(fi.f));
100      }
101      else if (new_exp > 15) {
102         /* map this value to infinity */
103         /* m = 0; - already set */
104         e = 31;
105      }
106      else {
107         /* The float32 lies in the range
108          *   [min_normal16, max_normal16 + max_step16)
109          * and is rounded to a nearby float16 value. The result will be
110          * either normal or infinite.
111          */
112         e = new_exp + 15;
113         m = _mesa_lroundevenf(flt_m / (float) (1 << 13));
114      }
115   }
116
117   assert(0 <= m && m <= 1024);
118   if (m == 1024) {
119      /* The float32 was rounded upwards into the range of the next exponent,
120       * so bump the exponent. This correctly handles the case where f32
121       * should be rounded up to float16 infinity.
122       */
123      ++e;
124      m = 0;
125   }
126
127   result = (s << 15) | (e << 10) | m;
128   return result;
129}
130
131uint16_t
132_mesa_float_to_float16_rtz_slow(float val)
133{
134    return _mesa_float_to_half_rtz_slow(val);
135}
136
137/**
138 * Convert a 2-byte half float to a 4-byte float.
139 * Based on code from:
140 * http://www.opengl.org/discussion_boards/ubb/Forum3/HTML/008786.html
141 */
142float
143_mesa_half_to_float_slow(uint16_t val)
144{
145   union fi infnan;
146   union fi magic;
147   union fi f32;
148
149   infnan.ui = 0x8f << 23;
150   infnan.f = 65536.0f;
151   magic.ui  = 0xef << 23;
152
153   /* Exponent / Mantissa */
154   f32.ui = (val & 0x7fff) << 13;
155
156   /* Adjust */
157   f32.f *= magic.f;
158   /* XXX: The magic mul relies on denorms being available */
159
160   /* Inf / NaN */
161   if (f32.f >= infnan.f)
162      f32.ui |= 0xff << 23;
163
164   /* Sign */
165   f32.ui |= (uint32_t)(val & 0x8000) << 16;
166
167   return f32.f;
168}
169
170/**
171  * Convert 0.0 to 0x00, 1.0 to 0xff.
172  * Values outside the range [0.0, 1.0] will give undefined results.
173  */
174uint8_t _mesa_half_to_unorm8(uint16_t val)
175{
176   const int m = val & 0x3ff;
177   const int e = (val >> 10) & 0x1f;
178   ASSERTED const int s = (val >> 15) & 0x1;
179
180   /* v = round_to_nearest(1.mmmmmmmmmm * 2^(e-15) * 255)
181    *   = round_to_nearest((1.mmmmmmmmmm * 255) * 2^(e-15))
182    *   = round_to_nearest((1mmmmmmmmmm * 255) * 2^(e-25))
183    *   = round_to_zero((1mmmmmmmmmm * 255) * 2^(e-25) + 0.5)
184    *   = round_to_zero(((1mmmmmmmmmm * 255) * 2^(e-24) + 1) / 2)
185    *
186    * This happens to give the correct answer for zero/subnormals too
187    */
188   assert(s == 0 && val <= FP16_ONE); /* check 0 <= this <= 1 */
189   /* (implies e <= 15, which means the bit-shifts below are safe) */
190
191   uint32_t v = ((1 << 10) | m) * 255;
192   v = ((v >> (24 - e)) + 1) >> 1;
193   return v;
194}
195
196/**
197  * Takes a uint16_t, divides by 65536, converts the infinite-precision
198  * result to fp16 with round-to-zero. Used by the ASTC decoder.
199  */
200uint16_t _mesa_uint16_div_64k_to_half(uint16_t v)
201{
202   /* Zero or subnormal. Set the mantissa to (v << 8) and return. */
203   if (v < 4)
204      return v << 8;
205
206   /* Count the leading 0s in the uint16_t */
207#ifdef HAVE___BUILTIN_CLZ
208   int n = __builtin_clz(v) - 16;
209#else
210   int n = 16;
211   for (int i = 15; i >= 0; i--) {
212      if (v & (1 << i)) {
213         n = 15 - i;
214         break;
215      }
216   }
217#endif
218
219   /* Shift the mantissa up so bit 16 is the hidden 1 bit,
220    * mask it off, then shift back down to 10 bits
221    */
222   int m = ( ((uint32_t)v << (n + 1)) & 0xffff ) >> 6;
223
224   /*  (0{n} 1 X{15-n}) * 2^-16
225    * = 1.X * 2^(15-n-16)
226    * = 1.X * 2^(14-n - 15)
227    * which is the FP16 form with e = 14 - n
228    */
229   int e = 14 - n;
230
231   assert(e >= 1 && e <= 30);
232   assert(m >= 0 && m < 0x400);
233
234   return (e << 10) | m;
235}
236