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