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
33typedef 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 */
53uint16_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 */
134float
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  */
186uint8_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  */
212uint16_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}
248