source: trunk/src/math/half.cpp @ 2183

Last change on this file since 2183 was 2183, checked in by sam, 7 years ago

build: fix the WTFPL site URL in all code comments.

  • Property svn:keywords set to Id
File size: 8.1 KB
Line 
1//
2// Lol Engine
3//
4// Copyright: (c) 2010-2012 Sam Hocevar <sam@hocevar.net>
5//   This program is free software; you can redistribute it and/or
6//   modify it under the terms of the Do What The Fuck You Want To
7//   Public License, Version 2, as published by Sam Hocevar. See
8//   http://www.wtfpl.net/ for more details.
9//
10
11#if defined HAVE_CONFIG_H
12#   include "config.h"
13#endif
14
15#if defined __CELLOS_LV2__
16#   if defined __SNC__
17#       include <ppu_altivec_internals.h>
18#   else
19#       include <altivec.h>
20#   endif
21#endif
22
23#include "core.h"
24
25using namespace std;
26
27namespace lol
28{
29
30/* These macros implement a finite iterator useful to build lookup
31 * tables. For instance, S64(0) will call S1(x) for all values of x
32 * between 0 and 63.
33 * Due to the exponential behaviour of the calls, the stress on the
34 * compiler may be important. */
35#define S4(x)    S1((x)),   S1((x)+1),     S1((x)+2),     S1((x)+3)
36#define S16(x)   S4((x)),   S4((x)+4),     S4((x)+8),     S4((x)+12)
37#define S64(x)   S16((x)),  S16((x)+16),   S16((x)+32),   S16((x)+48)
38#define S256(x)  S64((x)),  S64((x)+64),   S64((x)+128),  S64((x)+192)
39#define S1024(x) S256((x)), S256((x)+256), S256((x)+512), S256((x)+768)
40
41/* Lookup table-based algorithm from “Fast Half Float Conversions”
42 * by Jeroen van der Zijp, November 2008. No rounding is performed,
43 * and some NaN values may be incorrectly converted to Inf (because
44 * the lowest order bits in the mantissa are ignored). */
45static inline uint16_t float_to_half_nobranch(uint32_t x)
46{
47    static uint16_t const basetable[512] =
48    {
49#define S1(i) (((i) < 103) ? 0x0000u : \
50               ((i) < 113) ? 0x0400u >> (0x1f & (113 - (i))) : \
51               ((i) < 143) ? ((i) - 112) << 10 : 0x7c00u)
52        S256(0),
53#undef S1
54#define S1(i) (uint16_t)(0x8000u | basetable[i])
55        S256(0),
56#undef S1
57    };
58
59    static uint8_t const shifttable[512] =
60    {
61#define S1(i) (((i) < 103) ? 24 : \
62               ((i) < 113) ? 126 - (i) : \
63               ((i) < 143 || (i) == 255) ? 13 : 24)
64        S256(0), S256(0),
65#undef S1
66    };
67
68    uint16_t bits = basetable[(x >> 23) & 0x1ff];
69    bits |= (x & 0x007fffff) >> shifttable[(x >> 23) & 0x1ff];
70    return bits;
71}
72
73/* This method is faster than the OpenEXR implementation (very often
74 * used, eg. in Ogre), with the additional benefit of rounding, inspired
75 * by James Tursa’s half-precision code. */
76static inline uint16_t float_to_half_branch(uint32_t x)
77{
78    uint16_t bits = (x >> 16) & 0x8000; /* Get the sign */
79    uint16_t m = (x >> 12) & 0x07ff; /* Keep one extra bit for rounding */
80    unsigned int e = (x >> 23) & 0xff; /* Using int is faster here */
81
82    /* If zero, or denormal, or exponent underflows too much for a denormal
83     * half, return signed zero. */
84    if (e < 103)
85        return bits;
86
87    /* If NaN, return NaN. If Inf or exponent overflow, return Inf. */
88    if (e > 142)
89    {
90        bits |= 0x7c00u;
91        /* If exponent was 0xff and one mantissa bit was set, it means NaN,
92         * not Inf, so make sure we set one mantissa bit too. */
93        bits |= e == 255 && (x & 0x007fffffu);
94        return bits;
95    }
96
97    /* If exponent underflows but not too much, return a denormal */
98    if (e < 113)
99    {
100        m |= 0x0800u;
101        /* Extra rounding may overflow and set mantissa to 0 and exponent
102         * to 1, which is OK. */
103        bits |= (m >> (114 - e)) + ((m >> (113 - e)) & 1);
104        return bits;
105    }
106
107    bits |= ((e - 112) << 10) | (m >> 1);
108    /* Extra rounding. An overflow will set mantissa to 0 and increment
109     * the exponent, which is OK. */
110    bits += m & 1;
111    return bits;
112}
113
114/* We use this magic table, inspired by De Bruijn sequences, to compute a
115 * branchless integer log2. The actual value fetched is 24-log2(x+1) for x
116 * in 1, 3, 7, f, 1f, 3f, 7f, ff, 1fe, 1ff, 3fc, 3fd, 3fe, 3ff. See
117 * http://lol.zoy.org/blog/2012/4/3/beyond-de-bruijn for an explanation
118 * of how the value 0x5a1a1a2u was obtained. */
119static uint32_t const shifttable[16] =
120{
121    23, 22, 21, 15, 0, 20, 18, 14, 14, 16, 19, 0, 17, 0, 0, 0,
122};
123static uint32_t const shiftmagic = 0x5a1a1a2u;
124
125/* Lookup table-based algorithm from “Fast Half Float Conversions”
126 * by Jeroen van der Zijp, November 2008. Tables are generated using
127 * the C++ preprocessor, thanks to a branchless implementation also
128 * used in half_to_float_branch(). This code is very fast when performing
129 * conversions on arrays of values. */
130static inline uint32_t half_to_float_nobranch(uint16_t x)
131{
132#define M3(i) ((i) | ((i) >> 1))
133#define M7(i) (M3(i) | (M3(i) >> 2))
134#define MF(i) (M7(i) | (M7(i) >> 4))
135#define E(i) shifttable[(uint32_t)((uint64_t)MF(i) * shiftmagic) >> 28]
136
137    static uint32_t const mantissatable[2048] =
138    {
139#define S1(i) (((i) == 0) ? 0 : ((125 - E(i)) << 23) + ((i) << E(i)))
140        S1024(0),
141#undef S1
142#define S1(i) (0x38000000u + ((i) << 13))
143        S1024(0),
144#undef S1
145    };
146
147    static uint32_t const exponenttable[64] =
148    {
149#define S1(i) (((i) == 0) ? 0 : \
150               ((i) < 31) ? ((uint32_t)(i) << 23) : \
151               ((i) == 31) ? 0x47800000u : \
152               ((i) == 32) ? 0x80000000u : \
153               ((i) < 63) ? (0x80000000u | (((i) - 32) << 23)) : 0xc7800000)
154        S64(0),
155#undef S1
156    };
157
158    static int const offsettable[64] =
159    {
160#define S1(i) (((i) == 0 || (i) == 32) ? 0 : 1024)
161        S64(0),
162#undef S1
163    };
164
165    return mantissatable[offsettable[x >> 10] + (x & 0x3ff)]
166            + exponenttable[x >> 10];
167}
168
169/* This algorithm is similar to the OpenEXR implementation, except it
170 * uses branchless code in the denormal path. This is slower than the
171 * table version, but will be more friendly to the cache for occasional
172 * uses. */
173static inline uint32_t half_to_float_branch(uint16_t x)
174{
175    uint32_t s = (x & 0x8000u) << 16;
176
177    if ((x & 0x7fffu) == 0)
178        return (uint32_t)x << 16;
179
180    uint32_t e = x & 0x7c00u;
181    uint32_t m = x & 0x03ffu;
182
183    if (e == 0)
184    {
185        /* m has 10 significant bits but replicating the leading bit to
186         * 8 positions instead of 16 works just as well because of our
187         * handcrafted shiftmagic table. */
188        uint32_t v = m | (m >> 1);
189        v |= v >> 2;
190        v |= v >> 4;
191
192        e = shifttable[(v * shiftmagic) >> 28];
193
194        /* We don't have to remove the 10th mantissa bit because it gets
195         * added to our underestimated exponent. */
196        return s | (((125 - e) << 23) + (m << e));
197    }
198
199    if (e == 0x7c00u)
200    {
201        /* The amd64 pipeline likes the if() better than a ternary operator
202         * or any other trick I could find. --sam */
203        if (m == 0)
204            return s | 0x7f800000u;
205        return s | 0x7fc00000u;
206    }
207
208    return s | (((e >> 10) + 112) << 23) | (m << 13);
209}
210
211/* Constructor from float. Uses the non-branching version because benchmarks
212 * indicate it is about 80% faster on amd64, and 20% faster on the PS3. The
213 * penalty of loading the lookup tables does not seem important. */
214half half::makefast(float f)
215{
216    union { float f; uint32_t x; } u = { f };
217    return makebits(float_to_half_nobranch(u.x));
218}
219
220/* Constructor from float with better precision. */
221half half::makeaccurate(float f)
222{
223    union { float f; uint32_t x; } u = { f };
224    return makebits(float_to_half_branch(u.x));
225}
226
227/* Cast to float. Uses the branching version because loading the tables
228 * for only one value is going to be cache-expensive. */
229float half::tofloat(half h)
230{
231    union { float f; uint32_t x; } u;
232    u.x = half_to_float_branch(h.bits);
233    return u.f;
234}
235
236size_t half::convert(half *dst, float const *src, size_t nelem)
237{
238    for (size_t i = 0; i < nelem; i++)
239    {
240        union { float f; uint32_t x; } u;
241        u.f = *src++;
242        *dst++ = makebits(float_to_half_nobranch(u.x));
243    }
244
245    return nelem;
246}
247
248size_t half::convert(float *dst, half const *src, size_t nelem)
249{
250    for (size_t i = 0; i < nelem; i++)
251    {
252        union { float f; uint32_t x; } u;
253#if !defined __CELLOS_LV2__
254        /* This code is really too slow on the PS3, even with the denormal
255         * handling stripped off. */
256        u.x = half_to_float_nobranch((*src++).bits);
257#else
258        u.x = half_to_float_branch((*src++).bits);
259#endif
260        *dst++ = u.f;
261    }
262
263    return nelem;
264}
265
266} /* namespace lol */
267
Note: See TracBrowser for help on using the repository browser.