common/alnumeric.h

   1 #ifndef AL_NUMERIC_H
   2 #define AL_NUMERIC_H
   3
   4 #include <cstddef>
   5 #include <cstdint>
   6 #ifdef HAVE_INTRIN_H
   7 #include <intrin.h>
   8 #endif
   9 #ifdef HAVE_SSE_INTRINSICS
  10 #include <xmmintrin.h>
  11 #endif
  12
  13 #include "opthelpers.h"
  14
  15
  16 inline constexpr int64_t operator "" _i64(unsigned long long int n) noexcept { return static_cast<int64_t>(n); }
  17 inline constexpr uint64_t operator "" _u64(unsigned long long int n) noexcept { return static_cast<uint64_t>(n); }
  18
  19
  20 constexpr inline float minf(float a, float b) noexcept
  21 { return ((a > b) ? b : a); }
  22 constexpr inline float maxf(float a, float b) noexcept
  23 { return ((a > b) ? a : b); }
  24 constexpr inline float clampf(float val, float min, float max) noexcept
  25 { return minf(max, maxf(min, val)); }
  26
  27 constexpr inline double mind(double a, double b) noexcept
  28 { return ((a > b) ? b : a); }
  29 constexpr inline double maxd(double a, double b) noexcept
  30 { return ((a > b) ? a : b); }
  31 constexpr inline double clampd(double val, double min, double max) noexcept
  32 { return mind(max, maxd(min, val)); }
  33
  34 constexpr inline unsigned int minu(unsigned int a, unsigned int b) noexcept
  35 { return ((a > b) ? b : a); }
  36 constexpr inline unsigned int maxu(unsigned int a, unsigned int b) noexcept
  37 { return ((a > b) ? a : b); }
  38 constexpr inline unsigned int clampu(unsigned int val, unsigned int min, unsigned int max) noexcept
  39 { return minu(max, maxu(min, val)); }
  40
  41 constexpr inline int mini(int a, int b) noexcept
  42 { return ((a > b) ? b : a); }
  43 constexpr inline int maxi(int a, int b) noexcept
  44 { return ((a > b) ? a : b); }
  45 constexpr inline int clampi(int val, int min, int max) noexcept
  46 { return mini(max, maxi(min, val)); }
  47
  48 constexpr inline int64_t mini64(int64_t a, int64_t b) noexcept
  49 { return ((a > b) ? b : a); }
  50 constexpr inline int64_t maxi64(int64_t a, int64_t b) noexcept
  51 { return ((a > b) ? a : b); }
  52 constexpr inline int64_t clampi64(int64_t val, int64_t min, int64_t max) noexcept
  53 { return mini64(max, maxi64(min, val)); }
  54
  55 constexpr inline uint64_t minu64(uint64_t a, uint64_t b) noexcept
  56 { return ((a > b) ? b : a); }
  57 constexpr inline uint64_t maxu64(uint64_t a, uint64_t b) noexcept
  58 { return ((a > b) ? a : b); }
  59 constexpr inline uint64_t clampu64(uint64_t val, uint64_t min, uint64_t max) noexcept
  60 { return minu64(max, maxu64(min, val)); }
  61
  62 constexpr inline size_t minz(size_t a, size_t b) noexcept
  63 { return ((a > b) ? b : a); }
  64 constexpr inline size_t maxz(size_t a, size_t b) noexcept
  65 { return ((a > b) ? a : b); }
  66 constexpr inline size_t clampz(size_t val, size_t min, size_t max) noexcept
  67 { return minz(max, maxz(min, val)); }
  68
  69
  70 /** Find the next power-of-2 for non-power-of-2 numbers. */
  71 inline uint32_t NextPowerOf2(uint32_t value) noexcept
  72 {
  73     if(value > 0)
  74     {
  75         value--;
  76         value |= value>>1;
  77         value |= value>>2;
  78         value |= value>>4;
  79         value |= value>>8;
  80         value |= value>>16;
  81     }
  82     return value+1;
  83 }
  84
  85 /** Round up a value to the next multiple. */
  86 inline size_t RoundUp(size_t value, size_t r) noexcept
  87 {
  88     value += r-1;
  89     return value - (value%r);
  90 }
  91
  92
  93 /* Define CTZ macros (count trailing zeros), and POPCNT macros (population
  94  * count/count 1 bits), for 32- and 64-bit integers. The CTZ macros' results
  95  * are *UNDEFINED* if the value is 0.
  96  */
  97 #ifdef __GNUC__
  98
  99 #define POPCNT32 __builtin_popcount
 100 #define CTZ32 __builtin_ctz
 101 #if SIZEOF_LONG == 8
 102 #define POPCNT64 __builtin_popcountl
 103 #define CTZ64 __builtin_ctzl
 104 #else
 105 #define POPCNT64 __builtin_popcountll
 106 #define CTZ64 __builtin_ctzll
 107 #endif
 108
 109 #else
 110
 111 /* There be black magics here. The popcnt method is derived from
 112  * https://graphics.stanford.edu/~seander/bithacks.html#CountBitsSetParallel
 113  * while the ctz-utilizing-popcnt algorithm is shown here
 114  * http://www.hackersdelight.org/hdcodetxt/ntz.c.txt
 115  * as the ntz2 variant. These likely aren't the most efficient methods, but
 116  * they're good enough if the GCC built-ins aren't available.
 117  */
 118 inline int fallback_popcnt32(uint32_t v)
 119 {
 120     v = v - ((v >> 1) & 0x55555555u);
 121     v = (v & 0x33333333u) + ((v >> 2) & 0x33333333u);
 122     v = (v + (v >> 4)) & 0x0f0f0f0fu;
 123     return (int)((v * 0x01010101u) >> 24);
 124 }
 125 #define POPCNT32 fallback_popcnt32
 126 inline int fallback_popcnt64(uint64_t v)
 127 {
 128     v = v - ((v >> 1) & 0x5555555555555555_u64);
 129     v = (v & 0x3333333333333333_u64) + ((v >> 2) & 0x3333333333333333_u64);
 130     v = (v + (v >> 4)) & 0x0f0f0f0f0f0f0f0f_u64;
 131     return (int)((v * 0x0101010101010101_u64) >> 56);
 132 }
 133 #define POPCNT64 fallback_popcnt64
 134
 135 #if defined(HAVE_BITSCANFORWARD64_INTRINSIC)
 136
 137 inline int msvc64_ctz32(uint32_t v)
 138 {
 139     unsigned long idx = 32;
 140     _BitScanForward(&idx, v);
 141     return (int)idx;
 142 }
 143 #define CTZ32 msvc64_ctz32
 144 inline int msvc64_ctz64(uint64_t v)
 145 {
 146     unsigned long idx = 64;
 147     _BitScanForward64(&idx, v);
 148     return (int)idx;
 149 }
 150 #define CTZ64 msvc64_ctz64
 151
 152 #elif defined(HAVE_BITSCANFORWARD_INTRINSIC)
 153
 154 inline int msvc_ctz32(uint32_t v)
 155 {
 156     unsigned long idx = 32;
 157     _BitScanForward(&idx, v);
 158     return (int)idx;
 159 }
 160 #define CTZ32 msvc_ctz32
 161 inline int msvc_ctz64(uint64_t v)
 162 {
 163     unsigned long idx = 64;
 164     if(!_BitScanForward(&idx, (uint32_t)(v&0xffffffff)))
 165     {
 166         if(_BitScanForward(&idx, (uint32_t)(v>>32)))
 167             idx += 32;
 168     }
 169     return (int)idx;
 170 }
 171 #define CTZ64 msvc_ctz64
 172
 173 #else
 174
 175 inline int fallback_ctz32(uint32_t value)
 176 { return POPCNT32(~value & (value - 1)); }
 177 #define CTZ32 fallback_ctz32
 178 inline int fallback_ctz64(uint64_t value)
 179 { return POPCNT64(~value & (value - 1)); }
 180 #define CTZ64 fallback_ctz64
 181
 182 #endif
 183 #endif
 184
 185
 186 /**
 187  * Fast float-to-int conversion. No particular rounding mode is assumed; the
 188  * IEEE-754 default is round-to-nearest with ties-to-even, though an app could
 189  * change it on its own threads. On some systems, a truncating conversion may
 190  * always be the fastest method.
 191  */
 192 inline int fastf2i(float f) noexcept
 193 {
 194 #if defined(HAVE_SSE_INTRINSICS)
 195     return _mm_cvt_ss2si(_mm_set_ss(f));
 196
 197 #elif defined(_MSC_VER) && defined(_M_IX86_FP)
 198
 199     int i;
 200     __asm fld f
 201     __asm fistp i
 202     return i;
 203
 204 #elif (defined(__GNUC__) || defined(__clang__)) && (defined(__i386__) || defined(__x86_64__))
 205
 206     int i;
 207 #ifdef __SSE_MATH__
 208     __asm__("cvtss2si %1, %0" : "=r"(i) : "x"(f));
 209 #else
 210     __asm__ __volatile__("fistpl %0" : "=m"(i) : "t"(f) : "st");
 211 #endif
 212     return i;
 213
 214 #else
 215
 216     return static_cast<int>(f);
 217 #endif
 218 }
 219 inline unsigned int fastf2u(float f) noexcept
 220 { return static_cast<unsigned int>(fastf2i(f)); }
 221
 222 /** Converts float-to-int using standard behavior (truncation). */
 223 inline int float2int(float f) noexcept
 224 {
 225 #if defined(HAVE_SSE_INTRINSICS)
 226     return _mm_cvtt_ss2si(_mm_set_ss(f));
 227
 228 #elif ((defined(__GNUC__) || defined(__clang__)) && (defined(__i386__) || defined(__x86_64__)) && \
 229        !defined(__SSE_MATH__)) || (defined(_MSC_VER) && defined(_M_IX86_FP) && _M_IX86_FP == 0)
 230     int sign, shift, mant;
 231     union {
 232         float f;
 233         int i;
 234     } conv;
 235
 236     conv.f = f;
 237     sign = (conv.i>>31) | 1;
 238     shift = ((conv.i>>23)&0xff) - (127+23);
 239
 240     /* Over/underflow */
 241     if UNLIKELY(shift >= 31 || shift < -23)
 242         return 0;
 243
 244     mant = (conv.i&0x7fffff) | 0x800000;
 245     if LIKELY(shift < 0)
 246         return (mant >> -shift) * sign;
 247     return (mant << shift) * sign;
 248
 249 #else
 250
 251     return static_cast<int>(f);
 252 #endif
 253 }
 254 inline unsigned int float2uint(float f) noexcept
 255 { return static_cast<unsigned int>(float2int(f)); }
 256
 257 /** Converts double-to-int using standard behavior (truncation). */
 258 inline int double2int(double d) noexcept
 259 {
 260 #if defined(HAVE_SSE_INTRINSICS)
 261     return _mm_cvttsd_si32(_mm_set_sd(d));
 262
 263 #elif ((defined(__GNUC__) || defined(__clang__)) && (defined(__i386__) || defined(__x86_64__)) && \
 264        !defined(__SSE2_MATH__)) || (defined(_MSC_VER) && defined(_M_IX86_FP) && _M_IX86_FP < 2)
 265
 266     int sign, shift;
 267     int64_t mant;
 268     union {
 269         double d;
 270         int64_t i64;
 271     } conv;
 272
 273     conv.d = d;
 274     sign = (conv.i64 >> 63) | 1;
 275     shift = ((conv.i64 >> 52) & 0x7ff) - (1023 + 52);
 276
 277     /* Over/underflow */
 278     if UNLIKELY(shift >= 63 || shift < -52)
 279         return 0;
 280
 281     mant = (conv.i64 & 0xfffffffffffff_i64) | 0x10000000000000_i64;
 282     if LIKELY(shift < 0)
 283         return (int)(mant >> -shift) * sign;
 284     return (int)(mant << shift) * sign;
 285
 286 #else
 287
 288     return static_cast<int>(d);
 289 #endif
 290 }
 291
 292 /**
 293  * Rounds a float to the nearest integral value, according to the current
 294  * rounding mode. This is essentially an inlined version of rintf, although
 295  * makes fewer promises (e.g. -0 or -0.25 rounded to 0 may result in +0).
 296  */
 297 inline float fast_roundf(float f) noexcept
 298 {
 299 #if (defined(__GNUC__) || defined(__clang__)) && (defined(__i386__) || defined(__x86_64__)) && \
 300     !defined(__SSE_MATH__)
 301
 302     float out;
 303     __asm__ __volatile__("frndint" : "=t"(out) : "0"(f));
 304     return out;
 305
 306 #else
 307
 308     /* Integral limit, where sub-integral precision is not available for
 309      * floats.
 310      */
 311     static constexpr float ilim[2] = {
 312          8388608.0f /*  0x1.0p+23 */,
 313         -8388608.0f /* -0x1.0p+23 */
 314     };
 315     unsigned int sign, expo;
 316     union {
 317         float f;
 318         unsigned int i;
 319     } conv;
 320
 321     conv.f = f;
 322     sign = (conv.i>>31)&0x01;
 323     expo = (conv.i>>23)&0xff;
 324
 325     if UNLIKELY(expo >= 150/*+23*/)
 326     {
 327         /* An exponent (base-2) of 23 or higher is incapable of sub-integral
 328          * precision, so it's already an integral value. We don't need to worry
 329          * about infinity or NaN here.
 330          */
 331         return f;
 332     }
 333     /* Adding the integral limit to the value (with a matching sign) forces a
 334      * result that has no sub-integral precision, and is consequently forced to
 335      * round to an integral value. Removing the integral limit then restores
 336      * the initial value rounded to the integral. The compiler should not
 337      * optimize this out because of non-associative rules on floating-point
 338      * math (as long as you don't use -fassociative-math,
 339      * -funsafe-math-optimizations, -ffast-math, or -Ofast, in which case this
 340      * may break).
 341      */
 342     f += ilim[sign];
 343     return f - ilim[sign];
 344 #endif
 345 }
 346
 347 #endif /* AL_NUMERIC_H */