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 87void 88lws_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 163void 164lws_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