Acts/Utilities/QuickMath.hpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2024 CERN for the benefit of the Acts project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
0008
0009 #pragma once
0010
0011 #include <cmath>
0012 #include <cstdint>
0013 #include <limits>
0014 #include <type_traits>
0015
0016 namespace Acts {
0017
0018 /// @brief Fast inverse square root function
0019 /// Taken from https://en.wikipedia.org/wiki/Fast_inverse_square_root
0020 /// @param x the number to calculate the inverse square root of
0021 /// @param iterations the number of newton iterations to perform
0022 template <typename T>
0023 constexpr T fastInverseSqrt(T x, int iterations = 1) noexcept {
0024   static_assert(std::numeric_limits<T>::is_iec559 &&
0025                 "Only IEEE 754 is supported");
0026   static_assert(std::is_same_v<T, float> || std::is_same_v<T, double>,
0027                 "Only float and double are supported");
0028   using I = std::conditional_t<std::is_same_v<T, float>, std::uint32_t,
0029                                std::uint64_t>;
0030
0031   constexpr I magic =
0032       std::is_same_v<T, float> ? 0x5f3759df : 0x5fe6eb50c7b537a9;
0033
0034   union {
0035     T f;
0036     I i;
0037   } u = {x};
0038
0039   u.i = magic - (u.i >> 1);
0040
0041   for (int i = 0; i < iterations; ++i) {
0042     u.f *= 1.5f - 0.5f * x * u.f * u.f;
0043   }
0044
0045   return u.f;
0046 }
0047
0048 /// @brief Fast power function @see `std::pow`
0049 /// Taken from
0050 /// https://martin.ankerl.com/2012/01/25/optimized-approximative-pow-in-c-and-cpp
0051 /// @param a the base
0052 /// @param b the exponent
0053 constexpr double fastPow(double a, double b) {
0054   // enable only on IEEE 754
0055   static_assert(std::numeric_limits<double>::is_iec559);
0056
0057   constexpr std::int64_t magic = 0x3FEF127F00000000;
0058
0059   union {
0060     double f;
0061     std::int64_t i;
0062   } u = {a};
0063
0064   u.i = static_cast<std::int64_t>(b * (u.i - magic) + magic);
0065
0066   return u.f;
0067 }
0068
0069 /// @brief Fast power function more precise than @see `fastPow`
0070 /// Taken from
0071 /// https://martin.ankerl.com/2012/01/25/optimized-approximative-pow-in-c-and-cpp
0072 /// @param a the base
0073 /// @param b the exponent
0074 constexpr double fastPowMorePrecise(double a, double b) {
0075   // binary exponentiation
0076   double r = 1.0;
0077   int exp = std::abs(static_cast<int>(b));
0078   double base = b > 0 ? a : 1 / a;
0079   while (exp != 0) {
0080     if ((exp & 1) != 0) {
0081       r *= base;
0082     }
0083     base *= base;
0084     exp >>= 1;
0085   }
0086
0087   return r * fastPow(a, b - static_cast<int>(b));
0088 }
0089
0090 }  // namespace Acts