Cleanup: comments (long lines) in blenlib
[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
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)) {
69     return 0.0f;
70   }
71   else if (UNLIKELY(f < 0.0f)) {
72     return -(float)(exp(log(-f) / 3.0));
73   }
74   else {
75     return (float)(exp(log(f) / 3.0));
76   }
77 }
78
79 MINLINE double sqrt3d(double d)
80 {
81   if (UNLIKELY(d == 0.0)) {
82     return 0.0;
83   }
84   else if (UNLIKELY(d < 0.0)) {
85     return -exp(log(-d) / 3.0);
86   }
87   else {
88     return exp(log(d) / 3.0);
89   }
90 }
91
92 MINLINE float sqrtf_signed(float f)
93 {
94   return (f >= 0.0f) ? sqrtf(f) : -sqrtf(-f);
95 }
96
97 MINLINE float saacos(float fac)
98 {
99   if (UNLIKELY(fac <= -1.0f)) {
100     return (float)M_PI;
101   }
102   else if (UNLIKELY(fac >= 1.0f)) {
103     return 0.0f;
104   }
105   else {
106     return acosf(fac);
107   }
108 }
109
110 MINLINE float saasin(float fac)
111 {
112   if (UNLIKELY(fac <= -1.0f)) {
113     return (float)-M_PI / 2.0f;
114   }
115   else if (UNLIKELY(fac >= 1.0f)) {
116     return (float)M_PI / 2.0f;
117   }
118   else {
119     return asinf(fac);
120   }
121 }
122
123 MINLINE float sasqrt(float fac)
124 {
125   if (UNLIKELY(fac <= 0.0f)) {
126     return 0.0f;
127   }
128   else {
129     return sqrtf(fac);
130   }
131 }
132
133 MINLINE float saacosf(float fac)
134 {
135   if (UNLIKELY(fac <= -1.0f)) {
136     return (float)M_PI;
137   }
138   else if (UNLIKELY(fac >= 1.0f)) {
139     return 0.0f;
140   }
141   else {
142     return acosf(fac);
143   }
144 }
145
146 MINLINE float saasinf(float fac)
147 {
148   if (UNLIKELY(fac <= -1.0f)) {
149     return (float)-M_PI / 2.0f;
150   }
151   else if (UNLIKELY(fac >= 1.0f)) {
152     return (float)M_PI / 2.0f;
153   }
154   else {
155     return asinf(fac);
156   }
157 }
158
159 MINLINE float sasqrtf(float fac)
160 {
161   if (UNLIKELY(fac <= 0.0f)) {
162     return 0.0f;
163   }
164   else {
165     return sqrtf(fac);
166   }
167 }
168
169 MINLINE float interpf(float target, float origin, float fac)
170 {
171   return (fac * target) + (1.0f - fac) * origin;
172 }
173
174 /* used for zoom values*/
175 MINLINE float power_of_2(float val)
176 {
177   return (float)pow(2.0, ceil(log((double)val) / M_LN2));
178 }
179
180 MINLINE int is_power_of_2_i(int n)
181 {
182   return (n & (n - 1)) == 0;
183 }
184
185 MINLINE int power_of_2_max_i(int n)
186 {
187   if (is_power_of_2_i(n)) {
188     return n;
189   }
190
191   do {
192     n = n & (n - 1);
193   } while (!is_power_of_2_i(n));
194
195   return n * 2;
196 }
197
198 MINLINE int power_of_2_min_i(int n)
199 {
200   while (!is_power_of_2_i(n)) {
201     n = n & (n - 1);
202   }
203
204   return n;
205 }
206
207 MINLINE unsigned int power_of_2_max_u(unsigned int x)
208 {
209   x -= 1;
210   x |= (x >> 1);
211   x |= (x >> 2);
212   x |= (x >> 4);
213   x |= (x >> 8);
214   x |= (x >> 16);
215   return x + 1;
216 }
217
218 MINLINE unsigned power_of_2_min_u(unsigned x)
219 {
220   x |= (x >> 1);
221   x |= (x >> 2);
222   x |= (x >> 4);
223   x |= (x >> 8);
224   x |= (x >> 16);
225   return x - (x >> 1);
226 }
227
228 /* rounding and clamping */
229
230 #define _round_clamp_fl_impl(arg, ty, min, max) \
231   { \
232     float r = floorf(arg + 0.5f); \
233     if (UNLIKELY(r <= (float)min)) \
234       return (ty)min; \
235     else if (UNLIKELY(r >= (float)max)) \
236       return (ty)max; \
237     else \
238       return (ty)r; \
239   }
240
241 #define _round_clamp_db_impl(arg, ty, min, max) \
242   { \
243     double r = floor(arg + 0.5); \
244     if (UNLIKELY(r <= (double)min)) \
245       return (ty)min; \
246     else if (UNLIKELY(r >= (double)max)) \
247       return (ty)max; \
248     else \
249       return (ty)r; \
250   }
251
252 #define _round_fl_impl(arg, ty) \
253   { \
254     return (ty)floorf(arg + 0.5f); \
255   }
256 #define _round_db_impl(arg, ty) \
257   { \
258     return (ty)floor(arg + 0.5); \
259   }
260
261 MINLINE signed char round_fl_to_char(float a){_round_fl_impl(a, signed char)} MINLINE
262     unsigned char round_fl_to_uchar(float a){_round_fl_impl(a, unsigned char)} MINLINE
263     short round_fl_to_short(float a){_round_fl_impl(a, short)} MINLINE
264     unsigned short round_fl_to_ushort(float a){_round_fl_impl(a, unsigned short)} MINLINE
265     int round_fl_to_int(float a){_round_fl_impl(a, int)} MINLINE
266     unsigned int round_fl_to_uint(float a){_round_fl_impl(a, unsigned int)}
267
268 MINLINE signed char round_db_to_char(double a){_round_db_impl(a, signed char)} MINLINE
269     unsigned char round_db_to_uchar(double a){_round_db_impl(a, unsigned char)} MINLINE
270     short round_db_to_short(double a){_round_db_impl(a, short)} MINLINE
271     unsigned short round_db_to_ushort(double a){_round_db_impl(a, unsigned short)} MINLINE
272     int round_db_to_int(double a){_round_db_impl(a, int)} MINLINE
273     unsigned int round_db_to_uint(double a)
274 {
275   _round_db_impl(a, unsigned int)
276 }
277
278 #undef _round_fl_impl
279 #undef _round_db_impl
280
281 MINLINE signed char round_fl_to_char_clamp(float a){
282     _round_clamp_fl_impl(a, signed char, SCHAR_MIN, SCHAR_MAX)} MINLINE
283     unsigned char round_fl_to_uchar_clamp(float a){
284         _round_clamp_fl_impl(a, unsigned char, 0, UCHAR_MAX)} MINLINE
285     short round_fl_to_short_clamp(float a){
286         _round_clamp_fl_impl(a, short, SHRT_MIN, SHRT_MAX)} MINLINE
287     unsigned short round_fl_to_ushort_clamp(float a){
288         _round_clamp_fl_impl(a, unsigned short, 0, USHRT_MAX)} MINLINE
289     int round_fl_to_int_clamp(float a){_round_clamp_fl_impl(a, int, INT_MIN, INT_MAX)} MINLINE
290     unsigned int round_fl_to_uint_clamp(float a){
291         _round_clamp_fl_impl(a, unsigned int, 0, UINT_MAX)}
292
293 MINLINE signed char round_db_to_char_clamp(double a){
294     _round_clamp_db_impl(a, signed char, SCHAR_MIN, SCHAR_MAX)} MINLINE
295     unsigned char round_db_to_uchar_clamp(double a){
296         _round_clamp_db_impl(a, unsigned char, 0, UCHAR_MAX)} MINLINE
297     short round_db_to_short_clamp(double a){
298         _round_clamp_db_impl(a, short, SHRT_MIN, SHRT_MAX)} MINLINE
299     unsigned short round_db_to_ushort_clamp(double a){
300         _round_clamp_db_impl(a, unsigned short, 0, USHRT_MAX)} MINLINE
301     int round_db_to_int_clamp(double a){_round_clamp_db_impl(a, int, INT_MIN, INT_MAX)} MINLINE
302     unsigned int round_db_to_uint_clamp(double a)
303 {
304   _round_clamp_db_impl(a, unsigned int, 0, UINT_MAX)
305 }
306
307 #undef _round_clamp_fl_impl
308 #undef _round_clamp_db_impl
309
310 /* integer division that rounds 0.5 up, particularly useful for color blending
311  * with integers, to avoid gradual darkening when rounding down */
312 MINLINE int divide_round_i(int a, int b)
313 {
314   return (2 * a + b) / (2 * b);
315 }
316
317 /**
318  * Integer division that floors negative result.
319  * \note This works like Python's int division.
320  */
321 MINLINE int divide_floor_i(int a, int b)
322 {
323   int d = a / b;
324   int r = a % b; /* Optimizes into a single division. */
325   return r ? d - ((a < 0) ^ (b < 0)) : d;
326 }
327
328 /**
329  * modulo that handles negative numbers, works the same as Python's.
330  */
331 MINLINE int mod_i(int i, int n)
332 {
333   return (i % n + n) % n;
334 }
335
336 MINLINE float min_ff(float a, float b)
337 {
338   return (a < b) ? a : b;
339 }
340 MINLINE float max_ff(float a, float b)
341 {
342   return (a > b) ? a : b;
343 }
344
345 MINLINE int min_ii(int a, int b)
346 {
347   return (a < b) ? a : b;
348 }
349 MINLINE int max_ii(int a, int b)
350 {
351   return (b < a) ? a : b;
352 }
353
354 MINLINE float min_fff(float a, float b, float c)
355 {
356   return min_ff(min_ff(a, b), c);
357 }
358 MINLINE float max_fff(float a, float b, float c)
359 {
360   return max_ff(max_ff(a, b), c);
361 }
362
363 MINLINE int min_iii(int a, int b, int c)
364 {
365   return min_ii(min_ii(a, b), c);
366 }
367 MINLINE int max_iii(int a, int b, int c)
368 {
369   return max_ii(max_ii(a, b), c);
370 }
371
372 MINLINE float min_ffff(float a, float b, float c, float d)
373 {
374   return min_ff(min_fff(a, b, c), d);
375 }
376 MINLINE float max_ffff(float a, float b, float c, float d)
377 {
378   return max_ff(max_fff(a, b, c), d);
379 }
380
381 MINLINE int min_iiii(int a, int b, int c, int d)
382 {
383   return min_ii(min_iii(a, b, c), d);
384 }
385 MINLINE int max_iiii(int a, int b, int c, int d)
386 {
387   return max_ii(max_iii(a, b, c), d);
388 }
389
390 MINLINE size_t min_zz(size_t a, size_t b)
391 {
392   return (a < b) ? a : b;
393 }
394 MINLINE size_t max_zz(size_t a, size_t b)
395 {
396   return (b < a) ? a : b;
397 }
398
399 MINLINE int clamp_i(int value, int min, int max)
400 {
401   return min_ii(max_ii(value, min), max);
402 }
403
404 MINLINE float clamp_f(float value, float min, float max)
405 {
406   if (value > max) {
407     return max;
408   }
409   else if (value < min) {
410     return min;
411   }
412   return value;
413 }
414
415 MINLINE size_t clamp_z(size_t value, size_t min, size_t max)
416 {
417   return min_zz(max_zz(value, min), max);
418 }
419
420 /**
421  * Almost-equal for IEEE floats, using absolute difference method.
422  *
423  * \param max_diff: the maximum absolute difference.
424  */
425 MINLINE int compare_ff(float a, float b, const float max_diff)
426 {
427   return fabsf(a - b) <= max_diff;
428 }
429
430 /**
431  * Almost-equal for IEEE floats, using their integer representation
432  * (mixing ULP and absolute difference methods).
433  *
434  * \param max_diff: is the maximum absolute difference (allows to take care of the near-zero area,
435  * where relative difference methods cannot really work).
436  * \param max_ulps: is the 'maximum number of floats + 1'
437  * allowed between \a a and \a b to consider them equal.
438  *
439  * \see https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
440  */
441 MINLINE int compare_ff_relative(float a, float b, const float max_diff, const int max_ulps)
442 {
443   union {
444     float f;
445     int i;
446   } ua, ub;
447
448   BLI_assert(sizeof(float) == sizeof(int));
449   BLI_assert(max_ulps < (1 << 22));
450
451   if (fabsf(a - b) <= max_diff) {
452     return 1;
453   }
454
455   ua.f = a;
456   ub.f = b;
457
458   /* Important to compare sign from integers, since (-0.0f < 0) is false
459    * (though this shall not be an issue in common cases)... */
460   return ((ua.i < 0) != (ub.i < 0)) ? 0 : (abs(ua.i - ub.i) <= max_ulps) ? 1 : 0;
461 }
462
463 MINLINE float signf(float f)
464 {
465   return (f < 0.f) ? -1.f : 1.f;
466 }
467
468 MINLINE int signum_i_ex(float a, float eps)
469 {
470   if (a > eps) {
471     return 1;
472   }
473   if (a < -eps) {
474     return -1;
475   }
476   else {
477     return 0;
478   }
479 }
480
481 MINLINE int signum_i(float a)
482 {
483   if (a > 0.0f) {
484     return 1;
485   }
486   if (a < 0.0f) {
487     return -1;
488   }
489   else {
490     return 0;
491   }
492 }
493
494 /** Returns number of (base ten) *significant* digits of integer part of given float
495  * (negative in case of decimal-only floats, 0.01 returns -1 e.g.). */
496 MINLINE int integer_digits_f(const float f)
497 {
498   return (f == 0.0f) ? 0 : (int)floor(log10(fabs(f))) + 1;
499 }
500
501 /** Returns number of (base ten) *significant* digits of integer part of given double
502  * (negative in case of decimal-only floats, 0.01 returns -1 e.g.). */
503 MINLINE int integer_digits_d(const double d)
504 {
505   return (d == 0.0) ? 0 : (int)floor(log10(fabs(d))) + 1;
506 }
507
508 MINLINE int integer_digits_i(const int i)
509 {
510   return (int)log10((double)i) + 1;
511 }
512
513 /* Internal helpers for SSE2 implementation.
514  *
515  * NOTE: Are to be called ONLY from inside `#ifdef __SSE2__` !!!
516  */
517
518 #ifdef __SSE2__
519
520 /* Calculate initial guess for arg^exp based on float representation
521  * This method gives a constant bias, which can be easily compensated by
522  * multiplicating with bias_coeff.
523  * Gives better results for exponents near 1 (e. g. 4/5).
524  * exp = exponent, encoded as uint32_t
525  * e2coeff = 2^(127/exponent - 127) * bias_coeff^(1/exponent), encoded as
526  * uint32_t
527  *
528  * We hope that exp and e2coeff gets properly inlined
529  */
530 MALWAYS_INLINE __m128 _bli_math_fastpow(const int exp, const int e2coeff, const __m128 arg)
531 {
532   __m128 ret;
533   ret = _mm_mul_ps(arg, _mm_castsi128_ps(_mm_set1_epi32(e2coeff)));
534   ret = _mm_cvtepi32_ps(_mm_castps_si128(ret));
535   ret = _mm_mul_ps(ret, _mm_castsi128_ps(_mm_set1_epi32(exp)));
536   ret = _mm_castsi128_ps(_mm_cvtps_epi32(ret));
537   return ret;
538 }
539
540 /* Improve x ^ 1.0f/5.0f solution with Newton-Raphson method */
541 MALWAYS_INLINE __m128 _bli_math_improve_5throot_solution(const __m128 old_result, const __m128 x)
542 {
543   __m128 approx2 = _mm_mul_ps(old_result, old_result);
544   __m128 approx4 = _mm_mul_ps(approx2, approx2);
545   __m128 t = _mm_div_ps(x, approx4);
546   __m128 summ = _mm_add_ps(_mm_mul_ps(_mm_set1_ps(4.0f), old_result), t); /* fma */
547   return _mm_mul_ps(summ, _mm_set1_ps(1.0f / 5.0f));
548 }
549
550 /* Calculate powf(x, 2.4). Working domain: 1e-10 < x < 1e+10 */
551 MALWAYS_INLINE __m128 _bli_math_fastpow24(const __m128 arg)
552 {
553   /* max, avg and |avg| errors were calculated in gcc without FMA instructions
554    * The final precision should be better than powf in glibc */
555
556   /* Calculate x^4/5, coefficient 0.994 was constructed manually to minimize
557    * avg error.
558    */
559   /* 0x3F4CCCCD = 4/5 */
560   /* 0x4F55A7FB = 2^(127/(4/5) - 127) * 0.994^(1/(4/5)) */
561   /* error max = 0.17, avg = 0.0018, |avg| = 0.05 */
562   __m128 x = _bli_math_fastpow(0x3F4CCCCD, 0x4F55A7FB, arg);
563   __m128 arg2 = _mm_mul_ps(arg, arg);
564   __m128 arg4 = _mm_mul_ps(arg2, arg2);
565   /* error max = 0.018        avg = 0.0031    |avg| = 0.0031  */
566   x = _bli_math_improve_5throot_solution(x, arg4);
567   /* error max = 0.00021    avg = 1.6e-05    |avg| = 1.6e-05 */
568   x = _bli_math_improve_5throot_solution(x, arg4);
569   /* error max = 6.1e-07    avg = 5.2e-08    |avg| = 1.1e-07 */
570   x = _bli_math_improve_5throot_solution(x, arg4);
571   return _mm_mul_ps(x, _mm_mul_ps(x, x));
572 }
573
574 /* Calculate powf(x, 1.0f / 2.4) */
575 MALWAYS_INLINE __m128 _bli_math_fastpow512(const __m128 arg)
576 {
577   /* 5/12 is too small, so compute the 4th root of 20/12 instead.
578    * 20/12 = 5/3 = 1 + 2/3 = 2 - 1/3. 2/3 is a suitable argument for fastpow.
579    * weighting coefficient: a^-1/2 = 2 a; a = 2^-2/3
580    */
581   __m128 xf = _bli_math_fastpow(0x3f2aaaab, 0x5eb504f3, arg);
582   __m128 xover = _mm_mul_ps(arg, xf);
583   __m128 xfm1 = _mm_rsqrt_ps(xf);
584   __m128 x2 = _mm_mul_ps(arg, arg);
585   __m128 xunder = _mm_mul_ps(x2, xfm1);
586   /* sqrt2 * over + 2 * sqrt2 * under */
587   __m128 xavg = _mm_mul_ps(_mm_set1_ps(1.0f / (3.0f * 0.629960524947437f) * 0.999852f),
588                            _mm_add_ps(xover, xunder));
589   xavg = _mm_mul_ps(xavg, _mm_rsqrt_ps(xavg));
590   xavg = _mm_mul_ps(xavg, _mm_rsqrt_ps(xavg));
591   return xavg;
592 }
593
594 MALWAYS_INLINE __m128 _bli_math_blend_sse(const __m128 mask, const __m128 a, const __m128 b)
595 {
596   return _mm_or_ps(_mm_and_ps(mask, a), _mm_andnot_ps(mask, b));
597 }
598
599 #endif /* __SSE2__ */
600
601 /* Low level conversion functions */
602 MINLINE unsigned char unit_float_to_uchar_clamp(float val)
603 {
604   return (unsigned char)((
605       (val <= 0.0f) ? 0 : ((val > (1.0f - 0.5f / 255.0f)) ? 255 : ((255.0f * val) + 0.5f))));
606 }
607 #define unit_float_to_uchar_clamp(val) \
608   ((CHECK_TYPE_INLINE(val, float)), unit_float_to_uchar_clamp(val))
609
610 MINLINE unsigned short unit_float_to_ushort_clamp(float val)
611 {
612   return (unsigned short)((val >= 1.0f - 0.5f / 65535) ?
613                               65535 :
614                               (val <= 0.0f) ? 0 : (val * 65535.0f + 0.5f));
615 }
616 #define unit_float_to_ushort_clamp(val) \
617   ((CHECK_TYPE_INLINE(val, float)), unit_float_to_ushort_clamp(val))
618
619 MINLINE unsigned char unit_ushort_to_uchar(unsigned short val)
620 {
621   return (unsigned char)(((val) >= 65535 - 128) ? 255 : ((val) + 128) >> 8);
622 }
623 #define unit_ushort_to_uchar(val) \
624   ((CHECK_TYPE_INLINE(val, unsigned short)), unit_ushort_to_uchar(val))
625
626 #define unit_float_to_uchar_clamp_v3(v1, v2) \
627   { \
628     (v1)[0] = unit_float_to_uchar_clamp((v2[0])); \
629     (v1)[1] = unit_float_to_uchar_clamp((v2[1])); \
630     (v1)[2] = unit_float_to_uchar_clamp((v2[2])); \
631   } \
632   ((void)0)
633 #define unit_float_to_uchar_clamp_v4(v1, v2) \
634   { \
635     (v1)[0] = unit_float_to_uchar_clamp((v2[0])); \
636     (v1)[1] = unit_float_to_uchar_clamp((v2[1])); \
637     (v1)[2] = unit_float_to_uchar_clamp((v2[2])); \
638     (v1)[3] = unit_float_to_uchar_clamp((v2[3])); \
639   } \
640   ((void)0)
641
642 #endif /* __MATH_BASE_INLINE_C__ */