1 /******************************************************************************
2 *
3 * Filename: ieeehalfprecision.c
4 * Programmer: James Tursa
5 * Version: 1.0
6 * Date: March 3, 2009
7 * Copyright: (c) 2009 by James Tursa, All Rights Reserved
8 *
9 * This code uses the BSD License:
10 *
11 * Redistribution and use in source and binary forms, with or without
12 * modification, are permitted provided that the following conditions are
13 * met:
14 *
15 * * Redistributions of source code must retain the above copyright
16 * notice, this list of conditions and the following disclaimer.
17 * * Redistributions in binary form must reproduce the above copyright
18 * notice, this list of conditions and the following disclaimer in
19 * the documentation and/or other materials provided with the distribution
20 *
21 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
22 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
24 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
25 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
26 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
27 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
28 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
29 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
30 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31 * POSSIBILITY OF SUCH DAMAGE.
32 *
33 * This file contains C code to convert between IEEE double, single, and half
34 * precision floating point formats. The intended use is for standalone C code
35 * that does not rely on MATLAB mex.h. The bit pattern for the half precision
36 * floating point format is stored in a 16-bit unsigned int variable. The half
37 * precision bit pattern definition is:
38 *
39 * 1 bit sign bit
40 * 5 bits exponent, biased by 15
41 * 10 bits mantissa, hidden leading bit, normalized to 1.0
42 *
43 * Special floating point bit patterns recognized and supported:
44 *
45 * All exponent bits zero:
46 * - If all mantissa bits are zero, then number is zero (possibly signed)
47 * - Otherwise, number is a denormalized bit pattern
48 *
49 * All exponent bits set to 1:
50 * - If all mantissa bits are zero, then number is +Infinity or -Infinity
51 * - Otherwise, number is NaN (Not a Number)
52 *
53 * For the denormalized cases, note that 2^(-24) is the smallest number that can
54 * be represented in half precision exactly. 2^(-25) will convert to 2^(-24)
55 * because of the rounding algorithm used, and 2^(-26) is too small and
56 * underflows to zero.
57 *
58 ******************************************************************************/
59
60 /*
61 changes by K. Rogovin:
62 - changed macros UINT16_TYPE, etc to types from stdint.h
63 (i.e. UINT16_TYPE-->uint16_t, INT16_TYPE-->int16_t, etc)
64
65 - removed double conversion routines.
66
67 - changed run time checks of endianness to compile time macro.
68
69 - removed return value from routines
70
71 - changed source parameter type from * to const *
72
73 - changed pointer types from void ot uint16_t and uint32_t
74 */
75
76 /*
77 * andy@warmcat.com:
78 *
79 * - clean style and indenting
80 * - convert to single operation
81 * - export as lws_
82 */
83
84 #include <string.h>
85 #include <stdint.h>
86
87 void
lws_singles2halfp(uint16_t * hp,uint32_t x)88 lws_singles2halfp(uint16_t *hp, uint32_t x)
89 {
90 uint32_t xs, xe, xm;
91 uint16_t hs, he, hm;
92 int hes;
93
94 if (!(x & 0x7FFFFFFFu)) {
95 /* Signed zero */
96 *hp = (uint16_t)(x >> 16);
97
98 return;
99 }
100
101 xs = x & 0x80000000u; // Pick off sign bit
102 xe = x & 0x7F800000u; // Pick off exponent bits
103 xm = x & 0x007FFFFFu; // Pick off mantissa bits
104
105 if (xe == 0) { // Denormal will underflow, return a signed zero
106 *hp = (uint16_t) (xs >> 16);
107 return;
108 }
109
110 if (xe == 0x7F800000u) { // Inf or NaN (all the exponent bits are set)
111 if (!xm) { // If mantissa is zero ...
112 *hp = (uint16_t) ((xs >> 16) | 0x7C00u); // Signed Inf
113 return;
114 }
115
116 *hp = (uint16_t) 0xFE00u; // NaN, only 1st mantissa bit set
117
118 return;
119 }
120
121 /* Normalized number */
122
123 hs = (uint16_t) (xs >> 16); // Sign bit
124 /* Exponent unbias the single, then bias the halfp */
125 hes = ((int)(xe >> 23)) - 127 + 15;
126
127 if (hes >= 0x1F) { // Overflow
128 *hp = (uint16_t) ((xs >> 16) | 0x7C00u); // Signed Inf
129 return;
130 }
131
132 if (hes <= 0) { // Underflow
133 if ((14 - hes) > 24)
134 /*
135 * Mantissa shifted all the way off & no
136 * rounding possibility
137 */
138 hm = (uint16_t) 0u; // Set mantissa to zero
139 else {
140 xm |= 0x00800000u; // Add the hidden leading bit
141 hm = (uint16_t) (xm >> (14 - hes)); // Mantissa
142 if ((xm >> (13 - hes)) & 1u) // Check for rounding
143 /* Round, might overflow into exp bit,
144 * but this is OK */
145 hm = (uint16_t)(hm + 1u);
146 }
147 /* Combine sign bit and mantissa bits, biased exponent is 0 */
148 *hp = hs | hm;
149
150 return;
151 }
152
153 he = (uint16_t)(hes << 10); // Exponent
154 hm = (uint16_t)(xm >> 13); // Mantissa
155
156 if (xm & 0x00001000u) // Check for rounding
157 /* Round, might overflow to inf, this is OK */
158 *hp = (uint16_t)((hs | he | hm) + (uint16_t)1u);
159 else
160 *hp = hs | he | hm; // No rounding
161 }
162
163 void
lws_halfp2singles(uint32_t * xp,uint16_t h)164 lws_halfp2singles(uint32_t *xp, uint16_t h)
165 {
166 uint16_t hs, he, hm;
167 uint32_t xs, xe, xm;
168 int32_t xes;
169 int e;
170
171 if (!(h & 0x7FFFu)) { // Signed zero
172 *xp = ((uint32_t)h) << 16; // Return the signed zero
173
174 return;
175 }
176
177 hs = h & 0x8000u; // Pick off sign bit
178 he = h & 0x7C00u; // Pick off exponent bits
179 hm = h & 0x03FFu; // Pick off mantissa bits
180
181 if (!he) { // Denormal will convert to normalized
182 e = -1;
183
184 /* figure out how much extra to adjust the exponent */
185 do {
186 e++;
187 hm = (uint16_t)(hm << 1);
188 /* Shift until leading bit overflows into exponent */
189 } while (!(hm & 0x0400u));
190
191 xs = ((uint32_t) hs) << 16; // Sign bit
192
193 /* Exponent unbias the halfp, then bias the single */
194 xes = ((int32_t)(he >> 10)) - 15 + 127 - e;
195 xe = (uint32_t)(xes << 23); // Exponent
196 xm = ((uint32_t)(hm & 0x03FFu)) << 13; // Mantissa
197
198 *xp = xs | xe | xm;
199
200 return;
201 }
202
203 if (he == 0x7C00u) { /* Inf or NaN (all the exponent bits are set) */
204 if (!hm) { /* If mantissa is zero ...
205 * Signed Inf
206 */
207 *xp = (((uint32_t)hs) << 16) | ((uint32_t)0x7F800000u);
208
209 return;
210 }
211
212 /* ... NaN, only 1st mantissa bit set */
213 *xp = (uint32_t)0xFFC00000u;
214
215 return;
216 }
217
218 /* Normalized number */
219
220 xs = ((uint32_t)hs) << 16; // Sign bit
221 /* Exponent unbias the halfp, then bias the single */
222 xes = ((int32_t)(he >> 10)) - 15 + 127;
223 xe = (uint32_t)(xes << 23); // Exponent
224 xm = ((uint32_t)hm) << 13; // Mantissa
225
226 /* Combine sign bit, exponent bits, and mantissa bits */
227 *xp = xs | xe | xm;
228 }
229