fix: potential floating point over/underflow (#1709)

fixes potential floating point over/underflow cases. discovered by regex search `(.+)\s+\*\s+\1\s+[<>=]`

Core/include/Acts/Geometry/GeometryObjectSorter.hpp

@@ -43,31 +43,31 @@ class ObjectSorterT {
     switch (m_binningValue) {
       // compare on x
       case binX: {
-        return (one.x() < two.x());
+        return one.x() < two.x();
       }
       // compare on y
       case binY: {
-        return (one.y() < two.y());
+        return one.y() < two.y();
       }
       // compare on z
       case binZ: {
-        return (one.z() < two.z());
+        return one.z() < two.z();
       }
       // compare on r
       case binR: {
-        return (perp(one) < perp(two));
+        return perp(one) < perp(two);
       }
       // compare on phi
       case binPhi: {
-        return (phi(one) < phi(two));
+        return phi(one) < phi(two);
       }
       // compare on eta
       case binEta: {
-        return (eta(one) < eta(two));
+        return eta(one) < eta(two);
       }
       // default for the moment
       default: {
-        return (one.norm() < two.norm());
+        return one.norm() < two.norm();
       }
     }
   }
@@ -110,43 +110,43 @@ class DistanceSorterT {
       case binX: {
         double diffOneX = one.x() - m_reference.x();
         double diffTwoX = two.x() - m_reference.x();
-        return (diffOneX * diffOneX < diffTwoX * diffTwoX);
+        return std::abs(diffOneX) < std::abs(diffTwoX);
       }
       // compare on diff y
       case binY: {
         double diffOneY = one.y() - m_reference.y();
         double diffTwoY = two.y() - m_reference.y();
-        return (diffOneY * diffOneY < diffTwoY * diffTwoY);
+        return std::abs(diffOneY) < std::abs(diffTwoY);
       }
       // compare on diff z
       case binZ: {
         double diffOneZ = one.z() - m_reference.z();
         double diffTwoZ = two.z() - m_reference.z();
-        return (diffOneZ * diffOneZ < diffTwoZ * diffTwoZ);
+        return std::abs(diffOneZ) < std::abs(diffTwoZ);
       }
       // compare on r
       case binR: {
         double diffOneR = perp(one) - m_refR;
         double diffTwoR = perp(two) - m_refR;
-        return (diffOneR * diffOneR < diffTwoR * diffTwoR);
+        return std::abs(diffOneR) < std::abs(diffTwoR);
       }
       // compare on phi /// @todo add cyclic value
       case binPhi: {
         double diffOnePhi = phi(one) - m_refPhi;
         double diffTwoPhi = phi(two) - m_refPhi;
-        return (diffOnePhi * diffOnePhi < diffTwoPhi * diffTwoPhi);
+        return std::abs(diffOnePhi) < std::abs(diffTwoPhi);
       }
       // compare on eta
       case binEta: {
         double diffOneEta = eta(one) - m_refEta;
         double diffTwoEta = eta(two) - m_refEta;
-        return (diffOneEta * diffOneEta < diffTwoEta * diffTwoEta);
+        return std::abs(diffOneEta) < std::abs(diffTwoEta);
       }
       // default for the moment
       default: {
         T diffOne(one - m_reference);
         T diffTwo(two - m_reference);
-        return (diffOne.mag2() < diffTwo.mag2());
+        return diffOne.mag2() < diffTwo.mag2();
       }
     }
   }

Core/include/Acts/Surfaces/detail/PlanarHelper.hpp

@@ -37,7 +37,7 @@ inline Intersection3D intersect(const Transform3& transform,
     ActsScalar path = (pnormal.dot((pcenter - position))) / (denom);
     // Is valid hence either on surface or reachable
     Intersection3D::Status status =
-        (path * path < s_onSurfaceTolerance * s_onSurfaceTolerance)
+        std::abs(path) < std::abs(s_onSurfaceTolerance)
             ? Intersection3D::Status::onSurface
             : Intersection3D::Status::reachable;
     // Return the intersection

Core/src/Geometry/CylinderVolumeHelper.cpp

@@ -496,10 +496,13 @@ bool Acts::CylinderVolumeHelper::estimateAndCheckDimension(
       "Estimate/check CylinderVolumeBounds from/w.r.t. enclosed "
       "layers + envelope covers");
   // the z from the layers w and w/o envelopes
-  double zEstFromLayerEnv = 0.5 * ((layerZmax) + (layerZmin));
-  double halflengthFromLayer = 0.5 * std::abs((layerZmax) - (layerZmin));
+  double zEstFromLayerEnv = 0.5 * (layerZmax + layerZmin);
+  double halflengthFromLayer = 0.5 * std::abs(layerZmax - layerZmin);
 
-  bool concentric = (zEstFromLayerEnv * zEstFromLayerEnv < 0.001);
+  // using `sqrt(0.001) = 0.031622777` because previously it compared to
+  // `zEstFromLayerEnv * zEstFromLayerEnv < 0.001` which was changed due to
+  // underflow/overflow concerns
+  bool concentric = std::abs(zEstFromLayerEnv) < 0.031622777;
 
   bool idTrf = transform.isApprox(Transform3::Identity());
 

Core/src/Geometry/TrackingVolume.cpp

@@ -608,8 +608,7 @@ Acts::TrackingVolume::compatibleLayers(
         auto atIntersection =
             tLayer->surfaceOnApproach(gctx, position, direction, options);
         auto path = atIntersection.intersection.pathLength;
-        bool withinLimit =
-            (path * path <= options.pathLimit * options.pathLimit);
+        bool withinLimit = std::abs(path) <= std::abs(options.pathLimit);
         // Intersection is ok - take it (move to surface on appraoch)
         if (atIntersection &&
             (atIntersection.object != options.targetSurface) && withinLimit) {

Core/src/Surfaces/ConeSurface.cpp

@@ -292,7 +292,7 @@ Acts::SurfaceIntersection Acts::ConeSurface::intersect(
   // Check the validity of the first solution
   Vector3 solution1 = position + qe.first * direction;
   Intersection3D::Status status1 =
-      (qe.first * qe.first < s_onSurfaceTolerance * s_onSurfaceTolerance)
+      std::abs(qe.first) < std::abs(s_onSurfaceTolerance)
           ? Intersection3D::Status::onSurface
           : Intersection3D::Status::reachable;
 
@@ -303,7 +303,7 @@ Acts::SurfaceIntersection Acts::ConeSurface::intersect(
   // Check the validity of the second solution
   Vector3 solution2 = position + qe.first * direction;
   Intersection3D::Status status2 =
-      (qe.second * qe.second < s_onSurfaceTolerance * s_onSurfaceTolerance)
+      std::abs(qe.second) < std::abs(s_onSurfaceTolerance)
           ? Intersection3D::Status::onSurface
           : Intersection3D::Status::reachable;
   if (bcheck and not isOnSurface(gctx, solution2, direction, bcheck)) {
@@ -320,7 +320,7 @@ Acts::SurfaceIntersection Acts::ConeSurface::intersect(
                 (status1 == Intersection3D::Status::missed and
                  status2 == Intersection3D::Status::missed);
   // Check and (eventually) go with the first solution
-  if ((check1 and qe.first * qe.first < qe.second * qe.second) or
+  if ((check1 and (std::abs(qe.first) < std::abs(qe.second))) or
       status2 == Intersection3D::Status::missed) {
     // And add the alternative
     if (qe.solutions > 1) {

Core/src/Surfaces/CylinderSurface.cpp

@@ -224,7 +224,7 @@ Acts::SurfaceIntersection Acts::CylinderSurface::intersect(
   // Check the validity of the first solution
   Vector3 solution1 = position + qe.first * direction;
   Intersection3D::Status status1 =
-      qe.first * qe.first < s_onSurfaceTolerance * s_onSurfaceTolerance
+      std::abs(qe.first) < std::abs(s_onSurfaceTolerance)
           ? Intersection3D::Status::onSurface
           : Intersection3D::Status::reachable;
 
@@ -247,7 +247,8 @@ Acts::SurfaceIntersection Acts::CylinderSurface::intersect(
       double cZ = vecLocal.dot(tMatrix.block<3, 1>(0, 2));
       double tolerance = s_onSurfaceTolerance + bcheck.tolerance()[eBoundLoc1];
       double hZ = cBounds.get(CylinderBounds::eHalfLengthZ) + tolerance;
-      return (cZ * cZ < hZ * hZ) ? status : Intersection3D::Status::missed;
+      return std::abs(cZ) < std::abs(hZ) ? status
+                                         : Intersection3D::Status::missed;
     }
     return (isOnSurface(gctx, solution, direction, bcheck)
                 ? status
@@ -264,7 +265,7 @@ Acts::SurfaceIntersection Acts::CylinderSurface::intersect(
   // Check the validity of the second solution
   Vector3 solution2 = position + qe.second * direction;
   Intersection3D::Status status2 =
-      qe.second * qe.second < s_onSurfaceTolerance * s_onSurfaceTolerance
+      std::abs(qe.second) < std::abs(s_onSurfaceTolerance)
           ? Intersection3D::Status::onSurface
           : Intersection3D::Status::reachable;
   // Check first solution for boundary compatiblity
@@ -275,7 +276,7 @@ Acts::SurfaceIntersection Acts::CylinderSurface::intersect(
                 (status1 == Intersection3D::Status::missed and
                  status2 == Intersection3D::Status::missed);
   // Check and (eventually) go with the first solution
-  if ((check1 and qe.first * qe.first < qe.second * qe.second) or
+  if ((check1 and (std::abs(qe.first) < std::abs(qe.second))) or
       status2 == Intersection3D::Status::missed) {
     // And add the alternative
     cIntersection.alternative = second;

Core/src/Surfaces/DiscSurface.cpp

@@ -85,7 +85,7 @@ Acts::Result<Acts::Vector2> Acts::DiscSurface::globalToLocal(
     const Vector3& /*gmom*/, double tolerance) const {
   // transport it to the globalframe
   Vector3 loc3Dframe = (transform(gctx).inverse()) * position;
-  if (loc3Dframe.z() * loc3Dframe.z() > tolerance * tolerance) {
+  if (std::abs(loc3Dframe.z()) > std::abs(tolerance)) {
     return Result<Vector2>::failure(SurfaceError::GlobalPositionNotOnSurface);
   }
   return Result<Acts::Vector2>::success({perp(loc3Dframe), phi(loc3Dframe)});

Core/src/Surfaces/IntersectionHelper2D.cpp

@@ -55,7 +55,7 @@ Acts::detail::IntersectionHelper2D::intersectEllipse(ActsScalar Rx,
     ActsScalar solD = std::copysign(toSolD.norm(), toSolD.dot(dir));
     ActsScalar altD = std::copysign(toAltD.norm(), toAltD.dot(dir));
 
-    if (solD * solD < altD * altD) {
+    if (std::abs(solD) < std::abs(altD)) {
       return {Intersection2D(sol, solD, Intersection2D::Status::reachable),
               Intersection2D(alt, altD, Intersection2D::Status::reachable)};
     }

Core/src/Surfaces/LineSurface.cpp

@@ -141,10 +141,10 @@ Acts::SurfaceIntersection Acts::LineSurface::intersect(
   double denom = 1 - eaTeb * eaTeb;
   // validity parameter
   Intersection3D::Status status = Intersection3D::Status::unreachable;
-  if (denom * denom > s_onSurfaceTolerance * s_onSurfaceTolerance) {
+  if (std::abs(denom) > std::abs(s_onSurfaceTolerance)) {
     double u = (mab.dot(ea) - mab.dot(eb) * eaTeb) / denom;
     // Check if we are on the surface already
-    status = (u * u < s_onSurfaceTolerance * s_onSurfaceTolerance)
+    status = std::abs(u) < std::abs(s_onSurfaceTolerance)
                  ? Intersection3D::Status::onSurface
                  : Intersection3D::Status::reachable;
     Vector3 result = (ma + u * ea);
@@ -156,7 +156,7 @@ Acts::SurfaceIntersection Acts::LineSurface::intersect(
       double cZ = vecLocal.dot(eb);
       double hZ =
           m_bounds->get(LineBounds::eHalfLengthZ) + s_onSurfaceTolerance;
-      if ((cZ * cZ > hZ * hZ) or
+      if ((std::abs(cZ) > std::abs(hZ)) or
           ((vecLocal - cZ * eb).norm() >
            m_bounds->get(LineBounds::eR) + s_onSurfaceTolerance)) {
         status = Intersection3D::Status::missed;

Core/src/Surfaces/PlaneSurface.cpp

@@ -87,7 +87,7 @@ Acts::Result<Acts::Vector2> Acts::PlaneSurface::globalToLocal(
     const GeometryContext& gctx, const Vector3& position,
     const Vector3& /*unused*/, double tolerance) const {
   Vector3 loc3Dframe = transform(gctx).inverse() * position;
-  if (loc3Dframe.z() * loc3Dframe.z() > tolerance * tolerance) {
+  if (std::abs(loc3Dframe.z()) > std::abs(tolerance)) {
     return Result<Vector2>::failure(SurfaceError::GlobalPositionNotOnSurface);
   }
   return Result<Vector2>::success({loc3Dframe.x(), loc3Dframe.y()});

Tests/UnitTests/Core/Surfaces/CylinderSurfaceTests.cpp

@@ -180,7 +180,7 @@ BOOST_AUTO_TEST_CASE(CylinderSurfaceProperties) {
   // And it's path should be further away then the primary solution
   double pn = sfIntersection.intersection.pathLength;
   double pa = sfIntersection.alternative.pathLength;
-  BOOST_CHECK(pn * pn < pa * pa);
+  BOOST_CHECK(std::abs(pn) < std::abs(pa));
   BOOST_CHECK_EQUAL(sfIntersection.object, cylinderSurfaceObject.get());
 
   //