Source/WTF/wtf/MathExtras.h

   1 /*
   2  * Copyright (C) 2006-2018 Apple Inc. All rights reserved.
   3  *
   4  * Redistribution and use in source and binary forms, with or without
   5  * modification, are permitted provided that the following conditions
   6  * are met:
   7  * 1. Redistributions of source code must retain the above copyright
   8  *    notice, this list of conditions and the following disclaimer.
   9  * 2. Redistributions in binary form must reproduce the above copyright
  10  *    notice, this list of conditions and the following disclaimer in the
  11  *    documentation and/or other materials provided with the distribution.
  12  *
  13  * THIS SOFTWARE IS PROVIDED BY APPLE INC. ``AS IS'' AND ANY
  14  * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  15  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
  16  * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL APPLE INC. OR
  17  * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
  18  * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
  19  * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
  20  * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
  21  * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  22  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
  23  * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  24  */
  25
  26 #ifndef WTF_MathExtras_h
  27 #define WTF_MathExtras_h
  28
  29 #include <algorithm>
  30 #include <cmath>
  31 #include <float.h>
  32 #include <limits>
  33 #include <stdint.h>
  34 #include <stdlib.h>
  35 #include <wtf/StdLibExtras.h>
  36
  37 #if OS(OPENBSD)
  38 #include <sys/types.h>
  39 #include <machine/ieee.h>
  40 #endif
  41
  42 #ifndef M_PI
  43 const double piDouble = 3.14159265358979323846;
  44 const float piFloat = 3.14159265358979323846f;
  45 #else
  46 const double piDouble = M_PI;
  47 const float piFloat = static_cast<float>(M_PI);
  48 #endif
  49
  50 #ifndef M_PI_2
  51 const double piOverTwoDouble = 1.57079632679489661923;
  52 const float piOverTwoFloat = 1.57079632679489661923f;
  53 #else
  54 const double piOverTwoDouble = M_PI_2;
  55 const float piOverTwoFloat = static_cast<float>(M_PI_2);
  56 #endif
  57
  58 #ifndef M_PI_4
  59 const double piOverFourDouble = 0.785398163397448309616;
  60 const float piOverFourFloat = 0.785398163397448309616f;
  61 #else
  62 const double piOverFourDouble = M_PI_4;
  63 const float piOverFourFloat = static_cast<float>(M_PI_4);
  64 #endif
  65
  66 #ifndef M_SQRT2
  67 const double sqrtOfTwoDouble = 1.41421356237309504880;
  68 const float sqrtOfTwoFloat = 1.41421356237309504880f;
  69 #else
  70 const double sqrtOfTwoDouble = M_SQRT2;
  71 const float sqrtOfTwoFloat = static_cast<float>(M_SQRT2);
  72 #endif
  73
  74 #if COMPILER(MSVC)
  75
  76 // Work around a bug in Win, where atan2(+-infinity, +-infinity) yields NaN instead of specific values.
  77 extern "C" inline double wtf_atan2(double x, double y)
  78 {
  79     double posInf = std::numeric_limits<double>::infinity();
  80     double negInf = -std::numeric_limits<double>::infinity();
  81     double nan = std::numeric_limits<double>::quiet_NaN();
  82
  83     double result = nan;
  84
  85     if (x == posInf && y == posInf)
  86         result = piOverFourDouble;
  87     else if (x == posInf && y == negInf)
  88         result = 3 * piOverFourDouble;
  89     else if (x == negInf && y == posInf)
  90         result = -piOverFourDouble;
  91     else if (x == negInf && y == negInf)
  92         result = -3 * piOverFourDouble;
  93     else
  94         result = ::atan2(x, y);
  95
  96     return result;
  97 }
  98
  99 #define atan2(x, y) wtf_atan2(x, y)
 100
 101 #endif // COMPILER(MSVC)
 102
 103 inline double deg2rad(double d)  { return d * piDouble / 180.0; }
 104 inline double rad2deg(double r)  { return r * 180.0 / piDouble; }
 105 inline double deg2grad(double d) { return d * 400.0 / 360.0; }
 106 inline double grad2deg(double g) { return g * 360.0 / 400.0; }
 107 inline double turn2deg(double t) { return t * 360.0; }
 108 inline double deg2turn(double d) { return d / 360.0; }
 109 inline double rad2grad(double r) { return r * 200.0 / piDouble; }
 110 inline double grad2rad(double g) { return g * piDouble / 200.0; }
 111
 112 inline float deg2rad(float d)  { return d * piFloat / 180.0f; }
 113 inline float rad2deg(float r)  { return r * 180.0f / piFloat; }
 114 inline float deg2grad(float d) { return d * 400.0f / 360.0f; }
 115 inline float grad2deg(float g) { return g * 360.0f / 400.0f; }
 116 inline float turn2deg(float t) { return t * 360.0f; }
 117 inline float deg2turn(float d) { return d / 360.0f; }
 118 inline float rad2grad(float r) { return r * 200.0f / piFloat; }
 119 inline float grad2rad(float g) { return g * piFloat / 200.0f; }
 120
 121 // std::numeric_limits<T>::min() returns the smallest positive value for floating point types
 122 template<typename T> constexpr inline T defaultMinimumForClamp() { return std::numeric_limits<T>::min(); }
 123 template<> constexpr inline float defaultMinimumForClamp() { return -std::numeric_limits<float>::max(); }
 124 template<> constexpr inline double defaultMinimumForClamp() { return -std::numeric_limits<double>::max(); }
 125 template<typename T> constexpr inline T defaultMaximumForClamp() { return std::numeric_limits<T>::max(); }
 126
 127 template<typename T> inline T clampTo(double value, T min = defaultMinimumForClamp<T>(), T max = defaultMaximumForClamp<T>())
 128 {
 129     if (value >= static_cast<double>(max))
 130         return max;
 131     if (value <= static_cast<double>(min))
 132         return min;
 133     return static_cast<T>(value);
 134 }
 135 template<> inline long long int clampTo(double, long long int, long long int); // clampTo does not support long long ints.
 136
 137 inline int clampToInteger(double value)
 138 {
 139     return clampTo<int>(value);
 140 }
 141
 142 inline unsigned clampToUnsigned(double value)
 143 {
 144     return clampTo<unsigned>(value);
 145 }
 146
 147 inline float clampToFloat(double value)
 148 {
 149     return clampTo<float>(value);
 150 }
 151
 152 inline int clampToPositiveInteger(double value)
 153 {
 154     return clampTo<int>(value, 0);
 155 }
 156
 157 inline int clampToInteger(float value)
 158 {
 159     return clampTo<int>(value);
 160 }
 161
 162 template<typename T>
 163 inline int clampToInteger(T x)
 164 {
 165     static_assert(std::numeric_limits<T>::is_integer, "T must be an integer.");
 166
 167     const T intMax = static_cast<unsigned>(std::numeric_limits<int>::max());
 168
 169     if (x >= intMax)
 170         return std::numeric_limits<int>::max();
 171     return static_cast<int>(x);
 172 }
 173
 174 // Explicitly accept 64bit result when clamping double value.
 175 // Keep in mind that double can only represent 53bit integer precisely.
 176 template<typename T> constexpr inline T clampToAccepting64(double value, T min = defaultMinimumForClamp<T>(), T max = defaultMaximumForClamp<T>())
 177 {
 178     return (value >= static_cast<double>(max)) ? max : ((value <= static_cast<double>(min)) ? min : static_cast<T>(value));
 179 }
 180
 181 inline bool isWithinIntRange(float x)
 182 {
 183     return x > static_cast<float>(std::numeric_limits<int>::min()) && x < static_cast<float>(std::numeric_limits<int>::max());
 184 }
 185
 186 inline float normalizedFloat(float value)
 187 {
 188     if (value > 0 && value < std::numeric_limits<float>::min())
 189         return std::numeric_limits<float>::min();
 190     if (value < 0 && value > -std::numeric_limits<float>::min())
 191         return -std::numeric_limits<float>::min();
 192     return value;
 193 }
 194
 195 template<typename T> constexpr bool hasOneBitSet(T value)
 196 {
 197     return !((value - 1) & value) && value;
 198 }
 199
 200 template<typename T> constexpr bool hasZeroOrOneBitsSet(T value)
 201 {
 202     return !((value - 1) & value);
 203 }
 204
 205 template<typename T> constexpr bool hasTwoOrMoreBitsSet(T value)
 206 {
 207     return !hasZeroOrOneBitsSet(value);
 208 }
 209
 210 template <typename T> inline unsigned getLSBSet(T value)
 211 {
 212     typedef typename std::make_unsigned<T>::type UnsignedT;
 213     unsigned result = 0;
 214
 215     UnsignedT unsignedValue = static_cast<UnsignedT>(value);
 216     while (unsignedValue >>= 1)
 217         ++result;
 218
 219     return result;
 220 }
 221
 222 template<typename T> inline T divideRoundedUp(T a, T b)
 223 {
 224     return (a + b - 1) / b;
 225 }
 226
 227 template<typename T> inline T timesThreePlusOneDividedByTwo(T value)
 228 {
 229     // Mathematically equivalent to:
 230     //   (value * 3 + 1) / 2;
 231     // or:
 232     //   (unsigned)ceil(value * 1.5));
 233     // This form is not prone to internal overflow.
 234     return value + (value >> 1) + (value & 1);
 235 }
 236
 237 template<typename T> inline bool isNotZeroAndOrdered(T value)
 238 {
 239     return value > 0.0 || value < 0.0;
 240 }
 241
 242 template<typename T> inline bool isZeroOrUnordered(T value)
 243 {
 244     return !isNotZeroAndOrdered(value);
 245 }
 246
 247 template<typename T> inline bool isGreaterThanNonZeroPowerOfTwo(T value, unsigned power)
 248 {
 249     // The crazy way of testing of index >= 2 ** power
 250     // (where I use ** to denote pow()).
 251     return !!((value >> 1) >> (power - 1));
 252 }
 253
 254 template<typename T> constexpr inline bool isLessThan(const T& a, const T& b) { return a < b; }
 255 template<typename T> constexpr inline bool isLessThanEqual(const T& a, const T& b) { return a <= b; }
 256 template<typename T> constexpr inline bool isGreaterThan(const T& a, const T& b) { return a > b; }
 257 template<typename T> constexpr inline bool isGreaterThanEqual(const T& a, const T& b) { return a >= b; }
 258
 259 #ifndef UINT64_C
 260 #if COMPILER(MSVC)
 261 #define UINT64_C(c) c ## ui64
 262 #else
 263 #define UINT64_C(c) c ## ull
 264 #endif
 265 #endif
 266
 267 #if COMPILER(MINGW64) && (!defined(__MINGW64_VERSION_RC) || __MINGW64_VERSION_RC < 1)
 268 inline double wtf_pow(double x, double y)
 269 {
 270     // MinGW-w64 has a custom implementation for pow.
 271     // This handles certain special cases that are different.
 272     if ((x == 0.0 || std::isinf(x)) && std::isfinite(y)) {
 273         double f;
 274         if (modf(y, &f) != 0.0)
 275             return ((x == 0.0) ^ (y > 0.0)) ? std::numeric_limits<double>::infinity() : 0.0;
 276     }
 277
 278     if (x == 2.0) {
 279         int yInt = static_cast<int>(y);
 280         if (y == yInt)
 281             return ldexp(1.0, yInt);
 282     }
 283
 284     return pow(x, y);
 285 }
 286 #define pow(x, y) wtf_pow(x, y)
 287 #endif // COMPILER(MINGW64) && (!defined(__MINGW64_VERSION_RC) || __MINGW64_VERSION_RC < 1)
 288
 289
 290 // decompose 'number' to its sign, exponent, and mantissa components.
 291 // The result is interpreted as:
 292 //     (sign ? -1 : 1) * pow(2, exponent) * (mantissa / (1 << 52))
 293 inline void decomposeDouble(double number, bool& sign, int32_t& exponent, uint64_t& mantissa)
 294 {
 295     ASSERT(std::isfinite(number));
 296
 297     sign = std::signbit(number);
 298
 299     uint64_t bits = WTF::bitwise_cast<uint64_t>(number);
 300     exponent = (static_cast<int32_t>(bits >> 52) & 0x7ff) - 0x3ff;
 301     mantissa = bits & 0xFFFFFFFFFFFFFull;
 302
 303     // Check for zero/denormal values; if so, adjust the exponent,
 304     // if not insert the implicit, omitted leading 1 bit.
 305     if (exponent == -0x3ff)
 306         exponent = mantissa ? -0x3fe : 0;
 307     else
 308         mantissa |= 0x10000000000000ull;
 309 }
 310
 311 // Calculate d % 2^{64}.
 312 inline void doubleToInteger(double d, unsigned long long& value)
 313 {
 314     if (std::isnan(d) || std::isinf(d))
 315         value = 0;
 316     else {
 317         // -2^{64} < fmodValue < 2^{64}.
 318         double fmodValue = fmod(trunc(d), std::numeric_limits<unsigned long long>::max() + 1.0);
 319         if (fmodValue >= 0) {
 320             // 0 <= fmodValue < 2^{64}.
 321             // 0 <= value < 2^{64}. This cast causes no loss.
 322             value = static_cast<unsigned long long>(fmodValue);
 323         } else {
 324             // -2^{64} < fmodValue < 0.
 325             // 0 < fmodValueInUnsignedLongLong < 2^{64}. This cast causes no loss.
 326             unsigned long long fmodValueInUnsignedLongLong = static_cast<unsigned long long>(-fmodValue);
 327             // -1 < (std::numeric_limits<unsigned long long>::max() - fmodValueInUnsignedLongLong) < 2^{64} - 1.
 328             // 0 < value < 2^{64}.
 329             value = std::numeric_limits<unsigned long long>::max() - fmodValueInUnsignedLongLong + 1;
 330         }
 331     }
 332 }
 333
 334 namespace WTF {
 335
 336 // From http://graphics.stanford.edu/~seander/bithacks.html#RoundUpPowerOf2
 337 inline constexpr uint32_t roundUpToPowerOfTwo(uint32_t v)
 338 {
 339     v--;
 340     v |= v >> 1;
 341     v |= v >> 2;
 342     v |= v >> 4;
 343     v |= v >> 8;
 344     v |= v >> 16;
 345     v++;
 346     return v;
 347 }
 348
 349 inline constexpr unsigned maskForSize(unsigned size)
 350 {
 351     if (!size)
 352         return 0;
 353     return roundUpToPowerOfTwo(size) - 1;
 354 }
 355
 356 inline unsigned fastLog2(unsigned i)
 357 {
 358     unsigned log2 = 0;
 359     if (i & (i - 1))
 360         log2 += 1;
 361     if (i >> 16)
 362         log2 += 16, i >>= 16;
 363     if (i >> 8)
 364         log2 += 8, i >>= 8;
 365     if (i >> 4)
 366         log2 += 4, i >>= 4;
 367     if (i >> 2)
 368         log2 += 2, i >>= 2;
 369     if (i >> 1)
 370         log2 += 1;
 371     return log2;
 372 }
 373
 374 inline unsigned fastLog2(uint64_t value)
 375 {
 376     unsigned high = static_cast<unsigned>(value >> 32);
 377     if (high)
 378         return fastLog2(high) + 32;
 379     return fastLog2(static_cast<unsigned>(value));
 380 }
 381
 382 template <typename T>
 383 inline typename std::enable_if<std::is_floating_point<T>::value, T>::type safeFPDivision(T u, T v)
 384 {
 385     // Protect against overflow / underflow.
 386     if (v < 1 && u > v * std::numeric_limits<T>::max())
 387         return std::numeric_limits<T>::max();
 388     if (v > 1 && u < v * std::numeric_limits<T>::min())
 389         return 0;
 390     return u / v;
 391 }
 392
 393 // Floating point numbers comparison:
 394 // u is "essentially equal" [1][2] to v if: | u - v | / |u| <= e and | u - v | / |v| <= e
 395 //
 396 // [1] Knuth, D. E. "Accuracy of Floating Point Arithmetic." The Art of Computer Programming. 3rd ed. Vol. 2.
 397 //     Boston: Addison-Wesley, 1998. 229-45.
 398 // [2] http://www.boost.org/doc/libs/1_34_0/libs/test/doc/components/test_tools/floating_point_comparison.html
 399 template <typename T>
 400 inline typename std::enable_if<std::is_floating_point<T>::value, bool>::type areEssentiallyEqual(T u, T v, T epsilon = std::numeric_limits<T>::epsilon())
 401 {
 402     if (u == v)
 403         return true;
 404
 405     const T delta = std::abs(u - v);
 406     return safeFPDivision(delta, std::abs(u)) <= epsilon && safeFPDivision(delta, std::abs(v)) <= epsilon;
 407 }
 408
 409 // Match behavior of Math.min, where NaN is returned if either argument is NaN.
 410 template <typename T>
 411 inline typename std::enable_if<std::is_floating_point<T>::value, T>::type nanPropagatingMin(T a, T b)
 412 {
 413     return std::isnan(a) || std::isnan(b) ? std::numeric_limits<T>::quiet_NaN() : std::min(a, b);
 414 }
 415
 416 // Match behavior of Math.max, where NaN is returned if either argument is NaN.
 417 template <typename T>
 418 inline typename std::enable_if<std::is_floating_point<T>::value, T>::type nanPropagatingMax(T a, T b)
 419 {
 420     return std::isnan(a) || std::isnan(b) ? std::numeric_limits<T>::quiet_NaN() : std::max(a, b);
 421 }
 422
 423 inline bool isIntegral(float value)
 424 {
 425     return static_cast<int>(value) == value;
 426 }
 427
 428 template<typename T>
 429 inline void incrementWithSaturation(T& value)
 430 {
 431     if (value != std::numeric_limits<T>::max())
 432         value++;
 433 }
 434
 435 template<typename T>
 436 inline T leftShiftWithSaturation(T value, unsigned shiftAmount, T max = std::numeric_limits<T>::max())
 437 {
 438     T result = value << shiftAmount;
 439     // We will have saturated if shifting right doesn't recover the original value.
 440     if (result >> shiftAmount != value)
 441         return max;
 442     if (result > max)
 443         return max;
 444     return result;
 445 }
 446
 447 // Check if two ranges overlap assuming that neither range is empty.
 448 template<typename T>
 449 inline bool nonEmptyRangesOverlap(T leftMin, T leftMax, T rightMin, T rightMax)
 450 {
 451     ASSERT(leftMin < leftMax);
 452     ASSERT(rightMin < rightMax);
 453
 454     return leftMax > rightMin && rightMax > leftMin;
 455 }
 456
 457 // Pass ranges with the min being inclusive and the max being exclusive. For example, this should
 458 // return false:
 459 //
 460 //     rangesOverlap(0, 8, 8, 16)
 461 template<typename T>
 462 inline bool rangesOverlap(T leftMin, T leftMax, T rightMin, T rightMax)
 463 {
 464     ASSERT(leftMin <= leftMax);
 465     ASSERT(rightMin <= rightMax);
 466
 467     // Empty ranges interfere with nothing.
 468     if (leftMin == leftMax)
 469         return false;
 470     if (rightMin == rightMax)
 471         return false;
 472
 473     return nonEmptyRangesOverlap(leftMin, leftMax, rightMin, rightMax);
 474 }
 475
 476 // This mask is not necessarily the minimal mask, specifically if size is
 477 // a power of 2. It has the advantage that it's fast to compute, however.
 478 inline uint32_t computeIndexingMask(uint32_t size)
 479 {
 480     return static_cast<uint64_t>(static_cast<uint32_t>(-1)) >> std::clz(size);
 481 }
 482
 483 constexpr unsigned preciseIndexMaskShiftForSize(unsigned size)
 484 {
 485     return size * 8 - 1;
 486 }
 487
 488 template<typename T>
 489 constexpr unsigned preciseIndexMaskShift()
 490 {
 491     return preciseIndexMaskShiftForSize(sizeof(T));
 492 }
 493
 494 template<typename T>
 495 T opaque(T pointer)
 496 {
 497 #if !OS(WINDOWS)
 498     asm("" : "+r"(pointer));
 499 #endif
 500     return pointer;
 501 }
 502
 503 // This masks the given pointer with 0xffffffffffffffff (ptrwidth) if `index <
 504 // length`. Otherwise, it masks the pointer with 0. Similar to Linux kernel's array_ptr.
 505 template<typename T>
 506 inline T* preciseIndexMaskPtr(uintptr_t index, uintptr_t length, T* value)
 507 {
 508     uintptr_t result = bitwise_cast<uintptr_t>(value) & static_cast<uintptr_t>(
 509         static_cast<intptr_t>(index - opaque(length)) >>
 510         static_cast<intptr_t>(preciseIndexMaskShift<T*>()));
 511     return bitwise_cast<T*>(result);
 512 }
 513
 514 template<typename VectorType, typename RandomFunc>
 515 void shuffleVector(VectorType& vector, size_t size, const RandomFunc& randomFunc)
 516 {
 517     for (size_t i = 0; i + 1 < size; ++i)
 518         std::swap(vector[i], vector[i + randomFunc(size - i)]);
 519 }
 520
 521 template<typename VectorType, typename RandomFunc>
 522 void shuffleVector(VectorType& vector, const RandomFunc& randomFunc)
 523 {
 524     shuffleVector(vector, vector.size(), randomFunc);
 525 }
 526
 527 inline unsigned clz32(uint32_t number)
 528 {
 529 #if COMPILER(GCC_COMPATIBLE)
 530     if (number)
 531         return __builtin_clz(number);
 532     return 32;
 533 #elif COMPILER(MSVC)
 534     // Visual Studio 2008 or upper have __lzcnt, but we can't detect Intel AVX at compile time.
 535     // So we use bit-scan-reverse operation to calculate clz.
 536     unsigned long ret = 0;
 537     if (_BitScanReverse(&ret, number))
 538         return 31 - ret;
 539     return 32;
 540 #else
 541     unsigned zeroCount = 0;
 542     for (int i = 31; i >= 0; i--) {
 543         if (!(number >> i))
 544             zeroCount++;
 545         else
 546             break;
 547     }
 548     return zeroCount;
 549 #endif
 550 }
 551
 552 inline unsigned clz64(uint64_t number)
 553 {
 554 #if COMPILER(GCC_COMPATIBLE)
 555     if (number)
 556         return __builtin_clzll(number);
 557     return 64;
 558 #elif COMPILER(MSVC) && !CPU(X86)
 559     // Visual Studio 2008 or upper have __lzcnt, but we can't detect Intel AVX at compile time.
 560     // So we use bit-scan-reverse operation to calculate clz.
 561     // _BitScanReverse64 is defined in X86_64 and ARM in MSVC supported environments.
 562     unsigned long ret = 0;
 563     if (_BitScanReverse64(&ret, number))
 564         return 63 - ret;
 565     return 64;
 566 #else
 567     unsigned zeroCount = 0;
 568     for (int i = 63; i >= 0; i--) {
 569         if (!(number >> i))
 570             zeroCount++;
 571         else
 572             break;
 573     }
 574     return zeroCount;
 575 #endif
 576 }
 577
 578 } // namespace WTF
 579
 580 using WTF::opaque;
 581 using WTF::preciseIndexMaskPtr;
 582 using WTF::preciseIndexMaskShift;
 583 using WTF::preciseIndexMaskShiftForSize;
 584 using WTF::shuffleVector;
 585 using WTF::clz32;
 586 using WTF::clz64;
 587
 588 #endif // #ifndef WTF_MathExtras_h