Author: luc
Date: Sun Apr 3 14:33:29 2011
New Revision: 1088316
URL: http://svn.apache.org/viewvc?rev=1088316&view=rev
Log:
Reduced cancellation errors in Vector3D.crossProduct
Jira: MATH-554
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math/geometry/Vector3D.java
commons/proper/math/trunk/src/site/xdoc/changes.xml
commons/proper/math/trunk/src/test/java/org/apache/commons/math/geometry/Vector3DTest.java
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math/geometry/Vector3D.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/geometry/Vector3D.java?rev=1088316&r1=1088315&r2=1088316&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math/geometry/Vector3D.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math/geometry/Vector3D.java
Sun Apr 3 14:33:29 2011
@@ -454,10 +454,41 @@ public class Vector3D implements Seriali
* @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.
Modified: commons/proper/math/trunk/src/site/xdoc/changes.xml
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/changes.xml?rev=1088316&r1=1088315&r2=1088316&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/changes.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/changes.xml Sun Apr 3 14:33:29 2011
@@ -52,6 +52,9 @@ The <action> type attribute can be add,u
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>
Modified:
commons/proper/math/trunk/src/test/java/org/apache/commons/math/geometry/Vector3DTest.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/geometry/Vector3DTest.java?rev=1088316&r1=1088315&r2=1088316&view=diff
==============================================================================
---
commons/proper/math/trunk/src/test/java/org/apache/commons/math/geometry/Vector3DTest.java
(original)
+++
commons/proper/math/trunk/src/test/java/org/apache/commons/math/geometry/Vector3DTest.java
Sun Apr 3 14:33:29 2011
@@ -153,6 +153,19 @@ public class Vector3DTest {
}
@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);
Assert.assertEquals(0, Vector3D.PLUS_I.getDelta(), 1.0e-10);
