side_by_triangle foils scaled-epsilon for equality by comparing determinant against zero

[Hubert Tong <hstong@…>]

In boost/geometry/strategies/cartesian/side_by_triangle.hpp:

/*    88*/              promoted_type const s 
/*    89*/                  = geometry::detail::determinant<promoted_type>
/*    90*/                      (
/*    91*/                          dx, dy, 
/*    92*/                          dpx, dpy
/*    93*/                      );
/*    94*/      
/*    95*/              promoted_type const zero = promoted_type();
/*    96*/              return math::equals(s, zero) ? 0 
/*    97*/                  : s > zero ? 1 
/*    98*/                  : -1;

The math::equals() interface implements, in boost/geometry/util/math.hpp, an equality comparison with an epsilon value—scaled by its arguments—when the common type is a floating-point type:

/*    52*/          static inline bool apply(Type const& a, Type const& b)
/*    53*/          {
/*    54*/                      if (a == b)
/*    55*/                      {
/*    56*/                              return true;
/*    57*/                      }
/*    58*/
/*    59*/              // See http://www.parashift.com/c++-faq-lite/newbie.html#faq-29.17,
/*    60*/              // FUTURE: replace by some boost tool or boost::test::close_at_tolerance
/*    61*/              return std::abs(a - b) <= std::numeric_limits<Type>::epsilon() * get_max(std::abs(a), std::abs(b), 1.0);

The scaling performed by the math::equals() interface fails to reliably reflect the scale of the original inputs to the determinant calculation (dx, dy, dpx, dpy) through no fault of its own; it was not passed the appropriate information.

In particular,

fma(x, y, -(x * y)); /* x*y - x*y */

is known to produce, on some platforms, the double value which is nearest to the round-off from (x * y). This round-off can legitimately be greater than DBL_EPSILON.

The determinant calculation in boost/geometry/arithmetic/determinant.hpp:

/*    42*/              return rt(ux) * rt(vy) - rt(uy) * rt(vx);

is liable to hit similar behaviour.

[awulkiew]

Thanks for the report! It took a while to address it, sorry about that.

For floating-point numbers, the result is now compared according to the machine epsilon scaled by the greatest value passed into the determinant() formula.

The results of the determinant are compared this way not only in side_by_triangle but also cartesian intersection strategy.

For the details see:

• ​https://github.com/boostorg/geometry/commit/b5887efd5e839e0e458bce207a34272cc3193a5c
• ​https://github.com/boostorg/geometry/commit/f59f6ffaceb6cb36f0b07ae24f01d8790527985c
• ​https://github.com/boostorg/geometry/commit/dbeb823fcbd8ef5c30f8c2a7f25f3ac2300a0f27

Any suggestions are welcome!