1 /**************************************************************************
3 * Copyright 2008 VMware, Inc.
6 * Permission is hereby granted, free of charge, to any person obtaining a
7 * copy of this software and associated documentation files (the
8 * "Software"), to deal in the Software without restriction, including
9 * without limitation the rights to use, copy, modify, merge, publish,
10 * distribute, sub license, and/or sell copies of the Software, and to
11 * permit persons to whom the Software is furnished to do so, subject to
12 * the following conditions:
14 * The above copyright notice and this permission notice (including the
15 * next paragraph) shall be included in all copies or substantial portions
18 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
19 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
20 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NON-INFRINGEMENT.
21 * IN NO EVENT SHALL VMWARE AND/OR ITS SUPPLIERS BE LIABLE FOR
22 * ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
23 * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
24 * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
26 **************************************************************************/
30 * Math utilities and approximations for common math functions.
31 * Reduced precision is usually acceptable in shaders...
33 * "fast" is used in the names of functions which are low-precision,
34 * or at least lower-precision than the normal C lib functions.
42 #include "pipe/p_compiler.h"
52 #include <strings.h> /* for ffs */
66 #define M_SQRT2 1.41421356237309504880
69 #define POW2_TABLE_SIZE_LOG2 9
70 #define POW2_TABLE_SIZE (1 << POW2_TABLE_SIZE_LOG2)
71 #define POW2_TABLE_OFFSET (POW2_TABLE_SIZE/2)
72 #define POW2_TABLE_SCALE ((float)(POW2_TABLE_SIZE/2))
73 extern float pow2_table[POW2_TABLE_SIZE];
77 * Initialize math module. This should be called before using any
78 * other functions in this module.
99 * Extract the IEEE float32 exponent.
102 util_get_float32_exponent(float x)
108 return ((f.ui >> 23) & 0xff) - 127;
113 * Fast version of 2^x
114 * Identity: exp2(a + b) = exp2(a) * exp2(b)
116 * Let fpart = x - ipart;
117 * So, exp2(x) = exp2(ipart) * exp2(fpart)
118 * Compute exp2(ipart) with i << ipart
119 * Compute exp2(fpart) with lookup table.
122 util_fast_exp2(float x)
129 return 3.402823466e+38f;
135 fpart = x - (float) ipart;
138 * epart.f = (float) (1 << ipart)
139 * but faster and without integer overflow for ipart > 31
141 epart.i = (ipart + 127 ) << 23;
143 mpart = pow2_table[POW2_TABLE_OFFSET + (int)(fpart * POW2_TABLE_SCALE)];
145 return epart.f * mpart;
150 * Fast approximation to exp(x).
153 util_fast_exp(float x)
155 const float k = 1.44269f; /* = log2(e) */
156 return util_fast_exp2(k * x);
160 #define LOG2_TABLE_SIZE_LOG2 16
161 #define LOG2_TABLE_SCALE (1 << LOG2_TABLE_SIZE_LOG2)
162 #define LOG2_TABLE_SIZE (LOG2_TABLE_SCALE + 1)
163 extern float log2_table[LOG2_TABLE_SIZE];
167 * Fast approximation to log2(x).
170 util_fast_log2(float x)
175 epart = (float)(((num.i & 0x7f800000) >> 23) - 127);
176 /* mpart = log2_table[mantissa*LOG2_TABLE_SCALE + 0.5] */
177 mpart = log2_table[((num.i & 0x007fffff) + (1 << (22 - LOG2_TABLE_SIZE_LOG2))) >> (23 - LOG2_TABLE_SIZE_LOG2)];
178 return epart + mpart;
183 * Fast approximation to x^y.
186 util_fast_pow(float x, float y)
188 return util_fast_exp2(util_fast_log2(x) * y);
191 /* Note that this counts zero as a power of two.
193 static inline boolean
194 util_is_power_of_two( unsigned v )
196 return (v & (v-1)) == 0;
201 * Floor(x), returned as int.
209 af = (3 << 22) + 0.5 + (double) f;
210 bf = (3 << 22) + 0.5 - (double) f;
211 u.f = (float) af; ai = u.i;
212 u.f = (float) bf; bi = u.i;
213 return (ai - bi) >> 1;
218 * Round float to nearest int.
223 #if defined(PIPE_CC_GCC) && defined(PIPE_ARCH_X86)
225 __asm__ ("fistpl %0" : "=m" (r) : "t" (f) : "st");
227 #elif defined(PIPE_CC_MSVC) && defined(PIPE_ARCH_X86)
236 return (int) (f + 0.5f);
238 return (int) (f - 0.5f);
244 * Approximate floating point comparison
246 static inline boolean
247 util_is_approx(float a, float b, float tol)
249 return fabsf(b - a) <= tol;
254 * util_is_X_inf_or_nan = test if x is NaN or +/- Inf
255 * util_is_X_nan = test if x is NaN
256 * util_X_inf_sign = return +1 for +Inf, -1 for -Inf, or 0 for not Inf
258 * NaN can be checked with x != x, however this fails with the fast math flag
265 static inline boolean
266 util_is_inf_or_nan(float x)
270 return (tmp.ui & 0x7f800000) == 0x7f800000;
274 static inline boolean
279 return (tmp.ui & 0x7fffffff) > 0x7f800000;
284 util_inf_sign(float x)
288 if ((tmp.ui & 0x7fffffff) != 0x7f800000) {
292 return (x < 0) ? -1 : 1;
299 static inline boolean
300 util_is_double_inf_or_nan(double x)
304 return (tmp.ui & 0x7ff0000000000000ULL) == 0x7ff0000000000000ULL;
308 static inline boolean
309 util_is_double_nan(double x)
313 return (tmp.ui & 0x7fffffffffffffffULL) > 0x7ff0000000000000ULL;
318 util_double_inf_sign(double x)
322 if ((tmp.ui & 0x7fffffffffffffffULL) != 0x7ff0000000000000ULL) {
326 return (x < 0) ? -1 : 1;
333 static inline boolean
334 util_is_half_inf_or_nan(int16_t x)
336 return (x & 0x7c00) == 0x7c00;
340 static inline boolean
341 util_is_half_nan(int16_t x)
343 return (x & 0x7fff) > 0x7c00;
348 util_half_inf_sign(int16_t x)
350 if ((x & 0x7fff) != 0x7c00) {
354 return (x < 0) ? -1 : 1;
359 * Find first bit set in word. Least significant bit is 1.
360 * Return 0 if no bits set.
363 #define FFS_DEFINED 1
365 #if defined(_MSC_VER) && (_M_IX86 || _M_AMD64 || _M_IA64)
367 unsigned long ffs( unsigned long u )
370 if (_BitScanForward(&i, u))
375 #elif defined(PIPE_CC_MSVC) && defined(PIPE_ARCH_X86)
377 unsigned ffs( unsigned u )
391 #elif defined(__MINGW32__) || defined(PIPE_OS_ANDROID) || \
392 defined(HAVE___BUILTIN_FFS)
393 #define ffs __builtin_ffs
396 #ifdef HAVE___BUILTIN_FFSLL
397 #define ffsll __builtin_ffsll
400 ffsll(long long int val)
404 bit = ffs((unsigned) (val & 0xffffffff));
408 bit = ffs((unsigned) (val >> 32));
416 #endif /* FFS_DEFINED */
419 * Find first bit set in long long. Least significant bit is 1.
420 * Return 0 if no bits set.
422 #ifndef FFSLL_DEFINED
423 #define FFSLL_DEFINED 1
425 #if defined(__MINGW32__) || defined(PIPE_OS_ANDROID) || \
426 defined(HAVE___BUILTIN_FFSLL)
427 #define ffsll __builtin_ffsll
430 #endif /* FFSLL_DEFINED */
433 * Find last bit set in a word. The least significant bit is 1.
434 * Return 0 if no bits are set.
436 static inline unsigned
437 util_last_bit(unsigned u)
439 #if defined(HAVE___BUILTIN_CLZ)
440 return u == 0 ? 0 : 32 - __builtin_clz(u);
452 * Find last bit set in a word. The least significant bit is 1.
453 * Return 0 if no bits are set.
455 static inline unsigned
456 util_last_bit64(uint64_t u)
458 #if defined(HAVE___BUILTIN_CLZLL)
459 return u == 0 ? 0 : 64 - __builtin_clzll(u);
471 * Find last bit in a word that does not match the sign bit. The least
472 * significant bit is 1.
473 * Return 0 if no bits are set.
475 static inline unsigned
476 util_last_bit_signed(int i)
479 return util_last_bit(i);
481 return util_last_bit(~(unsigned)i);
484 /* Destructively loop over all of the bits in a mask as in:
487 * int i = u_bit_scan(&mymask);
488 * ... process element i
493 u_bit_scan(unsigned *mask)
495 int i = ffs(*mask) - 1;
502 u_bit_scan64(uint64_t *mask)
504 int i = ffsll(*mask) - 1;
505 *mask &= ~(1llu << i);
510 /* For looping over a bitmask when you want to loop over consecutive bits
511 * manually, for example:
514 * int start, count, i;
516 * u_bit_scan_consecutive_range(&mask, &start, &count);
518 * for (i = 0; i < count; i++)
519 * ... process element (start+i)
523 u_bit_scan_consecutive_range(unsigned *mask, int *start, int *count)
525 if (*mask == 0xffffffff) {
531 *start = ffs(*mask) - 1;
532 *count = ffs(~(*mask >> *start)) - 1;
533 *mask &= ~(((1u << *count) - 1) << *start);
537 u_bit_scan_consecutive_range64(uint64_t *mask, int *start, int *count)
539 if (*mask == ~0llu) {
545 *start = ffsll(*mask) - 1;
546 *count = ffsll(~(*mask >> *start)) - 1;
547 *mask &= ~(((1llu << *count) - 1) << *start);
550 /* Returns a bitfield in which the first count bits starting at start are
553 static inline unsigned
554 u_bit_consecutive(unsigned start, unsigned count)
556 assert(start + count <= 32);
559 return ((1u << count) - 1) << start;
565 static inline unsigned
583 * Convert ubyte to float in [0, 1].
584 * XXX a 256-entry lookup table would be slightly faster.
587 ubyte_to_float(ubyte ub)
589 return (float) ub * (1.0f / 255.0f);
594 * Convert float in [0,1] to ubyte in [0,255] with clamping.
597 float_to_ubyte(float f)
605 else if (tmp.i >= 0x3f800000 /* 1.0f */) {
609 tmp.f = tmp.f * (255.0f/256.0f) + 32768.0f;
610 return (ubyte) tmp.i;
615 byte_to_float_tex(int8_t b)
617 return (b == -128) ? -1.0F : b * 1.0F / 127.0F;
621 float_to_byte_tex(float f)
623 return (int8_t) (127.0F * f);
629 static inline unsigned
630 util_logbase2(unsigned n)
632 #if defined(HAVE___BUILTIN_CLZ)
633 return ((sizeof(unsigned) * 8 - 1) - __builtin_clz(n | 1));
636 if (n >= 1<<16) { n >>= 16; pos += 16; }
637 if (n >= 1<< 8) { n >>= 8; pos += 8; }
638 if (n >= 1<< 4) { n >>= 4; pos += 4; }
639 if (n >= 1<< 2) { n >>= 2; pos += 2; }
640 if (n >= 1<< 1) { pos += 1; }
647 * Returns the smallest power of two >= x
649 static inline unsigned
650 util_next_power_of_two(unsigned x)
652 #if defined(HAVE___BUILTIN_CLZ)
656 return (1 << ((sizeof(unsigned) * 8) - __builtin_clz(x - 1)));
663 if (util_is_power_of_two(x))
667 val = (val >> 1) | val;
668 val = (val >> 2) | val;
669 val = (val >> 4) | val;
670 val = (val >> 8) | val;
671 val = (val >> 16) | val;
679 * Return number of bits set in n.
681 static inline unsigned
682 util_bitcount(unsigned n)
684 #if defined(HAVE___BUILTIN_POPCOUNT)
685 return __builtin_popcount(n);
687 /* K&R classic bitcount.
689 * For each iteration, clear the LSB from the bitfield.
690 * Requires only one iteration per set bit, instead of
691 * one iteration per bit less than highest set bit.
694 for (bits = 0; n; bits++) {
702 static inline unsigned
703 util_bitcount64(uint64_t n)
705 #ifdef HAVE___BUILTIN_POPCOUNTLL
706 return __builtin_popcountll(n);
708 return util_bitcount(n) + util_bitcount(n >> 32);
715 * Algorithm taken from:
716 * http://stackoverflow.com/questions/9144800/c-reverse-bits-in-unsigned-integer
718 static inline unsigned
719 util_bitreverse(unsigned n)
721 n = ((n >> 1) & 0x55555555u) | ((n & 0x55555555u) << 1);
722 n = ((n >> 2) & 0x33333333u) | ((n & 0x33333333u) << 2);
723 n = ((n >> 4) & 0x0f0f0f0fu) | ((n & 0x0f0f0f0fu) << 4);
724 n = ((n >> 8) & 0x00ff00ffu) | ((n & 0x00ff00ffu) << 8);
725 n = ((n >> 16) & 0xffffu) | ((n & 0xffffu) << 16);
730 * Convert from little endian to CPU byte order.
733 #ifdef PIPE_ARCH_BIG_ENDIAN
734 #define util_le64_to_cpu(x) util_bswap64(x)
735 #define util_le32_to_cpu(x) util_bswap32(x)
736 #define util_le16_to_cpu(x) util_bswap16(x)
738 #define util_le64_to_cpu(x) (x)
739 #define util_le32_to_cpu(x) (x)
740 #define util_le16_to_cpu(x) (x)
743 #define util_cpu_to_le64(x) util_le64_to_cpu(x)
744 #define util_cpu_to_le32(x) util_le32_to_cpu(x)
745 #define util_cpu_to_le16(x) util_le16_to_cpu(x)
748 * Reverse byte order of a 32 bit word.
750 static inline uint32_t
751 util_bswap32(uint32_t n)
753 #if defined(HAVE___BUILTIN_BSWAP32)
754 return __builtin_bswap32(n);
757 ((n >> 8) & 0x0000ff00) |
758 ((n << 8) & 0x00ff0000) |
764 * Reverse byte order of a 64bit word.
766 static inline uint64_t
767 util_bswap64(uint64_t n)
769 #if defined(HAVE___BUILTIN_BSWAP64)
770 return __builtin_bswap64(n);
772 return ((uint64_t)util_bswap32((uint32_t)n) << 32) |
773 util_bswap32((n >> 32));
779 * Reverse byte order of a 16 bit word.
781 static inline uint16_t
782 util_bswap16(uint16_t n)
789 util_memcpy_cpu_to_le32(void * restrict dest, const void * restrict src, size_t n)
791 #ifdef PIPE_ARCH_BIG_ENDIAN
795 for (i = 0, e = n / 4; i < e; i++) {
796 uint32_t * restrict d = (uint32_t* restrict)dest;
797 const uint32_t * restrict s = (const uint32_t* restrict)src;
798 d[i] = util_bswap32(s[i]);
802 return memcpy(dest, src, n);
807 * Clamp X to [MIN, MAX].
808 * This is a macro to allow float, int, uint, etc. types.
810 #define CLAMP( X, MIN, MAX ) ( (X)<(MIN) ? (MIN) : ((X)>(MAX) ? (MAX) : (X)) )
812 #define MIN2( A, B ) ( (A)<(B) ? (A) : (B) )
813 #define MAX2( A, B ) ( (A)>(B) ? (A) : (B) )
815 #define MIN3( A, B, C ) ((A) < (B) ? MIN2(A, C) : MIN2(B, C))
816 #define MAX3( A, B, C ) ((A) > (B) ? MAX2(A, C) : MAX2(B, C))
818 #define MIN4( A, B, C, D ) ((A) < (B) ? MIN3(A, C, D) : MIN3(B, C, D))
819 #define MAX4( A, B, C, D ) ((A) > (B) ? MAX3(A, C, D) : MAX3(B, C, D))
823 * Align a value, only works pot alignemnts.
826 align(int value, int alignment)
828 return (value + alignment - 1) & ~(alignment - 1);
831 static inline uint64_t
832 align64(uint64_t value, unsigned alignment)
834 return (value + alignment - 1) & ~((uint64_t)alignment - 1);
838 * Works like align but on npot alignments.
841 util_align_npot(size_t value, size_t alignment)
843 if (value % alignment)
844 return value + (alignment - (value % alignment));
848 static inline unsigned
849 u_minify(unsigned value, unsigned levels)
851 return MAX2(1, value >> levels);
855 #define COPY_4V( DST, SRC ) \
857 (DST)[0] = (SRC)[0]; \
858 (DST)[1] = (SRC)[1]; \
859 (DST)[2] = (SRC)[2]; \
860 (DST)[3] = (SRC)[3]; \
866 #define COPY_4FV( DST, SRC ) COPY_4V(DST, SRC)
871 #define ASSIGN_4V( DST, V0, V1, V2, V3 ) \
881 static inline uint32_t
882 util_unsigned_fixed(float value, unsigned frac_bits)
884 return value < 0 ? 0 : (uint32_t)(value * (1<<frac_bits));
887 static inline int32_t
888 util_signed_fixed(float value, unsigned frac_bits)
890 return (int32_t)(value * (1<<frac_bits));
894 util_fpstate_get(void);
896 util_fpstate_set_denorms_to_zero(unsigned current_fpstate);
898 util_fpstate_set(unsigned fpstate);
906 #endif /* U_MATH_H */