Acts/Utilities/MathHelpers.hpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/.
0008
0009 #pragma once
0010
0011 #include <cassert>
0012 #include <cmath>
0013 #include <type_traits>
0014
0015 namespace Acts {
0016
0017 /// @brief Returns the absolute of a number
0018 ///        (Can be removed for c++ 23)
0019 /// @param n The number to take absolute value of
0020 /// @return The absolute value of the input
0021 template <typename T>
0022 constexpr T abs(const T n) {
0023   if (std::is_constant_evaluated()) {
0024     if constexpr (std::is_signed_v<T>) {
0025       if (n < 0) {
0026         return -n;
0027       }
0028     }
0029     return n;
0030   } else {
0031     return std::abs(n);
0032   }
0033 }
0034 /// @brief Copies the sign of a signed variable onto the copyTo input object
0035 ///        Return type & magnitude remain unaffected by this method which allows
0036 ///        usage for Vectors & other types providing the - operator.
0037 ///        By convention, the zero is assigned to a positive sign.
0038 /// @param copyTo: Variable to which the sign is copied to.
0039 /// @param sign: Variable from which the sign is taken.
0040 template <typename out_t, typename sign_t>
0041 constexpr out_t copySign(const out_t& copyTo, const sign_t& sign) {
0042   if constexpr (std::is_enum_v<sign_t>) {
0043     return copySign(copyTo, static_cast<std::underlying_type_t<sign_t>>(sign));
0044   } else {
0045     constexpr sign_t zero = 0;
0046     return sign >= zero ? copyTo : -copyTo;
0047   }
0048 }
0049
0050 /// @brief Calculates the ordinary power of the number x.
0051 /// @param x: Number to take the power from
0052 /// @param p: Power to take
0053 /// @return x raised to the power p
0054 template <typename T, std::integral P>
0055 constexpr T pow(T x, P p) {
0056   if (std::is_constant_evaluated()) {
0057     constexpr T one = 1;
0058     if constexpr (std::is_signed_v<P>) {
0059       if (p < 0 && abs(x) > std::numeric_limits<T>::epsilon()) {
0060         x = one / x;
0061         p = -p;
0062       }
0063     }
0064     using unsigned_p = std::make_unsigned_t<P>;
0065     return p == 0 ? one : x * pow(x, static_cast<unsigned_p>(p) - 1);
0066   } else {
0067     return static_cast<T>(std::pow(x, static_cast<T>(p)));
0068   }
0069 }
0070
0071 /// @brief Returns the square of the passed number
0072 /// @param x The number to square
0073 /// @return The square of the input
0074 template <typename T>
0075 constexpr auto square(T x) {
0076   return x * x;
0077 }
0078
0079 /// @brief Calculates the sum of squares of arguments
0080 /// @param args Variable number of arguments to square and sum
0081 /// @return Sum of squares of all arguments
0082 template <typename... T>
0083 constexpr auto hypotSquare(T... args) {
0084   return (square(args) + ...);
0085 }
0086
0087 /// @brief Fast hypotenuse calculation for multiple arguments
0088 /// @param args Variable number of arguments
0089 /// @return Square root of sum of squares of arguments
0090 template <typename... T>
0091 constexpr auto fastHypot(T... args) {
0092   return std::sqrt(hypotSquare(args...));
0093 }
0094
0095 /// @brief Calculates the sum of 1 + 2 + 3+ ... + N using the
0096 ///        Gaussian sum formula
0097 /// @param N: Number until which the sum runs
0098 /// @return Sum of integers from 1 to N
0099 template <std::integral T>
0100 constexpr T sumUpToN(const T N) {
0101   return N * (N + 1) / 2;
0102 }
0103 /// @brief Calculates the product of all integers
0104 ///        within the given integer range
0105 ///           (nLower)(nLower +1)(...)(upper-1)(upper)
0106 ///        If lowerN is bigger than upperN, the function
0107 ///        returns one
0108 /// @param lowerN: Lower range of the product calculation
0109 /// @param upperN: Upper range of the product calculation
0110 /// @return Factorial result
0111 template <std::unsigned_integral T>
0112 constexpr T product(const T lowerN, const T upperN) {
0113   if (lowerN == static_cast<T>(0)) {
0114     return 0;
0115   }
0116   T value{1};
0117   for (T iter = std::max(static_cast<T>(2), lowerN); iter <= upperN; ++iter) {
0118     assert(value * iter > value);
0119     value *= iter;
0120   }
0121   return value;
0122 }
0123 /// @brief Calculate the the factorial of an integer
0124 /// @param N: Number of which the factorial is to be calculated
0125 template <std::unsigned_integral T>
0126 constexpr T factorial(const T N) {
0127   return product<T>(1, N);
0128 }
0129 /// @brief Calculate the binomial coefficient
0130 ///              n        n!
0131 ///                 =  --------
0132 ///              k     k!(n-k)!
0133 /// @param n Upper value in binomial coefficient
0134 /// @param k Lower value in binomial coefficient
0135 /// @return Binomial coefficient n choose k
0136 template <std::unsigned_integral T>
0137 constexpr T binomial(const T n, const T k) {
0138   assert(k <= n);
0139   return product<T>(n - k + 1, n) / factorial<T>(k);
0140 }
0141
0142 }  // namespace Acts