Line data Source code
1 :
2 : /**
3 : @brief fast math library for float
4 : @author herumi
5 : @url http://homepage1.nifty.com/herumi/
6 : @note modified new BSD license
7 : http://opensource.org/licenses/BSD-3-Clause
8 :
9 : cl /Ox /Ob2 /arch:SSE2 /fp:fast bench.cpp -I../xbyak /EHsc /DNOMINMAX
10 : g++ -O3 -fomit-frame-pointer -fno-operator-names -march=core2 -mssse3 -mfpmath=sse -ffast-math -fexcess-precision=fast
11 : */
12 : /*
13 : function prototype list
14 :
15 : float fmath::exp(float);
16 : double fmath::expd(double);
17 : float fmath::log(float);
18 :
19 : __m128 fmath::exp_ps(__m128);
20 : __m256 fmath::exp_ps256(__m256);
21 : __m128 fmath::log_ps(__m128);
22 :
23 : double fmath::expd_v(double *, size_t n);
24 :
25 : if FMATH_USE_XBYAK is defined then Xbyak version are used
26 : */
27 :
28 : #ifndef ATK_UTILITY_FMATH_H
29 : #define ATK_UTILITY_FMATH_H
30 :
31 : #include <algorithm>
32 : #include <cmath>
33 : #include <cstddef>
34 : #include <cassert>
35 : #include <limits>
36 : #include <cstdint>
37 : #include <cfloat>
38 : #include <cstring> // for memcpy
39 :
40 : #if defined(_WIN32) && !defined(__GNUC__)
41 : # include <intrin.h>
42 : # ifndef MIE_ALIGN
43 : # define MIE_ALIGN(x) __declspec(align(x))
44 : # endif
45 : #elif defined(__GNUC__) && (defined(__x86_64__) || defined(__i386__))
46 : /* GCC >= 4.4 and non-GCC compilers */
47 : # include <x86intrin.h>
48 : #endif
49 : #ifndef MIE_ALIGN
50 : # define MIE_ALIGN(x) __attribute__((aligned(x)))
51 : #endif
52 :
53 : namespace fmath {
54 :
55 : namespace local {
56 :
57 : constexpr size_t EXP_TABLE_SIZE = 10;
58 : constexpr size_t EXPD_TABLE_SIZE = 11;
59 : constexpr size_t LOG_TABLE_SIZE = 12;
60 :
61 : union fi {
62 : float f;
63 : unsigned int i;
64 : };
65 :
66 : union di {
67 : double d;
68 : uint64_t i;
69 : };
70 :
71 18043717 : inline unsigned int mask(int x)
72 : {
73 18043717 : return (1U << x) - 1;
74 : }
75 :
76 497664 : inline uint64_t mask64(int x)
77 : {
78 497664 : return (1ULL << x) - 1;
79 : }
80 :
81 : template<class T>
82 : inline const T* cast_to(const void *p)
83 : {
84 : return reinterpret_cast<const T*>(p);
85 : }
86 :
87 : template<class T, size_t N>
88 : size_t NumOfArray(const T (&)[N]) { return N; }
89 :
90 : /*
91 : exp(88.722839F) = inf ; 0x42b17218
92 : exp(-87.33655F) = 1.175491e-038f(007fffe6) denormal ; 0xc2aeac50
93 : exp(-103.972081F) = 0 ; 0xc2cff1b5
94 : */
95 : template<size_t N = EXP_TABLE_SIZE>
96 : class ExpVar {
97 : public:
98 : enum {
99 : s = N,
100 : n = 1 << s,
101 : f88 = 0x42b00000 /* 88.0 */
102 : };
103 : float minX[8];
104 : float maxX[8];
105 : float a[8];
106 : float b[8];
107 : float f1[8];
108 : unsigned int i127s[8];
109 : unsigned int mask_s[8];
110 : unsigned int i7fffffff[8];
111 : unsigned int tbl[n];
112 : ExpVar()
113 : {
114 : float log_2 = std::log(2.0F);
115 : for (int i = 0; i < 8; i++) {
116 : maxX[i] = 88;
117 : minX[i] = -88;
118 : a[i] = n / log_2;
119 : b[i] = log_2 / n;
120 : f1[i] = 1.0F;
121 : i127s[i] = 127 << s;
122 : i7fffffff[i] = 0x7fffffff;
123 : mask_s[i] = mask(s);
124 : }
125 :
126 : for (int i = 0; i < n; i++) {
127 : float y = std::pow(2.0F, static_cast<float>(i) / n);
128 : fi fi;
129 : fi.f = y;
130 : tbl[i] = fi.i & mask(23);
131 : }
132 : }
133 : };
134 :
135 : template<size_t sbit_ = EXPD_TABLE_SIZE>
136 : class ExpdVar {
137 : public:
138 : enum {
139 : sbit = sbit_,
140 : s = 1UL << sbit,
141 : adj = (1UL << (sbit + 10)) - (1UL << sbit)
142 : };
143 : // A = 1, B = 1, C = 1/2, D = 1/6
144 : double C1[2]{}; // A
145 : double C2[2]{}; // D
146 : double C3[2]{}; // C/D
147 : uint64_t tbl[s]{};
148 : double a{};
149 : double ra{};
150 243 : ExpdVar()
151 : : a(s / ::log(2.0))
152 243 : , ra(1 / a)
153 : {
154 729 : for (int i = 0; i < 2; i++) {
155 486 : C1[i] = 1.0;
156 486 : C2[i] = 0.16666666685227835064;
157 486 : C3[i] = 3.0000000027955394;
158 : }
159 497907 : for (int i = 0; i < s; i++) {
160 : di di;
161 497664 : di.d = ::pow(2.0, i * (1.0 / s));
162 497664 : tbl[i] = di.i & mask64(52);
163 : }
164 243 : }
165 : };
166 :
167 : template<size_t N = LOG_TABLE_SIZE>
168 : class LogVar {
169 : public:
170 : enum {
171 : LEN = N - 1
172 : };
173 : unsigned int m1[4]; // 0
174 : unsigned int m2[4]; // 16
175 : unsigned int m3[4]; // 32
176 : float m4[4]; // 48
177 : unsigned int m5[4]; // 64
178 : struct {
179 : float app;
180 : float rev;
181 : } tbl[1 << LEN];
182 : float c_log2;
183 243 : LogVar()
184 243 : : c_log2(std::log(2.0F) / (1 << 23))
185 : {
186 243 : const double e = 1 / static_cast<double>(1 << 24);
187 243 : const double h = 1 / static_cast<double>(1 << LEN);
188 243 : constexpr size_t n = 1U << LEN;
189 497907 : for (size_t i = 0; i < n; i++) {
190 497664 : double x = 1 + static_cast<double>(i) / n;
191 497664 : double a = std::log(x);
192 497664 : tbl[i].app = static_cast<float>(a);
193 497664 : if (i < n - 1) {
194 497421 : double b = std::log(x + h - e);
195 497421 : tbl[i].rev = static_cast<float>((b - a) / ((h - e) * (1 << 23)));
196 : } else {
197 243 : tbl[i].rev = static_cast<float>(1 / (x * (1 << 23)));
198 : }
199 : }
200 1215 : for (int i = 0; i < 4; i++) {
201 972 : m1[i] = mask(8) << 23;
202 972 : m2[i] = mask(LEN) << (23 - LEN);
203 972 : m3[i] = mask(23 - LEN);
204 972 : m4[i] = c_log2;
205 972 : m5[i] = 127U << 23;
206 : }
207 243 : }
208 : };
209 : /* to define static variables in fmath.hpp */
210 : template<size_t EXP_N = EXP_TABLE_SIZE, size_t LOG_N = LOG_TABLE_SIZE, size_t EXPD_N = EXPD_TABLE_SIZE>
211 : struct C {
212 : static const ExpVar<EXP_N> expVar;
213 : static const LogVar<LOG_N> logVar;
214 : static const ExpdVar<EXPD_N> expdVar;
215 : };
216 :
217 : template<size_t EXP_N, size_t LOG_N, size_t EXPD_N>
218 : MIE_ALIGN(32) const ExpVar<EXP_N> C<EXP_N, LOG_N, EXPD_N>::expVar;
219 :
220 : template<size_t EXP_N, size_t LOG_N, size_t EXPD_N>
221 : MIE_ALIGN(32) const LogVar<LOG_N> C<EXP_N, LOG_N, EXPD_N>::logVar;
222 :
223 : template<size_t EXP_N, size_t LOG_N, size_t EXPD_N>
224 : MIE_ALIGN(32) const ExpdVar<EXPD_N> C<EXP_N, LOG_N, EXPD_N>::expdVar;
225 :
226 : } // fmath::local
227 :
228 : #if defined(ATK_NO_FMATH) || defined(__GNUC__) && !(defined(__x86_64__) || defined(__i386__))
229 : using std::exp;
230 : using std::log;
231 : #else
232 21459496 : inline double exp(double x)
233 : {
234 21459496 : if (x <= -708.39641853226408)
235 : {
236 3025571 : return 0;
237 : }
238 18433945 : if (x >= 709.78271289338397)
239 : {
240 393215 : return std::numeric_limits<double>::infinity();
241 : }
242 :
243 18040745 : const auto& c = local::C<>::expdVar;
244 18040745 : const double _b = static_cast<double>(uint64_t(3) << 51);
245 18040745 : __m128d b = _mm_load_sd(&_b);
246 18040745 : __m128d xx = _mm_load_sd(&x);
247 72162882 : __m128d d = _mm_add_sd(_mm_mul_sd(xx, _mm_load_sd(&c.a)), b);
248 18040745 : uint64_t di = _mm_cvtsi128_si32(_mm_castpd_si128(d));
249 18040745 : uint64_t iax = c.tbl[di & mask(c.sbit)];
250 72163282 : __m128d _t = _mm_sub_sd(_mm_mul_sd(_mm_sub_sd(d, b), _mm_load_sd(&c.ra)), xx);
251 18040845 : uint64_t u = ((di + c.adj) >> c.sbit) << 52;
252 : double t;
253 : _mm_store_sd(&t, _t);
254 18040845 : double y = (c.C3[0] - t) * (t * t) * c.C2[0] - t + c.C1[0];
255 : double did;
256 18040845 : u |= iax;
257 18040845 : memcpy(&did, &u, sizeof(did));
258 18040845 : return y * did;
259 : }
260 :
261 : inline float exp(float x)
262 : {
263 : return static_cast<float>(exp(static_cast<double>(x)));
264 : }
265 :
266 : inline float log(float x)
267 : {
268 : const auto& logVar = local::C<>::logVar;
269 : constexpr size_t logLen = logVar.LEN;
270 :
271 : local::fi fi;
272 : fi.f = x;
273 : int a = fi.i & (local::mask(8) << 23);
274 : unsigned int b1 = fi.i & (local::mask(logLen) << (23 - logLen));
275 : unsigned int b2 = fi.i & local::mask(23 - logLen);
276 : int idx = b1 >> (23 - logLen);
277 : float f = static_cast<float>(a - (127 << 23)) * logVar.c_log2 + logVar.tbl[idx].app + static_cast<float>(b2) * logVar.tbl[idx].rev;
278 : return f;
279 : }
280 :
281 : inline float log2(float x) { return fmath::log(x) * 1.442695F; }
282 :
283 40206489 : inline double log(double x)
284 : {
285 40206489 : return std::log(x);
286 : }
287 : #endif
288 :
289 16432752 : inline double pow(double x, double y)
290 : {
291 16432752 : return exp(y * fmath::log(x));
292 : }
293 :
294 : inline float pow(float x, float y)
295 : {
296 : return exp(y * fmath::log(x));
297 : }
298 :
299 : template<typename T>
300 11799100 : T log10(T x)
301 : {
302 11799100 : return log(x) / log(static_cast<T>(10));
303 : }
304 :
305 : /*
306 : for given y > 0
307 : get f_y(x) := pow(x, y) for x >= 0
308 : */
309 : class PowGenerator {
310 : enum {
311 : N = 11
312 : };
313 : float tbl0_[256]{};
314 : struct {
315 : float app;
316 : float rev;
317 : } tbl1_[1 << N]{};
318 : public:
319 : explicit PowGenerator(float y)
320 : {
321 : for (int i = 0; i < 256; i++) {
322 : tbl0_[i] = static_cast<float>(std::pow(2, (i - 127) * y));
323 : }
324 : const double e = 1 / static_cast<double>(1 << 24);
325 : const double h = 1 / static_cast<double>(1 << N);
326 : constexpr size_t n = 1U << N;
327 : for (size_t i = 0; i < n; i++) {
328 : double x = 1 + static_cast<double>(i) / n;
329 : double a = std::pow(x, static_cast<double>(y));
330 : tbl1_[i].app = static_cast<float>(a);
331 : double b = std::pow(x + h - e, static_cast<double>(y));
332 : tbl1_[i].rev = static_cast<float>((b - a) / (h - e) / (1 << 23));
333 : }
334 : }
335 : float get(float x) const
336 : {
337 : local::fi fi;
338 : fi.f = x;
339 : int a = (fi.i >> 23) & local::mask(8);
340 : unsigned int b = fi.i & local::mask(23);
341 : unsigned int b1 = b & (local::mask(N) << (23 - N));
342 : unsigned int b2 = b & local::mask(23 - N);
343 : float f;
344 : int idx = b1 >> (23 - N);
345 : f = tbl0_[a] * (tbl1_[idx].app + static_cast<float>(b2) * tbl1_[idx].rev);
346 : return f;
347 : }
348 : };
349 :
350 : // exp2(x) = pow(2, x)
351 : inline float exp2(float x) { return fmath::exp(x * 0.6931472F); }
352 :
353 : } // fmath
354 :
355 : #endif
|
