Merge branch 'master' into blender2.8
[blender.git] / source / blender / blenlib / intern / math_base_inline.c
1 /*
2  * ***** BEGIN GPL LICENSE BLOCK *****
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  *
18  * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
19  * All rights reserved.
20  *
21  * The Original Code is: some of this file.
22  *
23  * ***** END GPL LICENSE BLOCK *****
24  * */
25
26 /** \file blender/blenlib/intern/math_base_inline.c
27  *  \ingroup bli
28  */
29
30 #ifndef __MATH_BASE_INLINE_C__
31 #define __MATH_BASE_INLINE_C__
32
33 #include <float.h>
34 #include <stdio.h>
35 #include <stdlib.h>
36 #include <limits.h>
37
38 #ifdef __SSE2__
39 #  include <emmintrin.h>
40 #endif
41
42 #include "BLI_math_base.h"
43
44 /* copied from BLI_utildefines.h */
45 #ifdef __GNUC__
46 #  define UNLIKELY(x)     __builtin_expect(!!(x), 0)
47 #else
48 #  define UNLIKELY(x)     (x)
49 #endif
50
51 /* powf is really slow for raising to integer powers. */
52 MINLINE float pow2f(float x)
53 {
54         return x * x;
55 }
56 MINLINE float pow3f(float x)
57 {
58         return pow2f(x) * x;
59 }
60 MINLINE float pow4f(float x)
61 {
62         return pow2f(pow2f(x));
63 }
64 MINLINE float pow7f(float x)
65 {
66         return pow2f(pow3f(x)) * x;
67 }
68
69 MINLINE float sqrt3f(float f)
70 {
71         if      (UNLIKELY(f == 0.0f)) return 0.0f;
72         else if (UNLIKELY(f <  0.0f)) return -(float)(exp(log(-f) / 3.0));
73         else                          return  (float)(exp(log( f) / 3.0));
74 }
75
76 MINLINE double sqrt3d(double d)
77 {
78         if      (UNLIKELY(d == 0.0)) return 0.0;
79         else if (UNLIKELY(d <  0.0)) return -exp(log(-d) / 3.0);
80         else                         return  exp(log( d) / 3.0);
81 }
82
83 MINLINE float sqrtf_signed(float f)
84 {
85         return (f >= 0.0f) ? sqrtf(f) : -sqrtf(-f);
86 }
87
88 MINLINE float saacos(float fac)
89 {
90         if      (UNLIKELY(fac <= -1.0f)) return (float)M_PI;
91         else if (UNLIKELY(fac >=  1.0f)) return 0.0f;
92         else                             return acosf(fac);
93 }
94
95 MINLINE float saasin(float fac)
96 {
97         if      (UNLIKELY(fac <= -1.0f)) return (float)-M_PI / 2.0f;
98         else if (UNLIKELY(fac >=  1.0f)) return (float) M_PI / 2.0f;
99         else                             return asinf(fac);
100 }
101
102 MINLINE float sasqrt(float fac)
103 {
104         if (UNLIKELY(fac <= 0.0f)) return 0.0f;
105         else                       return sqrtf(fac);
106 }
107
108 MINLINE float saacosf(float fac)
109 {
110         if      (UNLIKELY(fac <= -1.0f)) return (float)M_PI;
111         else if (UNLIKELY(fac >=  1.0f)) return 0.0f;
112         else                             return acosf(fac);
113 }
114
115 MINLINE float saasinf(float fac)
116 {
117         if      (UNLIKELY(fac <= -1.0f)) return (float)-M_PI / 2.0f;
118         else if (UNLIKELY(fac >=  1.0f)) return (float) M_PI / 2.0f;
119         else                             return asinf(fac);
120 }
121
122 MINLINE float sasqrtf(float fac)
123 {
124         if (UNLIKELY(fac <= 0.0f)) return 0.0f;
125         else                       return sqrtf(fac);
126 }
127
128 MINLINE float interpf(float target, float origin, float fac)
129 {
130         return (fac * target) + (1.0f - fac) * origin;
131 }
132
133 /* used for zoom values*/
134 MINLINE float power_of_2(float val)
135 {
136         return (float)pow(2.0, ceil(log((double)val) / M_LN2));
137 }
138
139 MINLINE int is_power_of_2_i(int n)
140 {
141         return (n & (n - 1)) == 0;
142 }
143
144 MINLINE int power_of_2_max_i(int n)
145 {
146         if (is_power_of_2_i(n))
147                 return n;
148
149         do {
150                 n = n & (n - 1);
151         } while (!is_power_of_2_i(n));
152
153         return n * 2;
154 }
155
156 MINLINE int power_of_2_min_i(int n)
157 {
158         while (!is_power_of_2_i(n))
159                 n = n & (n - 1);
160
161         return n;
162 }
163
164 MINLINE unsigned int power_of_2_max_u(unsigned int x)
165 {
166         x -= 1;
167         x |= (x >>  1);
168         x |= (x >>  2);
169         x |= (x >>  4);
170         x |= (x >>  8);
171         x |= (x >> 16);
172         return x + 1;
173 }
174
175 MINLINE unsigned power_of_2_min_u(unsigned x)
176 {
177         x |= (x >>  1);
178         x |= (x >>  2);
179         x |= (x >>  4);
180         x |= (x >>  8);
181         x |= (x >> 16);
182         return x - (x >> 1);
183 }
184
185 /* rounding and clamping */
186
187 MINLINE int iroundf(float a)
188 {
189         return (int)floorf(a + 0.5f);
190 }
191
192 #define _round_clamp_fl_impl(arg, ty, min, max) { \
193         float r = floorf(arg + 0.5f); \
194         if      (UNLIKELY(r <= (float)min)) return (ty)min; \
195         else if (UNLIKELY(r >= (float)max)) return (ty)max; \
196         else return (ty)r; \
197 }
198
199 #define _round_clamp_db_impl(arg, ty, min, max) { \
200         double r = floor(arg + 0.5); \
201         if      (UNLIKELY(r <= (double)min)) return (ty)min; \
202         else if (UNLIKELY(r >= (double)max)) return (ty)max; \
203         else return (ty)r; \
204 }
205
206 MINLINE signed char    round_fl_to_char_clamp(float a) { _round_clamp_fl_impl(a, signed char, SCHAR_MIN, SCHAR_MAX) }
207 MINLINE unsigned char  round_fl_to_uchar_clamp(float a) { _round_clamp_fl_impl(a, unsigned char, 0, UCHAR_MAX) }
208 MINLINE short          round_fl_to_short_clamp(float a) { _round_clamp_fl_impl(a, short, SHRT_MIN, SHRT_MAX) }
209 MINLINE unsigned short round_fl_to_ushort_clamp(float a) { _round_clamp_fl_impl(a, unsigned short, 0, USHRT_MAX) }
210 MINLINE int            round_fl_to_int_clamp(float a) { _round_clamp_fl_impl(a, int, INT_MIN, INT_MAX) }
211 MINLINE unsigned int   round_fl_to_uint_clamp(float a) { _round_clamp_fl_impl(a, unsigned int, 0, UINT_MAX) }
212
213 MINLINE signed char    round_db_to_char_clamp(double a) { _round_clamp_db_impl(a, signed char, SCHAR_MIN, SCHAR_MAX) }
214 MINLINE unsigned char  round_db_to_uchar_clamp(double a) { _round_clamp_db_impl(a, unsigned char, 0, UCHAR_MAX) }
215 MINLINE short          round_db_to_short_clamp(double a) { _round_clamp_db_impl(a, short, SHRT_MIN, SHRT_MAX) }
216 MINLINE unsigned short round_db_to_ushort_clamp(double a) { _round_clamp_db_impl(a, unsigned short, 0, USHRT_MAX) }
217 MINLINE int            round_db_to_int_clamp(double a) { _round_clamp_db_impl(a, int, INT_MIN, INT_MAX) }
218 MINLINE unsigned int   round_db_to_uint_clamp(double a) { _round_clamp_db_impl(a, unsigned int, 0, UINT_MAX) }
219
220 #undef _round_clamp_fl_impl
221 #undef _round_clamp_db_impl
222
223 /* integer division that rounds 0.5 up, particularly useful for color blending
224  * with integers, to avoid gradual darkening when rounding down */
225 MINLINE int divide_round_i(int a, int b)
226 {
227         return (2 * a + b) / (2 * b);
228 }
229
230 /**
231  * Integer division that floors negative result.
232  * \note This works like Python's int division.
233  */
234 MINLINE int divide_floor_i(int a, int b)
235 {
236         int d = a / b;
237         int r = a % b;  /* Optimizes into a single division. */
238         return r ? d - ((a < 0) ^ (b < 0)) : d;
239 }
240
241 /**
242  * modulo that handles negative numbers, works the same as Python's.
243  */
244 MINLINE int mod_i(int i, int n)
245 {
246         return (i % n + n) % n;
247 }
248
249 MINLINE float min_ff(float a, float b)
250 {
251         return (a < b) ? a : b;
252 }
253 MINLINE float max_ff(float a, float b)
254 {
255         return (a > b) ? a : b;
256 }
257
258 MINLINE int min_ii(int a, int b)
259 {
260         return (a < b) ? a : b;
261 }
262 MINLINE int max_ii(int a, int b)
263 {
264         return (b < a) ? a : b;
265 }
266
267 MINLINE float min_fff(float a, float b, float c)
268 {
269         return min_ff(min_ff(a, b), c);
270 }
271 MINLINE float max_fff(float a, float b, float c)
272 {
273         return max_ff(max_ff(a, b), c);
274 }
275
276 MINLINE int min_iii(int a, int b, int c)
277 {
278         return min_ii(min_ii(a, b), c);
279 }
280 MINLINE int max_iii(int a, int b, int c)
281 {
282         return max_ii(max_ii(a, b), c);
283 }
284
285 MINLINE float min_ffff(float a, float b, float c, float d)
286 {
287         return min_ff(min_fff(a, b, c), d);
288 }
289 MINLINE float max_ffff(float a, float b, float c, float d)
290 {
291         return max_ff(max_fff(a, b, c), d);
292 }
293
294 MINLINE int min_iiii(int a, int b, int c, int d)
295 {
296         return min_ii(min_iii(a, b, c), d);
297 }
298 MINLINE int max_iiii(int a, int b, int c, int d)
299 {
300         return max_ii(max_iii(a, b, c), d);
301 }
302
303 /**
304  * Almost-equal for IEEE floats, using absolute difference method.
305  *
306  * \param max_diff the maximum absolute difference.
307  */
308 MINLINE int compare_ff(float a, float b, const float max_diff)
309 {
310         return fabsf(a - b) <= max_diff;
311 }
312
313 /**
314  * Almost-equal for IEEE floats, using their integer representation (mixing ULP and absolute difference methods).
315  *
316  * \param max_diff is the maximum absolute difference (allows to take care of the near-zero area,
317  *                 where relative difference methods cannot really work).
318  * \param max_ulps is the 'maximum number of floats + 1' allowed between \a a and \a b to consider them equal.
319  *
320  * \see https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
321  */
322 MINLINE int compare_ff_relative(float a, float b, const float max_diff, const int max_ulps)
323 {
324         union {float f; int i;} ua, ub;
325
326 #if 0  /* No BLI_assert in INLINE :/ */
327         BLI_assert(sizeof(float) == sizeof(int));
328         BLI_assert(max_ulps < (1 << 22));
329 #endif
330
331         if (fabsf(a - b) <= max_diff) {
332                 return 1;
333         }
334
335         ua.f = a;
336         ub.f = b;
337
338         /* Important to compare sign from integers, since (-0.0f < 0) is false
339          * (though this shall not be an issue in common cases)... */
340         return ((ua.i < 0) != (ub.i < 0)) ? 0 : (abs(ua.i - ub.i) <= max_ulps) ? 1 : 0;
341 }
342
343 MINLINE float signf(float f)
344 {
345         return (f < 0.f) ? -1.f : 1.f;
346 }
347
348 MINLINE int signum_i_ex(float a, float eps)
349 {
350         if (a >  eps) return  1;
351         if (a < -eps) return -1;
352         else          return  0;
353 }
354
355 MINLINE int signum_i(float a)
356 {
357         if (a > 0.0f) return  1;
358         if (a < 0.0f) return -1;
359         else          return  0;
360 }
361
362 /** Returns number of (base ten) *significant* digits of integer part of given float
363  * (negative in case of decimal-only floats, 0.01 returns -1 e.g.). */
364 MINLINE int integer_digits_f(const float f)
365 {
366         return (f == 0.0f) ? 0 : (int)floor(log10(fabs(f))) + 1;
367 }
368
369 /** Returns number of (base ten) *significant* digits of integer part of given double
370  * (negative in case of decimal-only floats, 0.01 returns -1 e.g.). */
371 MINLINE int integer_digits_d(const double d)
372 {
373         return (d == 0.0) ? 0 : (int)floor(log10(fabs(d))) + 1;
374 }
375
376
377 /* Internal helpers for SSE2 implementation.
378  *
379  * NOTE: Are to be called ONLY from inside `#ifdef __SSE2__` !!!
380  */
381
382 #ifdef __SSE2__
383
384 /* Calculate initial guess for arg^exp based on float representation
385  * This method gives a constant bias, which can be easily compensated by
386  * multiplicating with bias_coeff.
387  * Gives better results for exponents near 1 (e. g. 4/5).
388  * exp = exponent, encoded as uint32_t
389  * e2coeff = 2^(127/exponent - 127) * bias_coeff^(1/exponent), encoded as
390  * uint32_t
391  *
392  * We hope that exp and e2coeff gets properly inlined
393  */
394 MALWAYS_INLINE __m128 _bli_math_fastpow(const int exp,
395                                         const int e2coeff,
396                                         const __m128 arg)
397 {
398         __m128 ret;
399         ret = _mm_mul_ps(arg, _mm_castsi128_ps(_mm_set1_epi32(e2coeff)));
400         ret = _mm_cvtepi32_ps(_mm_castps_si128(ret));
401         ret = _mm_mul_ps(ret, _mm_castsi128_ps(_mm_set1_epi32(exp)));
402         ret = _mm_castsi128_ps(_mm_cvtps_epi32(ret));
403         return ret;
404 }
405
406 /* Improve x ^ 1.0f/5.0f solution with Newton-Raphson method */
407 MALWAYS_INLINE __m128 _bli_math_improve_5throot_solution(
408     const __m128 old_result,
409     const __m128 x)
410 {
411         __m128 approx2 = _mm_mul_ps(old_result, old_result);
412         __m128 approx4 = _mm_mul_ps(approx2, approx2);
413         __m128 t = _mm_div_ps(x, approx4);
414         __m128 summ = _mm_add_ps(_mm_mul_ps(_mm_set1_ps(4.0f), old_result), t); /* fma */
415         return _mm_mul_ps(summ, _mm_set1_ps(1.0f / 5.0f));
416 }
417
418 /* Calculate powf(x, 2.4). Working domain: 1e-10 < x < 1e+10 */
419 MALWAYS_INLINE __m128 _bli_math_fastpow24(const __m128 arg)
420 {
421         /* max, avg and |avg| errors were calculated in gcc without FMA instructions
422          * The final precision should be better than powf in glibc */
423
424         /* Calculate x^4/5, coefficient 0.994 was constructed manually to minimize
425          * avg error.
426          */
427         /* 0x3F4CCCCD = 4/5 */
428         /* 0x4F55A7FB = 2^(127/(4/5) - 127) * 0.994^(1/(4/5)) */
429         /* error max = 0.17     avg = 0.0018    |avg| = 0.05 */
430         __m128 x = _bli_math_fastpow(0x3F4CCCCD, 0x4F55A7FB, arg);
431         __m128 arg2 = _mm_mul_ps(arg, arg);
432         __m128 arg4 = _mm_mul_ps(arg2, arg2);
433         /* error max = 0.018        avg = 0.0031    |avg| = 0.0031  */
434         x = _bli_math_improve_5throot_solution(x, arg4);
435         /* error max = 0.00021    avg = 1.6e-05    |avg| = 1.6e-05 */
436         x = _bli_math_improve_5throot_solution(x, arg4);
437         /* error max = 6.1e-07    avg = 5.2e-08    |avg| = 1.1e-07 */
438         x = _bli_math_improve_5throot_solution(x, arg4);
439         return _mm_mul_ps(x, _mm_mul_ps(x, x));
440 }
441
442 /* Calculate powf(x, 1.0f / 2.4) */
443 MALWAYS_INLINE __m128 _bli_math_fastpow512(const __m128 arg)
444 {
445         /* 5/12 is too small, so compute the 4th root of 20/12 instead.
446          * 20/12 = 5/3 = 1 + 2/3 = 2 - 1/3. 2/3 is a suitable argument for fastpow.
447          * weighting coefficient: a^-1/2 = 2 a; a = 2^-2/3
448          */
449         __m128 xf = _bli_math_fastpow(0x3f2aaaab, 0x5eb504f3, arg);
450         __m128 xover = _mm_mul_ps(arg, xf);
451         __m128 xfm1 = _mm_rsqrt_ps(xf);
452         __m128 x2 = _mm_mul_ps(arg, arg);
453         __m128 xunder = _mm_mul_ps(x2, xfm1);
454         /* sqrt2 * over + 2 * sqrt2 * under */
455         __m128 xavg = _mm_mul_ps(_mm_set1_ps(1.0f / (3.0f * 0.629960524947437f) * 0.999852f),
456                                  _mm_add_ps(xover, xunder));
457         xavg = _mm_mul_ps(xavg, _mm_rsqrt_ps(xavg));
458         xavg = _mm_mul_ps(xavg, _mm_rsqrt_ps(xavg));
459         return xavg;
460 }
461
462 MALWAYS_INLINE __m128 _bli_math_blend_sse(const __m128 mask,
463                                           const __m128 a,
464                                           const __m128 b)
465 {
466         return _mm_or_ps(_mm_and_ps(mask, a), _mm_andnot_ps(mask, b));
467 }
468
469 #endif  /* __SSE2__ */
470
471 #endif /* __MATH_BASE_INLINE_C__ */