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 namespace Acts {
0015
0016 /// @brief Returns the absolute of a number
0017 ///        (Can be removed for c++ 23)
0018 /// @param n The number to take absolute value of
0019 /// @return The absolute value of the input
0020 template <typename T>
0021 constexpr T abs(const T n) {
0022   if (std::is_constant_evaluated()) {
0023     if constexpr (std::is_signed_v<T>) {
0024       if (n < 0) {
0025         return -n;
0026       }
0027     }
0028     return n;
0029   } else {
0030     return std::abs(n);
0031   }
0032 }
0033 /// @brief Copies the sign of a signed variable onto the copyTo input object
0034 ///        Return type & magnitude remain unaffected by this method which allows
0035 ///        usage for Vectors & other types providing the - operator.
0036 ///        By convention, the zero is assigned to a positive sign.
0037 /// @param copyTo: Variable to which the sign is copied to.
0038 /// @param sign: Variable from which the sign is taken.
0039 template <typename out_t, typename sign_t>
0040 constexpr out_t copySign(const out_t& copyTo, const sign_t& sign) {
0041   if constexpr (std::is_enum_v<sign_t>) {
0042     return copySign(copyTo, static_cast<std::underlying_type_t<sign_t>>(sign));
0043   } else {
0044     constexpr sign_t zero = 0;
0045     return sign >= zero ? copyTo : -copyTo;
0046   }
0047 }
0048
0049 /// @brief Calculates the ordinary power of the number x.
0050 /// @param x: Number to take the power from
0051 /// @param p: Power to take
0052 /// @return x raised to the power p
0053 template <typename T, std::integral P>
0054 constexpr T pow(T x, P p) {
0055   if (std::is_constant_evaluated()) {
0056     constexpr T one = 1;
0057     if constexpr (std::is_signed_v<P>) {
0058       if (p < 0 && abs(x) > std::numeric_limits<T>::epsilon()) {
0059         x = one / x;
0060         p = -p;
0061       }
0062     }
0063     using unsigned_p = std::make_unsigned_t<P>;
0064     return p == 0 ? one : x * pow(x, static_cast<unsigned_p>(p) - 1);
0065   } else {
0066     return static_cast<T>(std::pow(x, static_cast<T>(p)));
0067   }
0068 }
0069
0070 /// @brief Returns the square of the passed number
0071 /// @param x The number to square
0072 /// @return The square of the input
0073 template <typename T>
0074 constexpr auto square(T x) {
0075   return x * x;
0076 }
0077
0078 /// @brief Calculates the sum of squares of arguments
0079 /// @param args Variable number of arguments to square and sum
0080 /// @return Sum of squares of all arguments
0081 template <typename... T>
0082 constexpr auto hypotSquare(T... args) {
0083   return (square(args) + ...);
0084 }
0085
0086 /// @brief Fast hypotenuse calculation for multiple arguments
0087 /// @param args Variable number of arguments
0088 /// @return Square root of sum of squares of arguments
0089 template <typename... T>
0090 constexpr auto fastHypot(T... args) {
0091   return std::sqrt(hypotSquare(args...));
0092 }
0093
0094 /// @brief Calculates the sum of 1 + 2 + 3+ ... + N using the
0095 ///        Gaussian sum formula
0096 /// @param N: Number until which the sum runs
0097 /// @return Sum of integers from 1 to N
0098 template <std::integral T>
0099 constexpr T sumUpToN(const T N) {
0100   return N * (N + 1) / 2;
0101 }
0102 /// @brief Calculates the factorial of a number
0103 ///        N!= N*(N-1)....*3*2*1
0104 /// @param upperN: Upper factor until which the factorial is calculated
0105 /// @param lowerN: Optional argument to remove the first factors from the calculation
0106 /// @return Factorial result
0107 template <std::integral T>
0108 constexpr T factorial(const T upperN, const T lowerN = 1) {
0109   constexpr T one = 1;
0110   const T& limit = std::max(one, lowerN);
0111   return upperN >= limit ? upperN * factorial(upperN - 1, limit) : one;
0112 }
0113 /// @brief Calculate the binomial coefficient
0114 ///              n        n!
0115 ///                 =  --------
0116 ///              k     k!(n-k)!
0117 /// @param n Upper value in binomial coefficient
0118 /// @param k Lower value in binomial coefficient
0119 /// @return Binomial coefficient n choose k
0120 template <std::integral T>
0121 constexpr T binomial(const T n, const T k) {
0122   return factorial<T>(n, n - k + 1) / factorial<T>(k);
0123 }
0124
0125 }  // namespace Acts