fix cross product for complexes and add support for mixed real-complex cross products
diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index f318bfd..096438f 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h
@@ -399,8 +399,16 @@ /////////// Geometry module /////////// + #ifndef EIGEN_PARSED_BY_DOXYGEN + /// \internal helper struct to form the return type of the cross product + template<typename OtherDerived> struct cross_product_return_type { + typedef typename internal::scalar_product_traits<typename internal::traits<Derived>::Scalar,typename internal::traits<OtherDerived>::Scalar>::ReturnType Scalar; + typedef Matrix<Scalar,RowsAtCompileTime,ColsAtCompileTime> type; + }; + #endif // EIGEN_PARSED_BY_DOXYGEN template<typename OtherDerived> - PlainObject cross(const MatrixBase<OtherDerived>& other) const; + typename cross_product_return_type<OtherDerived>::type + cross(const MatrixBase<OtherDerived>& other) const; template<typename OtherDerived> PlainObject cross3(const MatrixBase<OtherDerived>& other) const; PlainObject unitOrthogonal(void) const;
diff --git a/Eigen/src/Geometry/OrthoMethods.h b/Eigen/src/Geometry/OrthoMethods.h index 64ca52e..52b4698 100644 --- a/Eigen/src/Geometry/OrthoMethods.h +++ b/Eigen/src/Geometry/OrthoMethods.h
@@ -35,7 +35,7 @@ */ template<typename Derived> template<typename OtherDerived> -inline typename MatrixBase<Derived>::PlainObject +inline typename MatrixBase<Derived>::template cross_product_return_type<OtherDerived>::type MatrixBase<Derived>::cross(const MatrixBase<OtherDerived>& other) const { EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(Derived,3) @@ -45,10 +45,10 @@ // optimize such a small temporary very well (even within a complex expression) const typename internal::nested<Derived,2>::type lhs(derived()); const typename internal::nested<OtherDerived,2>::type rhs(other.derived()); - return typename internal::plain_matrix_type<Derived>::type( - lhs.coeff(1) * rhs.coeff(2) - lhs.coeff(2) * rhs.coeff(1), - lhs.coeff(2) * rhs.coeff(0) - lhs.coeff(0) * rhs.coeff(2), - lhs.coeff(0) * rhs.coeff(1) - lhs.coeff(1) * rhs.coeff(0) + return typename cross_product_return_type<OtherDerived>::type( + internal::conj(lhs.coeff(1) * rhs.coeff(2) - lhs.coeff(2) * rhs.coeff(1)), + internal::conj(lhs.coeff(2) * rhs.coeff(0) - lhs.coeff(0) * rhs.coeff(2)), + internal::conj(lhs.coeff(0) * rhs.coeff(1) - lhs.coeff(1) * rhs.coeff(0)) ); } @@ -62,9 +62,9 @@ run(const VectorLhs& lhs, const VectorRhs& rhs) { return typename internal::plain_matrix_type<VectorLhs>::type( - lhs.coeff(1) * rhs.coeff(2) - lhs.coeff(2) * rhs.coeff(1), - lhs.coeff(2) * rhs.coeff(0) - lhs.coeff(0) * rhs.coeff(2), - lhs.coeff(0) * rhs.coeff(1) - lhs.coeff(1) * rhs.coeff(0), + internal::conj(lhs.coeff(1) * rhs.coeff(2) - lhs.coeff(2) * rhs.coeff(1)), + internal::conj(lhs.coeff(2) * rhs.coeff(0) - lhs.coeff(0) * rhs.coeff(2)), + internal::conj(lhs.coeff(0) * rhs.coeff(1) - lhs.coeff(1) * rhs.coeff(0)), 0 ); } @@ -121,16 +121,16 @@ if(Direction==Vertical) { eigen_assert(CrossReturnType::RowsAtCompileTime==3 && "the matrix must have exactly 3 rows"); - res.row(0) = _expression().row(1) * other.coeff(2) - _expression().row(2) * other.coeff(1); - res.row(1) = _expression().row(2) * other.coeff(0) - _expression().row(0) * other.coeff(2); - res.row(2) = _expression().row(0) * other.coeff(1) - _expression().row(1) * other.coeff(0); + res.row(0) = (_expression().row(1) * other.coeff(2) - _expression().row(2) * other.coeff(1)).conjugate(); + res.row(1) = (_expression().row(2) * other.coeff(0) - _expression().row(0) * other.coeff(2)).conjugate(); + res.row(2) = (_expression().row(0) * other.coeff(1) - _expression().row(1) * other.coeff(0)).conjugate(); } else { eigen_assert(CrossReturnType::ColsAtCompileTime==3 && "the matrix must have exactly 3 columns"); - res.col(0) = _expression().col(1) * other.coeff(2) - _expression().col(2) * other.coeff(1); - res.col(1) = _expression().col(2) * other.coeff(0) - _expression().col(0) * other.coeff(2); - res.col(2) = _expression().col(0) * other.coeff(1) - _expression().col(1) * other.coeff(0); + res.col(0) = (_expression().col(1) * other.coeff(2) - _expression().col(2) * other.coeff(1)).conjugate(); + res.col(1) = (_expression().col(2) * other.coeff(0) - _expression().col(0) * other.coeff(2)).conjugate(); + res.col(2) = (_expression().col(0) * other.coeff(1) - _expression().col(1) * other.coeff(0)).conjugate(); } return res; }
diff --git a/test/geo_orthomethods.cpp b/test/geo_orthomethods.cpp index 9152d99..020ae71 100644 --- a/test/geo_orthomethods.cpp +++ b/test/geo_orthomethods.cpp
@@ -33,6 +33,7 @@ template<typename Scalar> void orthomethods_3() { + typedef typename NumTraits<Scalar>::Real RealScalar; typedef Matrix<Scalar,3,3> Matrix3; typedef Matrix<Scalar,3,1> Vector3; @@ -44,6 +45,9 @@ // cross product VERIFY_IS_MUCH_SMALLER_THAN(v1.cross(v2).dot(v1), Scalar(1)); + VERIFY_IS_MUCH_SMALLER_THAN(v1.dot(v1.cross(v2)), Scalar(1)); + VERIFY_IS_MUCH_SMALLER_THAN(v1.cross(v2).dot(v2), Scalar(1)); + VERIFY_IS_MUCH_SMALLER_THAN(v2.dot(v1.cross(v2)), Scalar(1)); Matrix3 mat3; mat3 << v0.normalized(), (v0.cross(v1)).normalized(), @@ -68,6 +72,12 @@ v40.w() = v41.w() = v42.w() = 0; v42.template head<3>() = v40.template head<3>().cross(v41.template head<3>()); VERIFY_IS_APPROX(v40.cross3(v41), v42); + + // check mixed product + typedef Matrix<RealScalar, 3, 1> RealVector3; + RealVector3 rv1 = RealVector3::Random(); + VERIFY_IS_APPROX(v1.cross(rv1.template cast<Scalar>()), v1.cross(rv1)); + VERIFY_IS_APPROX(rv1.template cast<Scalar>().cross(v1), rv1.cross(v1)); } template<typename Scalar, int Size> void orthomethods(int size=Size) @@ -115,6 +125,7 @@ for(int i = 0; i < g_repeat; i++) { CALL_SUBTEST_1( orthomethods_3<float>() ); CALL_SUBTEST_2( orthomethods_3<double>() ); + CALL_SUBTEST_4( orthomethods_3<std::complex<double> >() ); CALL_SUBTEST_1( (orthomethods<float,2>()) ); CALL_SUBTEST_2( (orthomethods<double,2>()) ); CALL_SUBTEST_1( (orthomethods<float,3>()) );