[wip] Math: improve precision of cross().

While this makes 32-bit float cross product precision basically
equivalent to a 64-bit calculation casted back to 32-bit, it stays with
the speed halfway between the straightforward 32- and 64-bit
implementation.

However only on platforms that actually have a FMA instruction. For
example on Emscripten the code is TEN TIMES slower than the baseline
implementation, which is not an acceptable tradeoff -- there simply
using doubles to calculate the result is faster. And enabling the more
precise variant only on some platforms doesn't seem like a good idea for
portability.

Stashing this aside until I'm clearer what to do with this.

src/Magnum/Math/Test/Vector2Test.cpp

@@ -66,6 +66,7 @@ struct Vector2Test: Corrade::TestSuite::Tester {
 
     void access();
     void cross();
+    void crossCatastrophicCancellation();
     void axes();
     void scales();
     void perpendicular();
@@ -79,6 +80,7 @@ struct Vector2Test: Corrade::TestSuite::Tester {
 
 typedef Math::Vector3<Int> Vector3i;
 typedef Math::Vector2<Float> Vector2;
+typedef Math::Vector2<Double> Vector2d;
 typedef Math::Vector2<Int> Vector2i;
 
 Vector2Test::Vector2Test() {
@@ -92,6 +94,7 @@ Vector2Test::Vector2Test() {
 
               &Vector2Test::access,
               &Vector2Test::cross,
+              &Vector2Test::crossCatastrophicCancellation,
               &Vector2Test::axes,
               &Vector2Test::scales,
               &Vector2Test::perpendicular,
@@ -208,6 +211,18 @@ void Vector2Test::cross() {
     CORRADE_COMPARE(Math::cross<Int>({a, 0}, {b, 0}), Vector3i(0, 0, Math::cross(a, b)));
 }
 
+void Vector2Test::crossCatastrophicCancellation() {
+    /* Repro case and algorithm from
+       https://pharr.org/matt/blog/2019/11/03/difference-of-floats.html */
+
+    Vector2 a{33962.035f, 41563.4f};
+    Vector2 b{-24871.969f, -30438.8f};
+
+    Float expected = -75.1656f;
+    CORRADE_COMPARE(Float(Math::cross(Vector2d{a}, Vector2d{b})), expected);
+    CORRADE_COMPARE(Math::cross(a, b), expected);
+}
+
 void Vector2Test::axes() {
     constexpr Vector2 x = Vector2::xAxis(5.0f);
     constexpr Vector2 y = Vector2::yAxis(6.0f);

src/Magnum/Math/Test/Vector3Test.cpp

@@ -67,6 +67,7 @@ struct Vector3Test: Corrade::TestSuite::Tester {
 
     void access();
     void cross();
+    void crossCatastrophicCancellation();
     void axes();
     void scales();
     void twoComponent();
@@ -78,6 +79,7 @@ struct Vector3Test: Corrade::TestSuite::Tester {
 };
 
 typedef Math::Vector3<Float> Vector3;
+typedef Math::Vector3<Double> Vector3d;
 typedef Math::Vector3<Int> Vector3i;
 typedef Math::Vector2<Float> Vector2;
 
@@ -93,6 +95,7 @@ Vector3Test::Vector3Test() {
 
               &Vector3Test::access,
               &Vector3Test::cross,
+              &Vector3Test::crossCatastrophicCancellation,
               &Vector3Test::axes,
               &Vector3Test::scales,
               &Vector3Test::twoComponent,
@@ -228,6 +231,17 @@ void Vector3Test::cross() {
     CORRADE_COMPARE(Math::cross(a, b), Vector3i(-10, -3, 7));
 }
 
+void Vector3Test::crossCatastrophicCancellation() {
+    /* Repro case and algorithm from
+       https://pharr.org/matt/blog/2019/11/03/difference-of-floats.html */
+    Vector3 a{33962.035f, 41563.4f, 7706.415f};
+    Vector3 b{-24871.969f, -30438.8f, -5643.727f};
+
+    Vector3 expected{1556.0276f, -1257.5151f, -75.1656f};
+    CORRADE_COMPARE(Vector3{Math::cross(Vector3d{a}, Vector3d{b})}, expected);
+    CORRADE_COMPARE(Math::cross(a, b), expected);
+}
+
 void Vector3Test::axes() {
     constexpr Vector3 x = Vector3::xAxis(5.0f);
     constexpr Vector3 y = Vector3::yAxis(6.0f);

src/Magnum/Math/Vector2.h

@@ -33,6 +33,16 @@
 
 namespace Magnum { namespace Math {
 
+namespace Implementation {
+    /* https://pharr.org/matt/blog/2019/11/03/difference-of-floats.html */
+    template<class T> inline T differenceOfProducts(T a, T b, T c, T d) {
+        T cd = c*d;
+        T err = std::fma(-c, d, cd);
+        T dop = std::fma(a, b, -cd);
+        return dop + err;
+    }
+}
+
 /**
 @brief 2D cross product
 
@@ -44,11 +54,15 @@ are parallel or antiparallel and `1` when two *normalized* vectors are
 perpendicular. @f[
     \boldsymbol a \times \boldsymbol b = \boldsymbol a_\bot \cdot \boldsymbol b = a_xb_y - a_yb_x
 @f]
+
+The precision of the calculation is further improved using Kahan's algorithm
+as described in [an article by Matt Pharr](https://pharr.org/matt/blog/2019/11/03/difference-of-floats.html).
 @see @ref Vector2::perpendicular(),
     @ref dot(const Vector<size, T>&, const Vector<size, T>&)
  */
 template<class T> inline T cross(const Vector2<T>& a, const Vector2<T>& b) {
-    return a._data[0]*b._data[1] - a._data[1]*b._data[0];
+    return Implementation::differenceOfProducts(a._data[0], b._data[1],
+                                                a._data[1], b._data[0]);
 }
 
 /**

src/Magnum/Math/Vector3.h

@@ -41,13 +41,19 @@ parallel or antiparallel and length of @cpp 1 @ce when two *normalized* vectors
 are perpendicular. @f[
     \boldsymbol a \times \boldsymbol b = \begin{pmatrix}a_yb_z - a_zb_y \\ a_zb_x - a_xb_z \\ a_xb_y - a_yb_x \end{pmatrix}
 @f]
+
+The precision of the calculation is further improved using Kahan's algorithm
+as described in [an article by Matt Pharr](https://pharr.org/matt/blog/2019/11/03/difference-of-floats.html).
 @see @ref cross(const Vector2<T>&, const Vector2<T>&), @ref planeEquation()
 */
 template<class T> inline Vector3<T> cross(const Vector3<T>& a, const Vector3<T>& b) {
     return {
-        a._data[1]*b._data[2] - b._data[1]*a._data[2],
+        Implementation::differenceOfProducts(a._data[1], b._data[2],
-        a._data[2]*b._data[0] - b._data[2]*a._data[0],
+                                             a._data[2], b._data[1]),
-        a._data[0]*b._data[1] - b._data[0]*a._data[1]
+        Implementation::differenceOfProducts(a._data[2], b._data[0],
+                                             a._data[0], b._data[2]),
+        Implementation::differenceOfProducts(a._data[0], b._data[1],
+                                             a._data[1], b._data[0])
     };
 }