Math.h

Go to the documentation of this file.

 
#pragma once
 
#include <Eigen/Core>
#include <cmath>
#include <vector>
 
namespace autopas::utils::Math {
const double normalScale = 1. / std::sqrt(2 * M_PI);
 
template <class Int_t, std::enable_if_t<std::is_integral_v<Int_t>, bool> = true>
Int_t safeAdd(const Int_t &a, const Int_t &b, const Int_t &valUnderflow = std::numeric_limits<Int_t>::min(),
              const Int_t &valOverflow = std::numeric_limits<Int_t>::max()) {
  Int_t result;
  bool overflow = __builtin_add_overflow(a, b, &result);
  if (overflow) {
    // if both args are negative this is an underflow.
    if (a < 0 and b < 0) {
      result = valUnderflow;
    } else {
      result = valOverflow;
    }
  }
  return result;
}
 
template <class Float_t, std::enable_if_t<std::is_floating_point_v<Float_t>, bool> = true>
Float_t safeAdd(const Float_t &a, const Float_t &b, const Float_t &valUnderflow = -std::numeric_limits<Float_t>::max(),
                const Float_t &valOverflow = std::numeric_limits<Float_t>::max()) {
  Float_t result = a + b;
  if (std::isinf(result)) {
    if (result > 0) {
      result = valOverflow;
    } else {
      result = valUnderflow;
    }
  }
  return result;
}
 
template <class Int_t, std::enable_if_t<std::is_integral_v<Int_t>, bool> = true>
Int_t safeSub(const Int_t &a, const Int_t &b, const Int_t &valUnderflow = std::numeric_limits<Int_t>::min(),
              const Int_t &valOverflow = std::numeric_limits<Int_t>::max()) {
  Int_t result;
  bool overflow = __builtin_sub_overflow(a, b, &result);
  if (overflow) {
    // underflow is only possible if we subtract a positive number. Everything else is an overflow.
    if (b > 0) {
      result = valUnderflow;
    } else {
      result = valOverflow;
    }
  }
  return result;
}
 
template <class Float_t, std::enable_if_t<std::is_floating_point_v<Float_t>, bool> = true>
Float_t safeSub(const Float_t &a, const Float_t &b, const Float_t &valUnderflow = -std::numeric_limits<Float_t>::max(),
                const Float_t &valOverflow = std::numeric_limits<Float_t>::max()) {
  Float_t result = a - b;
  if (std::isinf(result)) {
    if (result > 0) {
      result = valOverflow;
    } else {
      result = valUnderflow;
    }
  }
  return result;
}
 
template <class Int_t, std::enable_if_t<std::is_integral_v<Int_t>, bool> = true>
Int_t safeMul(const Int_t &a, const Int_t &b, const Int_t &valUnderflow = std::numeric_limits<Int_t>::min(),
              const Int_t &valOverflow = std::numeric_limits<Int_t>::max()) {
  Int_t result;
  bool overflow = __builtin_mul_overflow(a, b, &result);
  if (overflow) {
    // if exactly one arg is negative this is an underflow.
    if ((a < 0) xor (b < 0)) {
      result = valUnderflow;
    } else {
      result = valOverflow;
    }
  }
  return result;
}
 
template <class Float_t, std::enable_if_t<std::is_floating_point_v<Float_t>, bool> = true>
Float_t safeMul(const Float_t &a, const Float_t &b, const Float_t &valUnderflow = -std::numeric_limits<Float_t>::max(),
                const Float_t &valOverflow = std::numeric_limits<Float_t>::max()) {
  Float_t result = a * b;
  if (std::isinf(result)) {
    if (result > 0) {
      result = valOverflow;
    } else {
      result = valUnderflow;
    }
  }
  return result;
}
 
template <size_t exponent, class T>
T pow(const T &base) {
  if (exponent == 0) {
    return 1;
  }
 
  T res = base;
  // the compiler should unroll this loop
  for (size_t i = 0; i < exponent - 1; ++i) {
    res *= base;
  }
  return res;
}
 
double normalPDF(double x);
 
double normalCDF(double x);
 
double sigmoid(double x);
 
constexpr double EPSILON_RELATIVE_EQUALITY = 1e-9;
 
constexpr unsigned int MAX_ULP_DISTANCE = 4;
 
namespace internal {
template <typename T>
constexpr size_t bit_size_v = sizeof(T) * 8;
 
template <std::floating_point FloatType>
struct int_t_impl {
  static_assert(bit_size_v<FloatType> == 32 || bit_size_v<FloatType> == 64,
                "int_t only supports 32-bit and 64-bit floating-point types yet. Please extend to 128-bit types");
  using type = std::conditional_t<bit_size_v<FloatType> == 32, std::int32_t, std::int64_t>;
};
}  // namespace internal
 
template <std::floating_point FloatType>
using int_t = typename internal::int_t_impl<FloatType>::type;
 
template <std::floating_point FloatType>
bool isInUlp(FloatType lhs, FloatType rhs, unsigned int ulpDistance = MAX_ULP_DISTANCE) {
  // In case the floats are equal in their bitwise representation, return true
  if (lhs == rhs) {
    return true;
  }
 
  // In case the signs mismatch, return false
  if (lhs < static_cast<FloatType>(0.0) && rhs > static_cast<FloatType>(0.0) ||
      lhs > static_cast<FloatType>(0.0) && rhs < static_cast<FloatType>(0.0)) {
    return false;
  }
 
  // Compute ULP distance by interpreting the floating point as an equivalent-sized integer
  return reinterpret_cast<int_t<FloatType> &>(rhs) - reinterpret_cast<int_t<FloatType> &>(lhs) <= ulpDistance;
}
 
template <std::floating_point FloatType>
bool isNearRel(FloatType a, FloatType b, double maxRelativeDifference = EPSILON_RELATIVE_EQUALITY) {
  const auto greaterNumber = std::max(std::abs(a), std::abs(b));
  const auto absoluteDifference = maxRelativeDifference * greaterNumber;
  const auto diff = std::abs(a - b);
  return diff <= absoluteDifference;
}
 
template <std::floating_point FloatType>
bool isNearAbs(FloatType a, FloatType b, double maxAbsoluteDifference) {
  return std::abs(a - b) <= maxAbsoluteDifference;
}
 
double roundFixed(double d, int fixedPrecision);
 
double roundFloating(double d, int floatingPrecision);
 
Eigen::VectorXd makeVectorXd(const std::vector<double> &elements);
 
Eigen::VectorXi makeVectorXi(const std::vector<int> &elements);
}  // namespace autopas::utils::Math