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 <concepts>
0014 #include <stdexcept>
0015 #include <type_traits>
0016
0017 namespace Acts {
0018
0019 /// Returns the absolute of a number
0020 /// @note Can be removed for C++23
0021 /// @param n The number to take absolute value of
0022 /// @return The absolute value of the input
0023 template <typename T>
0024 constexpr T abs(const T n) {
0025   if (std::is_constant_evaluated()) {
0026     if constexpr (std::is_signed_v<T>) {
0027       if (n < 0) {
0028         return -n;
0029       }
0030     }
0031     return n;
0032   } else {
0033     return std::abs(n);
0034   }
0035 }
0036 /// Copies the sign of a signed variable onto the copyTo input object. Return
0037 /// type & magnitude remain unaffected by this method which allows usage for
0038 /// Vectors & other types providing the - operator. By convention, the zero is
0039 /// assigned to a positive sign.
0040 /// @param copyTo Variable to which the sign is copied to.
0041 /// @param sign Variable from which the sign is taken.
0042 /// @return The copyTo variable with the sign of the sign parameter
0043 template <typename out_t, typename sign_t>
0044 constexpr out_t copySign(const out_t& copyTo, const sign_t& sign) {
0045   if constexpr (std::is_enum_v<sign_t>) {
0046     return copySign(copyTo, static_cast<std::underlying_type_t<sign_t>>(sign));
0047   } else {
0048     constexpr sign_t zero = 0;
0049     return sign >= zero ? copyTo : -copyTo;
0050   }
0051 }
0052
0053 /// Calculates the ordinary power of the number x.
0054 /// @param x Number to take the power from
0055 /// @param p Power to take
0056 /// @return x raised to the power p
0057 template <typename T, std::integral P>
0058 constexpr T pow(T x, P p) {
0059   if (std::is_constant_evaluated()) {
0060     constexpr T one = 1;
0061     if constexpr (std::is_signed_v<P>) {
0062       if (p < 0 && abs(x) > std::numeric_limits<T>::epsilon()) {
0063         x = one / x;
0064         p = -p;
0065       }
0066     }
0067     using unsigned_p = std::make_unsigned_t<P>;
0068     return p == 0 ? one : x * pow(x, static_cast<unsigned_p>(p) - 1);
0069   } else {
0070     return static_cast<T>(std::pow(x, static_cast<T>(p)));
0071   }
0072 }
0073
0074 /// Returns the square of the passed number
0075 /// @param x The number to square
0076 /// @return The square of the input
0077 template <typename T>
0078 constexpr auto square(T x) {
0079   return x * x;
0080 }
0081
0082 /// Calculates the sum of squares of arguments
0083 /// @param args Variable number of arguments to square and sum
0084 /// @return Sum of squares of all arguments
0085 template <typename... T>
0086 constexpr auto hypotSquare(T... args) {
0087   return (square(args) + ...);
0088 }
0089
0090 /// Fast hypotenuse calculation for multiple arguments
0091 /// @param args Variable number of arguments
0092 /// @return Square root of sum of squares of arguments
0093 template <typename... T>
0094 constexpr auto fastHypot(T... args) {
0095   return std::sqrt(hypotSquare(args...));
0096 }
0097 /// Overload for single argument
0098 /// @param arg Single argument for which the hypotenuse is calculated
0099 /// @return Absolute value of the argument
0100 template <typename T>
0101 constexpr auto fastHypot(T arg) {
0102   return std::abs(arg);
0103 }
0104
0105 /// Slow but more accurate hypotenuse calculation for multiple arguments
0106 /// @param args Variable number of arguments
0107 /// @return Square root of sum of squares of arguments
0108 template <typename... T>
0109 constexpr auto slowHypot(T... args) {
0110   decltype((args + ...)) r = 0;
0111   ((r = std::hypot(r, args)), ...);
0112   return r;
0113 }
0114 /// Overload for single argument
0115 /// @param arg Single argument for which the hypotenuse is calculated
0116 /// @return Absolute value of the argument
0117 template <typename T>
0118 constexpr auto slowHypot(T arg) {
0119   return std::abs(arg);
0120 }
0121 /// Overload for two arguments
0122 /// @param x First argument
0123 /// @param y Second argument
0124 /// @return Square root of sum of squares of x and y
0125 template <typename T>
0126 constexpr auto slowHypot(T x, T y) {
0127   return std::hypot(x, y);
0128 }
0129 /// Overload for three arguments
0130 /// @param x First argument
0131 /// @param y Second argument
0132 /// @param z Third argument
0133 /// @return Square root of sum of squares of x, y and z
0134 template <typename T>
0135 constexpr auto slowHypot(T x, T y, T z) {
0136   return std::hypot(x, y, z);
0137 }
0138
0139 /// Calculates the squared cathetus of arguments, i.e. the difference between
0140 /// the square of the hypotenuse and the square of the given arguments without
0141 /// checking for domain errors.
0142 /// @param hypotenuse The hypotenuse value
0143 /// @param args Variable number of arguments to calculate the cathetus for
0144 /// @return Difference between the square of the hypotenuse and the square of the given arguments
0145 template <typename T, typename... Args>
0146 constexpr auto cathetusSquareUnchecked(T hypotenuse, Args... args) {
0147   return square(hypotenuse) - hypotSquare(args...);
0148 }
0149
0150 /// Calculates the squared cathetus of arguments, i.e. the difference between
0151 /// the square of the hypotenuse and the square of the given arguments.
0152 /// @param hypotenuse The hypotenuse value
0153 /// @param args Variable number of arguments to calculate the cathetus for
0154 /// @return Difference between the square of the hypotenuse and the square of the
0155 ///   given arguments. Returns NaN for floating point types if the hypotenuse is
0156 ///   smaller than the cathetus.
0157 /// @throws std::domain_error if hypotenuse is smaller than the cathetus for
0158 ///   integral types
0159 template <typename T, typename... Args>
0160 constexpr auto cathetusSquare(T hypotenuse, Args... args) {
0161   const auto hypotenuseSquare = square(hypotenuse);
0162   const auto argsSquareSum = hypotSquare(args...);
0163   if (std::floating_point<T> && hypotenuseSquare < argsSquareSum) {
0164     return std::numeric_limits<T>::quiet_NaN();
0165   } else if (std::integral<T> && hypotenuseSquare < argsSquareSum) {
0166     throw std::domain_error(
0167         "hypotenuse must be greater than or equal to the cathetus");
0168   }
0169   return hypotenuseSquare - argsSquareSum;
0170 }
0171
0172 /// Fast cathetus calculation for multiple arguments
0173 /// @param hypotenuse The hypotenuse value
0174 /// @param args Variable number of arguments to calculate the cathetus for
0175 /// @return Square root of difference between the square of the hypotenuse and the square of the given arguments
0176 template <typename T, typename... Args>
0177 constexpr auto fastCathetus(T hypotenuse, Args... args) {
0178   const auto hypotArgs = fastHypot(args...);
0179   return std::sqrt((hypotenuse - hypotArgs) * (hypotenuse + hypotArgs));
0180 }
0181
0182 /// Slow but more accurate cathetus calculation for multiple arguments
0183 /// @note For 2 arguments, the fast and slow cathetus are identical as they both use the same formula.
0184 /// @param hypotenuse The hypotenuse value
0185 /// @param args Variable number of arguments to calculate the cathetus for
0186 /// @return Square root of difference between the square of the hypotenuse and the square of the given arguments
0187 template <typename T, typename... Args>
0188 constexpr auto slowCathetus(T hypotenuse, Args... args) {
0189   const auto hypotArgs = slowHypot(args...);
0190   return std::sqrt((hypotenuse - hypotArgs) * (hypotenuse + hypotArgs));
0191 }
0192
0193 /// Calculates the sum of 1 + 2 + 3+ ... + N using the Gaussian sum formula
0194 /// @param N Number until which the sum runs
0195 /// @return Sum of integers from 1 to N
0196 template <std::integral T>
0197 constexpr T sumUpToN(const T N) {
0198   return N * (N + 1) / 2;
0199 }
0200
0201 /// Calculates the product of all integers within the given integer range
0202 /// (nLower)(nLower+1)(...)(upper-1)(upper). If lowerN is bigger than upperN,
0203 /// the function returns one.
0204 /// @param lowerN Lower range of the product calculation
0205 /// @param upperN Upper range of the product calculation
0206 /// @return Product result
0207 template <std::unsigned_integral T>
0208 constexpr T product(const T lowerN, const T upperN) {
0209   if (lowerN == T{0}) {
0210     return T{0};
0211   }
0212
0213   T value{1};
0214   for (T iter = lowerN; iter <= upperN; ++iter) {
0215     if (std::is_constant_evaluated() &&
0216         value > std::numeric_limits<T>::max() / iter) {
0217       throw std::overflow_error("product overflow");
0218     }
0219     assert(value <= std::numeric_limits<T>::max() / iter);
0220
0221     value *= iter;
0222   }
0223
0224   return value;
0225 }
0226
0227 /// Calculate the the factorial of an integer
0228 /// @param N Number of which the factorial is to be calculated
0229 /// @return The factorial of N
0230 template <std::unsigned_integral T>
0231 constexpr T factorial(const T N) {
0232   return product<T>(1, N);
0233 }
0234
0235 /// Calculate the binomial coefficient
0236 ///              n        n!
0237 ///                 =  --------
0238 ///              k     k!(n-k)!
0239 /// @param n Upper value in binomial coefficient
0240 /// @param k Lower value in binomial coefficient
0241 /// @return Binomial coefficient n choose k
0242 template <std::unsigned_integral T>
0243 constexpr T binomial(const T n, const T k) {
0244   if (std::is_constant_evaluated() && k > n) {
0245     throw std::overflow_error("k must be <= n");
0246   }
0247   assert(k <= n && "k must be <= n");
0248
0249   return product<T>(n - k + 1, n) / factorial<T>(k);
0250 }
0251
0252 }  // namespace Acts