Changeset 872 for trunk/src/half.cpp
 Timestamp:
 Aug 29, 2011, 2:07:54 AM (9 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/half.cpp
r869 r872 20 20 { 21 21 22 /* These macros implement a finite iterator useful to build lookup 23 * tables. For instance, S64(0) will call S1(x) for all values of x 24 * between 0 and 63. 25 * Due to the exponential behaviour of the calls, the stress on the 26 * compiler may be important. */ 27 #define S4(x) S1((x)), S1((x)+1), S1((x)+2), S1((x)+3) 28 #define S16(x) S4((x)), S4((x)+4), S4((x)+8), S4((x)+12) 29 #define S64(x) S16((x)), S16((x)+16), S16((x)+32), S16((x)+48) 30 #define S256(x) S64((x)), S64((x)+64), S64((x)+128), S64((x)+192) 31 #define S1024(x) S256((x)), S256((x)+256), S256((x)+512), S256((x)+768) 32 22 33 /* Lookup tablebased algorithm from “Fast Half Float Conversions” 23 34 * by Jeroen van der Zijp, November 2008. No rounding is performed, 24 35 * and some NaN values may be incorrectly converted to Inf. */ 25 half half::makefast(float f) 26 { 27 #define S4(x) S1(4*(x)), S1(4*(x)+1), S1(4*(x)+2), S1(4*(x)+3) 28 #define S16(x) S4(4*(x)), S4(4*(x)+1), S4(4*(x)+2), S4(4*(x)+3) 29 #define S64(x) S16(4*(x)), S16(4*(x)+1), S16(4*(x)+2), S16(4*(x)+3) 30 #define S256(x) S64(4*(x)), S64(4*(x)+1), S64(4*(x)+2), S64(4*(x)+3) 31 36 static inline uint16_t float_to_half_nobranch(uint32_t x) 37 { 32 38 static uint16_t const basetable[512] = 33 39 { … … 53 59 }; 54 60 55 union { float f; uint32_t x; } u = { f }; 56 57 uint16_t bits = basetable[(u.x >> 23) & 0x1ff]; 58 bits = (u.x & 0x007fffff) >> shifttable[(u.x >> 23) & 0x1ff]; 59 return makebits(bits); 61 uint16_t bits = basetable[(x >> 23) & 0x1ff]; 62 bits = (x & 0x007fffff) >> shifttable[(x >> 23) & 0x1ff]; 63 return bits; 60 64 } 61 65 … … 63 67 * used, eg. in Ogre), with the additional benefit of rounding, inspired 64 68 * by James Tursa’s halfprecision code. */ 65 half half::makeslow(float f) 66 { 67 union { float f; uint32_t x; } u = { f }; 68 69 uint16_t bits = (u.x >> 16) & 0x8000; /* Get the sign */ 70 uint16_t m = (u.x >> 12) & 0x07ff; /* Keep one extra bit for rounding */ 71 unsigned int e = (u.x >> 23) & 0xff; /* Using int is faster here */ 69 static inline uint16_t float_to_half_branch(uint32_t x) 70 { 71 uint16_t bits = (x >> 16) & 0x8000; /* Get the sign */ 72 uint16_t m = (x >> 12) & 0x07ff; /* Keep one extra bit for rounding */ 73 unsigned int e = (x >> 23) & 0xff; /* Using int is faster here */ 72 74 73 75 /* If zero, or denormal, or exponent underflows too much for a denormal, 74 76 * return signed zero. */ 75 77 if (e < 103) 76 return makebits(bits);78 return bits; 77 79 78 80 /* If NaN, return NaN. If Inf or exponent overflow, return Inf. */ … … 82 84 /* If exponent was 0xff and one mantissa bit was set, it means NaN, 83 85 * not Inf, so make sure we set one mantissa bit too. */ 84 bits = e == 255 && ( u.x & 0x007fffffu);85 return makebits(bits);86 bits = e == 255 && (x & 0x007fffffu); 87 return bits; 86 88 } 87 89 … … 93 95 * to 1, which is OK. */ 94 96 bits = (m >> (114  e)) + ((m >> (113  e)) & 1); 95 return makebits(bits);97 return bits; 96 98 } 97 99 … … 100 102 * the exponent, which is OK. */ 101 103 bits += m & 1; 102 return makebits(bits); 103 } 104 105 half::operator float() const 106 { 107 union { float f; uint32_t x; } u; 108 109 uint32_t s = (m_bits & 0x8000u) << 16; 110 111 if ((m_bits & 0x7fffu) == 0) 112 { 113 u.x = (uint32_t)m_bits << 16; 114 return u.f; 115 } 116 117 uint32_t e = m_bits & 0x7c00u; 118 uint32_t m = m_bits & 0x03ffu; 104 return bits; 105 } 106 107 static int const shifttable[32] = 108 { 109 23, 14, 22, 0, 0, 0, 21, 0, 0, 0, 0, 0, 0, 0, 20, 0, 110 15, 0, 0, 0, 0, 0, 0, 16, 0, 0, 0, 17, 0, 18, 19, 0, 111 }; 112 static uint32_t const shiftmagic = 0x07c4acddu; 113 114 /* Lookup tablebased algorithm from “Fast Half Float Conversions” 115 * by Jeroen van der Zijp, November 2008. Tables are generated using 116 * the C++ preprocessor, thanks to a branchless implementation also 117 * used in half_to_float_branch(). This code is actually almost always 118 * slower than the branching one. */ 119 static inline uint32_t half_to_float_nobranch(uint16_t x) 120 { 121 #define M3(i) ((i)  ((i) >> 1)) 122 #define M7(i) (M3(i)  (M3(i) >> 2)) 123 #define MF(i) (M7(i)  (M7(i) >> 4)) 124 #define MFF(i) (MF(i)  (MF(i) >> 8)) 125 #define E(i) shifttable[(unsigned int)(MFF(i) * shiftmagic) >> 27] 126 127 static uint32_t const mantissatable[2048] = 128 { 129 #define S1(i) (((i) == 0) ? 0 : ((125  E(i)) << 23) + ((i) << E(i))) 130 S1024(0), 131 #undef S1 132 #define S1(i) (0x38000000u + ((i) << 13)) 133 S1024(0), 134 #undef S1 135 }; 136 137 static uint32_t const exponenttable[64] = 138 { 139 #define S1(i) (((i) == 0) ? 0 : \ 140 ((i) < 31) ? ((i) << 23) : \ 141 ((i) == 31) ? 0x47800000u : \ 142 ((i) == 32) ? 0x80000000u : \ 143 ((i) < 63) ? (0x80000000u + (((i)  32) << 23)) : 0xc7800000) 144 S64(0), 145 #undef S1 146 }; 147 148 static int const offsettable[64] = 149 { 150 #define S1(i) (((i) == 0  (i) == 32) ? 0 : 1024) 151 S64(0), 152 #undef S1 153 }; 154 155 return mantissatable[offsettable[x >> 10] + (x & 0x3ff)] 156 + exponenttable[x >> 10]; 157 } 158 159 /* This algorithm is similar to the OpenEXR implementation, except it 160 * uses branchless code in the denormal path. */ 161 static inline uint32_t half_to_float_branch(uint16_t x) 162 { 163 uint32_t s = (x & 0x8000u) << 16; 164 165 if ((x & 0x7fffu) == 0) 166 return (uint32_t)x << 16; 167 168 uint32_t e = x & 0x7c00u; 169 uint32_t m = x & 0x03ffu; 119 170 120 171 if (e == 0) 121 172 { 122 static int const shifttable[32] =123 {124 10, 1, 9, 0, 0, 0, 8, 0, 0, 0, 0, 0, 0, 0, 7, 0,125 2, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 4, 0, 5, 6, 0,126 };127 128 173 uint32_t v = m  (m >> 1); 129 174 v = v >> 2; … … 131 176 v = v >> 8; 132 177 133 e = shifttable[(v * 0x07C4ACDDU) >> 27]; 134 m <<= e; 178 e = shifttable[(v * shiftmagic) >> 27]; 135 179 136 180 /* We don't have to remove the 10th mantissa bit because it gets 137 181 * added to our underestimated exponent. */ 138 u.x = s  (((112  e) << 23) + (m << 13)); 139 return u.f; 182 return s  (((125  e) << 23) + (m << e)); 140 183 } 141 184 … … 145 188 * or any other trick I could find. sam */ 146 189 if (m == 0) 147 u.x = s  0x7f800000u; 148 else 149 u.x = s  0x7fc00000u; 150 151 return u.f; 152 } 153 154 u.x = s  (((e >> 10) + 112) << 23)  (m << 13); 155 190 return s  0x7f800000u; 191 return s  0x7fc00000u; 192 } 193 194 return s  (((e >> 10) + 112) << 23)  (m << 13); 195 } 196 197 half half::makefast(float f) 198 { 199 union { float f; uint32_t x; } u = { f }; 200 return makebits(float_to_half_nobranch(u.x)); 201 } 202 203 half half::makeslow(float f) 204 { 205 union { float f; uint32_t x; } u = { f }; 206 return makebits(float_to_half_branch(u.x)); 207 } 208 209 half::operator float() const 210 { 211 union { float f; uint32_t x; } u; 212 u.x = half_to_float_branch(bits); 156 213 return u.f; 157 214 }
Note: See TracChangeset
for help on using the changeset viewer.