Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

feat: safe divide #3253

Closed
wants to merge 4 commits into from
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
43 changes: 43 additions & 0 deletions Core/include/Acts/Utilities/AlgebraHelpers.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -249,4 +249,47 @@ constexpr T safeExp(T val) noexcept {
return std::exp(val);
}

/// Specialisation of the save divide epsilon to be used for safe division,
/// depending on the floating point type.
template <typename T>
struct SafeDivideEpsilon {};
template <>
struct SafeDivideEpsilon<float> {
static constexpr float value = 1e-7f;
};
template <>
struct SafeDivideEpsilon<double> {
static constexpr double value = 1e-15;
};
template <>
struct SafeDivideEpsilon<long double> {
static constexpr long double value = 1e-19L;
};

/// Perform safe division while avoiding division by zero or near-zero values.
///
/// @param numerator the numerator of the division.
/// @param denominator the denominator of the division.
///
/// @return numerator / denominator if denominator is not zero or near-zero.
/// Otherwise it will round the denominator towards +-epsilon to ensure a proper
/// result. Zero will be rounded up.
template <typename T>
constexpr T safeDivide(T numerator, T denominator) {
static_assert(std::is_floating_point<T>::value,
"safeDivide requires floating-point types");

constexpr T epsilon = SafeDivideEpsilon<T>::value;

if (std::abs(denominator) < epsilon) {
if (denominator >= 0) {
return numerator / epsilon;
} else {
return -numerator / epsilon;
}
}

return numerator / denominator;
}

} // namespace Acts
175 changes: 175 additions & 0 deletions Tests/UnitTests/Core/Utilities/AlgebraHelpersTests.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,8 @@ namespace Acts::Test {

BOOST_AUTO_TEST_SUITE(AlgebraHelpers)

BOOST_AUTO_TEST_SUITE(SafeInverse)

ACTS_LOCAL_LOGGER(Acts::getDefaultLogger("SafeInverse", logLevel))

BOOST_AUTO_TEST_CASE(SafeInverseSmallMatrix) {
Expand Down Expand Up @@ -119,4 +121,177 @@ BOOST_AUTO_TEST_CASE(SafeInverseFPELargeMatrix) {

BOOST_AUTO_TEST_SUITE_END()

BOOST_AUTO_TEST_SUITE(SafeDivide)

ACTS_LOCAL_LOGGER(Acts::getDefaultLogger("SafeDivide", logLevel))

BOOST_AUTO_TEST_CASE(SafeDivideFloat) {
using FloatType = float;

const FloatType largerThanEpsilon = 1e-6f;
const FloatType epsilon = 1e-7f;
const FloatType smallerThanEpsilon = 1e-8f;

const FloatType zero = 0.;
const FloatType a = 2.;

{
const FloatType denom = 5.;
BOOST_CHECK_EQUAL(Acts::safeDivide(a, denom), a / denom);
BOOST_CHECK_EQUAL(Acts::safeDivide(-a, denom), -a / denom);
BOOST_CHECK_EQUAL(Acts::safeDivide(a, -denom), -a / denom);
BOOST_CHECK_EQUAL(Acts::safeDivide(-a, -denom), a / denom);
BOOST_CHECK_EQUAL(Acts::safeDivide(zero, denom), zero);
BOOST_CHECK_EQUAL(Acts::safeDivide(zero, -denom), zero);
}
{
const FloatType denom = largerThanEpsilon;
BOOST_CHECK_EQUAL(Acts::safeDivide(a, denom), a / denom);
BOOST_CHECK_EQUAL(Acts::safeDivide(-a, denom), -a / denom);
BOOST_CHECK_EQUAL(Acts::safeDivide(a, -denom), -a / denom);
BOOST_CHECK_EQUAL(Acts::safeDivide(-a, -denom), a / denom);
BOOST_CHECK_EQUAL(Acts::safeDivide(zero, denom), zero);
BOOST_CHECK_EQUAL(Acts::safeDivide(zero, -denom), zero);
}
{
const FloatType denom = epsilon;
BOOST_CHECK_EQUAL(Acts::safeDivide(a, denom), a / denom);
BOOST_CHECK_EQUAL(Acts::safeDivide(-a, denom), -a / denom);
BOOST_CHECK_EQUAL(Acts::safeDivide(a, -denom), -a / denom);
BOOST_CHECK_EQUAL(Acts::safeDivide(-a, -denom), a / denom);
BOOST_CHECK_EQUAL(Acts::safeDivide(zero, denom), zero);
BOOST_CHECK_EQUAL(Acts::safeDivide(zero, -denom), zero);
}
{
const FloatType denom = smallerThanEpsilon;
BOOST_CHECK_EQUAL(Acts::safeDivide(a, denom), a / epsilon);
BOOST_CHECK_EQUAL(Acts::safeDivide(-a, denom), -a / epsilon);
BOOST_CHECK_EQUAL(Acts::safeDivide(a, -denom), -a / epsilon);
BOOST_CHECK_EQUAL(Acts::safeDivide(-a, -denom), a / epsilon);
BOOST_CHECK_EQUAL(Acts::safeDivide(zero, denom), zero);
BOOST_CHECK_EQUAL(Acts::safeDivide(zero, -denom), zero);
}
{
const FloatType denom = zero;
BOOST_CHECK_EQUAL(Acts::safeDivide(a, denom), a / epsilon);
BOOST_CHECK_EQUAL(Acts::safeDivide(-a, denom), -a / epsilon);

BOOST_CHECK_EQUAL(Acts::safeDivide(zero, denom), zero);
}
}

BOOST_AUTO_TEST_CASE(SafeDivideDouble) {
using FloatType = double;

const FloatType largerThanEpsilon = 1e-14;
const FloatType epsilon = 1e-15;
const FloatType smallerThanEpsilon = 1e-16;

const FloatType zero = 0.;
const FloatType a = 2.;

{
const FloatType denom = 5.;
BOOST_CHECK_EQUAL(Acts::safeDivide(a, denom), a / denom);
BOOST_CHECK_EQUAL(Acts::safeDivide(-a, denom), -a / denom);
BOOST_CHECK_EQUAL(Acts::safeDivide(a, -denom), -a / denom);
BOOST_CHECK_EQUAL(Acts::safeDivide(-a, -denom), a / denom);
BOOST_CHECK_EQUAL(Acts::safeDivide(zero, denom), zero);
BOOST_CHECK_EQUAL(Acts::safeDivide(zero, -denom), zero);
}
{
const FloatType denom = largerThanEpsilon;
BOOST_CHECK_EQUAL(Acts::safeDivide(a, denom), a / denom);
BOOST_CHECK_EQUAL(Acts::safeDivide(-a, denom), -a / denom);
BOOST_CHECK_EQUAL(Acts::safeDivide(a, -denom), -a / denom);
BOOST_CHECK_EQUAL(Acts::safeDivide(-a, -denom), a / denom);
BOOST_CHECK_EQUAL(Acts::safeDivide(zero, denom), zero);
BOOST_CHECK_EQUAL(Acts::safeDivide(zero, -denom), zero);
}
{
const FloatType denom = epsilon;
BOOST_CHECK_EQUAL(Acts::safeDivide(a, denom), a / denom);
BOOST_CHECK_EQUAL(Acts::safeDivide(-a, denom), -a / denom);
BOOST_CHECK_EQUAL(Acts::safeDivide(a, -denom), -a / denom);
BOOST_CHECK_EQUAL(Acts::safeDivide(-a, -denom), a / denom);
BOOST_CHECK_EQUAL(Acts::safeDivide(zero, denom), zero);
BOOST_CHECK_EQUAL(Acts::safeDivide(zero, -denom), zero);
}
{
const FloatType denom = smallerThanEpsilon;
BOOST_CHECK_EQUAL(Acts::safeDivide(a, denom), a / epsilon);
BOOST_CHECK_EQUAL(Acts::safeDivide(-a, denom), -a / epsilon);
BOOST_CHECK_EQUAL(Acts::safeDivide(a, -denom), -a / epsilon);
BOOST_CHECK_EQUAL(Acts::safeDivide(-a, -denom), a / epsilon);
BOOST_CHECK_EQUAL(Acts::safeDivide(zero, denom), zero);
BOOST_CHECK_EQUAL(Acts::safeDivide(zero, -denom), zero);
}
{
const FloatType denom = zero;
BOOST_CHECK_EQUAL(Acts::safeDivide(a, denom), a / epsilon);
BOOST_CHECK_EQUAL(Acts::safeDivide(-a, denom), -a / epsilon);

BOOST_CHECK_EQUAL(Acts::safeDivide(zero, denom), zero);
}
}

BOOST_AUTO_TEST_CASE(SafeDivideLongDouble) {
using FloatType = long double;

const FloatType largerThanEpsilon = 1e-18L;
const FloatType epsilon = 1e-19L;
const FloatType smallerThanEpsilon = 1e-20L;

const FloatType zero = 0.;
const FloatType a = 2.;

{
const FloatType denom = 5.;
BOOST_CHECK_EQUAL(Acts::safeDivide(a, denom), a / denom);
BOOST_CHECK_EQUAL(Acts::safeDivide(-a, denom), -a / denom);
BOOST_CHECK_EQUAL(Acts::safeDivide(a, -denom), -a / denom);
BOOST_CHECK_EQUAL(Acts::safeDivide(-a, -denom), a / denom);
BOOST_CHECK_EQUAL(Acts::safeDivide(zero, denom), zero);
BOOST_CHECK_EQUAL(Acts::safeDivide(zero, -denom), zero);
}
{
const FloatType denom = largerThanEpsilon;
BOOST_CHECK_EQUAL(Acts::safeDivide(a, denom), a / denom);
BOOST_CHECK_EQUAL(Acts::safeDivide(-a, denom), -a / denom);
BOOST_CHECK_EQUAL(Acts::safeDivide(a, -denom), -a / denom);
BOOST_CHECK_EQUAL(Acts::safeDivide(-a, -denom), a / denom);
BOOST_CHECK_EQUAL(Acts::safeDivide(zero, denom), zero);
BOOST_CHECK_EQUAL(Acts::safeDivide(zero, -denom), zero);
}
{
const FloatType denom = epsilon;
BOOST_CHECK_EQUAL(Acts::safeDivide(a, denom), a / denom);
BOOST_CHECK_EQUAL(Acts::safeDivide(-a, denom), -a / denom);
BOOST_CHECK_EQUAL(Acts::safeDivide(a, -denom), -a / denom);
BOOST_CHECK_EQUAL(Acts::safeDivide(-a, -denom), a / denom);
BOOST_CHECK_EQUAL(Acts::safeDivide(zero, denom), zero);
BOOST_CHECK_EQUAL(Acts::safeDivide(zero, -denom), zero);
}
{
const FloatType denom = smallerThanEpsilon;
BOOST_CHECK_EQUAL(Acts::safeDivide(a, denom), a / epsilon);
BOOST_CHECK_EQUAL(Acts::safeDivide(-a, denom), -a / epsilon);
BOOST_CHECK_EQUAL(Acts::safeDivide(a, -denom), -a / epsilon);
BOOST_CHECK_EQUAL(Acts::safeDivide(-a, -denom), a / epsilon);
BOOST_CHECK_EQUAL(Acts::safeDivide(zero, denom), zero);
BOOST_CHECK_EQUAL(Acts::safeDivide(zero, -denom), zero);
}
{
const FloatType denom = zero;
BOOST_CHECK_EQUAL(Acts::safeDivide(a, denom), a / epsilon);
BOOST_CHECK_EQUAL(Acts::safeDivide(-a, denom), -a / epsilon);

BOOST_CHECK_EQUAL(Acts::safeDivide(zero, denom), zero);
}
}

BOOST_AUTO_TEST_SUITE_END()

BOOST_AUTO_TEST_SUITE_END()

} // namespace Acts::Test
Loading