EIC code displayed by LXR master/​acts/​Core/​include/​Acts/​Utilities/​MathHelpers.hpp	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

File indexing completed on 2026-05-23 07:34:31

 // This file is part of the ACTS project.
 //
 // Copyright (C) 2016 CERN for the benefit of the ACTS project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
 
 #pragma once
 
 #include <cassert>
 #include <cmath>
 #include <concepts>
 #include <stdexcept>
 #include <type_traits>
 
 namespace Acts {
 
 /// Returns the absolute of a number
 /// @note Can be removed for C++23
 /// @param n The number to take absolute value of
 /// @return The absolute value of the input
 template <typename T>
 constexpr T abs(const T n) {
   if (std::is_constant_evaluated()) {
     if constexpr (std::is_signed_v<T>) {
       if (n < 0) {
         return -n;
       }
     }
     return n;
   } else {
     return std::abs(n);
   }
 }
 /// Copies the sign of a signed variable onto the copyTo input object. Return
 /// type & magnitude remain unaffected by this method which allows usage for
 /// Vectors & other types providing the - operator. By convention, the zero is
 /// assigned to a positive sign.
 /// @param copyTo Variable to which the sign is copied to.
 /// @param sign Variable from which the sign is taken.
 /// @return The copyTo variable with the sign of the sign parameter
 template <typename out_t, typename sign_t>
 constexpr out_t copySign(const out_t& copyTo, const sign_t& sign) {
   if constexpr (std::is_enum_v<sign_t>) {
     return copySign(copyTo, static_cast<std::underlying_type_t<sign_t>>(sign));
   } else {
     constexpr sign_t zero = 0;
     return sign >= zero ? copyTo : -copyTo;
   }
 }
 
 /// Calculates the ordinary power of the number x.
 /// @param x Number to take the power from
 /// @param p Power to take
 /// @return x raised to the power p
 template <typename T, std::integral P>
 constexpr T pow(T x, P p) {
   if (std::is_constant_evaluated()) {
     constexpr T one = 1;
     if constexpr (std::is_signed_v<P>) {
       if (p < 0 && abs(x) > std::numeric_limits<T>::epsilon()) {
         x = one / x;
         p = -p;
       }
     }
     using unsigned_p = std::make_unsigned_t<P>;
     return p == 0 ? one : x * pow(x, static_cast<unsigned_p>(p) - 1);
   } else {
     return static_cast<T>(std::pow(x, static_cast<T>(p)));
   }
 }
 
 /// Returns the square of the passed number
 /// @param x The number to square
 /// @return The square of the input
 template <typename T>
 constexpr auto square(T x) {
   return x * x;
 }
 
 /// Calculates the sum of squares of arguments
 /// @param args Variable number of arguments to square and sum
 /// @return Sum of squares of all arguments
 template <typename... T>
 constexpr auto hypotSquare(T... args) {
   return (square(args) + ...);
 }
 
 /// Fast hypotenuse calculation for multiple arguments
 /// @param args Variable number of arguments
 /// @return Square root of sum of squares of arguments
 template <typename... T>
 constexpr auto fastHypot(T... args) {
   return std::sqrt(hypotSquare(args...));
 }
 /// Overload for single argument
 /// @param arg Single argument for which the hypotenuse is calculated
 /// @return Absolute value of the argument
 template <typename T>
 constexpr auto fastHypot(T arg) {
   return std::abs(arg);
 }
 
 /// Slow but more accurate hypotenuse calculation for multiple arguments
 /// @param args Variable number of arguments
 /// @return Square root of sum of squares of arguments
 template <typename... T>
 constexpr auto slowHypot(T... args) {
   decltype((args + ...)) r = 0;
   ((r = std::hypot(r, args)), ...);
   return r;
 }
 /// Overload for single argument
 /// @param arg Single argument for which the hypotenuse is calculated
 /// @return Absolute value of the argument
 template <typename T>
 constexpr auto slowHypot(T arg) {
   return std::abs(arg);
 }
 /// Overload for two arguments
 /// @param x First argument
 /// @param y Second argument
 /// @return Square root of sum of squares of x and y
 template <typename T>
 constexpr auto slowHypot(T x, T y) {
   return std::hypot(x, y);
 }
 /// Overload for three arguments
 /// @param x First argument
 /// @param y Second argument
 /// @param z Third argument
 /// @return Square root of sum of squares of x, y and z
 template <typename T>
 constexpr auto slowHypot(T x, T y, T z) {
   return std::hypot(x, y, z);
 }
 
 /// Calculates the squared cathetus of arguments, i.e. the difference between
 /// the square of the hypotenuse and the square of the given arguments without
 /// checking for domain errors.
 /// @param hypotenuse The hypotenuse value
 /// @param args Variable number of arguments to calculate the cathetus for
 /// @return Difference between the square of the hypotenuse and the square of the given arguments
 template <typename T, typename... Args>
 constexpr auto cathetusSquareUnchecked(T hypotenuse, Args... args) {
   return square(hypotenuse) - hypotSquare(args...);
 }
 
 /// Calculates the squared cathetus of arguments, i.e. the difference between
 /// the square of the hypotenuse and the square of the given arguments.
 /// @param hypotenuse The hypotenuse value
 /// @param args Variable number of arguments to calculate the cathetus for
 /// @return Difference between the square of the hypotenuse and the square of the
 ///   given arguments. Returns NaN for floating point types if the hypotenuse is
 ///   smaller than the cathetus.
 /// @throws std::domain_error if hypotenuse is smaller than the cathetus for
 ///   integral types
 template <typename T, typename... Args>
 constexpr auto cathetusSquare(T hypotenuse, Args... args) {
   const auto hypotenuseSquare = square(hypotenuse);
   const auto argsSquareSum = hypotSquare(args...);
   if (std::floating_point<T> && hypotenuseSquare < argsSquareSum) {
     return std::numeric_limits<T>::quiet_NaN();
   } else if (std::integral<T> && hypotenuseSquare < argsSquareSum) {
     throw std::domain_error(
         "hypotenuse must be greater than or equal to the cathetus");
   }
   return hypotenuseSquare - argsSquareSum;
 }
 
 /// Fast cathetus calculation for multiple arguments
 /// @param hypotenuse The hypotenuse value
 /// @param args Variable number of arguments to calculate the cathetus for
 /// @return Square root of difference between the square of the hypotenuse and the square of the given arguments
 template <typename T, typename... Args>
 constexpr auto fastCathetus(T hypotenuse, Args... args) {
   const auto hypotArgs = fastHypot(args...);
   return std::sqrt((hypotenuse - hypotArgs) * (hypotenuse + hypotArgs));
 }
 
 /// Slow but more accurate cathetus calculation for multiple arguments
 /// @note For 2 arguments, the fast and slow cathetus are identical as they both use the same formula.
 /// @param hypotenuse The hypotenuse value
 /// @param args Variable number of arguments to calculate the cathetus for
 /// @return Square root of difference between the square of the hypotenuse and the square of the given arguments
 template <typename T, typename... Args>
 constexpr auto slowCathetus(T hypotenuse, Args... args) {
   const auto hypotArgs = slowHypot(args...);
   return std::sqrt((hypotenuse - hypotArgs) * (hypotenuse + hypotArgs));
 }
 
 /// Calculates the sum of 1 + 2 + 3+ ... + N using the Gaussian sum formula
 /// @param N Number until which the sum runs
 /// @return Sum of integers from 1 to N
 template <std::integral T>
 constexpr T sumUpToN(const T N) {
   return N * (N + 1) / 2;
 }
 
 /// Calculates the product of all integers within the given integer range
 /// (nLower)(nLower+1)(...)(upper-1)(upper). If lowerN is bigger than upperN,
 /// the function returns one.
 /// @param lowerN Lower range of the product calculation
 /// @param upperN Upper range of the product calculation
 /// @return Product result
 template <std::unsigned_integral T>
 constexpr T product(const T lowerN, const T upperN) {
   if (lowerN == T{0}) {
     return T{0};
   }
 
   T value{1};
   for (T iter = lowerN; iter <= upperN; ++iter) {
     if (std::is_constant_evaluated() &&
         value > std::numeric_limits<T>::max() / iter) {
       throw std::overflow_error("product overflow");
     }
     assert(value <= std::numeric_limits<T>::max() / iter);
 
     value *= iter;
   }
 
   return value;
 }
 
 /// Calculate the the factorial of an integer
 /// @param N Number of which the factorial is to be calculated
 /// @return The factorial of N
 template <std::unsigned_integral T>
 constexpr T factorial(const T N) {
   return product<T>(1, N);
 }
 
 /// Calculate the binomial coefficient
 ///              n        n!
 ///                 =  --------
 ///              k     k!(n-k)!
 /// @param n Upper value in binomial coefficient
 /// @param k Lower value in binomial coefficient
 /// @return Binomial coefficient n choose k
 template <std::unsigned_integral T>
 constexpr T binomial(const T n, const T k) {
   if (std::is_constant_evaluated() && k > n) {
     throw std::overflow_error("k must be <= n");
   }
   assert(k <= n && "k must be <= n");
 
   return product<T>(n - k + 1, n) / factorial<T>(k);
 }
 
 /// Calculate the sinc function, defined as sin(x)/x, with a special handling
 /// for x=0 to avoid numerical instability.
 /// @param x The input value for which to calculate the sinc function
 /// @return The value of sinc(x)
 inline double sinc(double x) {
   // Numerical limit for double to get a different number than 1 from the first
   // order Taylor expansion of sin(x)/x ~ 1-x*x/6 around x=0.
   static const double eps =
       std::sqrt(std::numeric_limits<double>::epsilon()) * 6;
   if (std::abs(x) < eps) {
     return 1.0;
   }
   // Otherwise std::sin(x) ~ x is stable for small x
   return std::sin(x) / x;
 }
 
 }  // namespace Acts

This page was automatically generated by the 2.3.7 LXR engine. The LXR team