Reduced cancellation errors in Vector3D.crossProduct

Jira: MATH-554

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1088316 13f79535-47bb-0310-9956-ffa450edef68

src/main/java/org/apache/commons/math/geometry/Vector3D.java

@@ -454,10 +454,41 @@ public class Vector3D implements Serializable {
    * @param v2 second vector
    * @return the cross product v1 ^ v2 as a new Vector
    */
-  public static Vector3D crossProduct(Vector3D v1, Vector3D v2) {
-    return new Vector3D(v1.y * v2.z - v1.z * v2.y,
-                        v1.z * v2.x - v1.x * v2.z,
-                        v1.x * v2.y - v1.y * v2.x);
+  public static Vector3D crossProduct(final Vector3D v1, final Vector3D v2) {
+
+      final double n1 = v1.getNormSq();
+      final double n2 = v2.getNormSq();
+      if ((n1 * n2) < MathUtils.SAFE_MIN) {
+          return ZERO;
+      }
+
+      // rescale both vectors without losing precision,
+      // to ensure their norm are the same order of magnitude
+      final int deltaExp = (FastMath.getExponent(n1) - FastMath.getExponent(n2)) / 4;
+      final double x1    = FastMath.scalb(v1.x, -deltaExp);
+      final double y1    = FastMath.scalb(v1.y, -deltaExp);
+      final double z1    = FastMath.scalb(v1.z, -deltaExp);
+      final double x2    = FastMath.scalb(v2.x,  deltaExp);
+      final double y2    = FastMath.scalb(v2.y,  deltaExp);
+      final double z2    = FastMath.scalb(v2.z,  deltaExp);
+
+      // we reduce cancellation errors by preconditioning,
+      // we replace v1 by v3 = v1 - rho v2 with rho chosen in order to compute
+      // v3 without loss of precision. See Kahan lecture
+      // "Computing Cross-Products and Rotations in 2- and 3-Dimensional Euclidean Spaces"
+      // available at http://www.cs.berkeley.edu/~wkahan/MathH110/Cross.pdf
+
+      // compute rho as an 8 bits approximation of v1.v2 / v2.v2
+      final double ratio = (x1 * x2 + y1 * y2 + z1 * z2) / FastMath.scalb(n2, 2 * deltaExp);
+      final double rho   = FastMath.rint(256 * ratio) / 256;
+
+      final double x3 = x1 - rho * x2;
+      final double y3 = y1 - rho * y2;
+      final double z3 = z1 - rho * z2;
+
+      // compute cross product from v3 and v2 instead of v1 and v2
+      return new Vector3D(y3 * z2 - z3 * y2, z3 * x2 - x3 * z2, x3 * y2 - y3 * x2);
+
   }
 
   /** Compute the distance between two vectors according to the L<sub>1</sub> norm.

src/site/xdoc/changes.xml

@@ -52,6 +52,9 @@ The <action> type attribute can be add,update,fix,remove.
     If the output is not quite correct, check for invisible trailing spaces!
      -->
     <release version="3.0" date="TBD" description="TBD">
+      <action dev="luc" type="fix" issue="MATH-554" >
+        Reduced cancellation errors in Vector3D.crossProduct
+      </action>
       <action dev="erans" type="fix" issue="MATH-552" due-to="James Bence">
         Fixed bug in "MultidimensionalCounter".
       </action>

src/test/java/org/apache/commons/math/geometry/Vector3DTest.java

@@ -152,6 +152,19 @@ public class Vector3DTest {
         Assert.assertTrue(FastMath.abs(Vector3D.dotProduct(v2, v3)) < 1.0e-12);
     }
 
+    @Test
+    public void testCrossProductCancellation() {
+        Vector3D v1 = new Vector3D(9070467121.0, 4535233560.0, 1);
+        Vector3D v2 = new Vector3D(9070467123.0, 4535233561.0, 1);
+        checkVector(Vector3D.crossProduct(v1, v2), -1, 2, 1);
+
+        double scale    = FastMath.scalb(1.0, 100);
+        Vector3D big1   = new Vector3D(scale, v1);
+        Vector3D small2 = new Vector3D(1 / scale, v2);
+        checkVector(Vector3D.crossProduct(big1, small2), -1, 2, 1);
+
+    }
+
     @Test
     public void testAngular() {
         Assert.assertEquals(0,           Vector3D.PLUS_I.getAlpha(), 1.0e-10);