Fix #522: improve matrix identity test

Provide floating-point error margins in testing whether or not a matrix
is the identity matrix. This is done on a component by component basis.
Previously, the matrix values were compared directly to constants '0'
and '1', so the test could fail on account of floating-point errors.

(squash with parent commit before merging)

HEN_HOUSE/egs++/egs_transformations.h

@@ -218,11 +218,68 @@ class EGS_EXPORT EGS_RotationMatrix {
  /*! \brief Returns \c true, if this object is approximately the unit matrix, \a false
  otherwise */
  bool isI() const {
- // compute square distance to unit matrix
- EGS_Float d2 = (rxx-1)*(rxx-1) + rxy*rxy + rxz*rxz +
- ryx*ryx + (ryy-1)*(ryy-1) + ryz*ryz +
- rzx*rzx + rzy*rzy + (rzz-1)*(rzz-1);
- return (d2 < epsilon*epsilon);
+
+ EGS_Float d2;
+ EGS_Float eps2 = epsilon*epsilon;
+
+ // compute square distance to unit matrix, component-by-component;
+ // return as quickly as possible.
+
+ // rxy
+ d2 = rxy*rxy;
+ if (d2 > eps2) {
+ return false;
+ }
+
+ // rxz
+ d2 = rxz*rxz;
+ if (d2 > eps2) {
+ return false;
+ }
+
+ // ryx
+ d2 = ryx*ryx;
+ if (d2 > eps2) {
+ return false;
+ }
+
+ // ryz
+ d2 = ryz*ryz;
+ if (d2 > eps2) {
+ return false;
+ }
+
+ // rzx
+ d2 = rzx*rzx;
+ if (d2 > eps2) {
+ return false;
+ }
+
+ // rzy
+ d2 = rzy*rzy;
+ if (d2 > eps2) {
+ return false;
+ }
+
+ // rxx
+ EGS_Float d2 = (rxx-1)*(rxx-1);
+ if (d2 > eps2) {
+ return false
+ }
+
+ // ryy
+ EGS_Float d2 = (ryy-1)*(ryy-1);
+ if (d2 > eps2) {
+ return false
+ }
+
+ // rzz
+ EGS_Float d2 = (rzz-1)*(rzz-1);
+ if (d2 > eps2) {
+ return false
+ }
+
+ return true;
  };
 
  /*! \brief Returns the rotated a vector \f$ R \cdot \vec{v}\f$ */