Add mixed dense/skew-symmetric arithmetic operators

libeigen/eigen!2308

Closes #2925

Eigen/src/Core/SkewSymmetricMatrix3.h

diff --git a/Eigen/src/Core/SkewSymmetricMatrix3.h b/Eigen/src/Core/SkewSymmetricMatrix3.h
index a6ad143..72b65e1 100644
--- a/Eigen/src/Core/SkewSymmetricMatrix3.h
+++ b/Eigen/src/Core/SkewSymmetricMatrix3.h
@@ -159,6 +159,61 @@
       const SkewSymmetricBase<OtherDerived>& other) const {
     return (vector() - other.vector()).asSkewSymmetric();
   }
+
+  // Return type of dense +/- skew. Scalar follows ScalarBinaryOpTraits, so
+  // mixed-scalar cases go through the same promotion machinery as the rest of
+  // Eigen (and are rejected at compile time unless a user-provided
+  // specialization makes them valid). Shape is always 3x3.
+  template <typename OtherDerived, typename BinaryOp>
+  using DenseSkewBinaryReturnType = Matrix<
+      typename ScalarBinaryOpTraits<typename internal::traits<OtherDerived>::Scalar, Scalar, BinaryOp>::ReturnType, 3,
+      3>;
+
+  template <typename OtherDerived>
+  using DenseSkewSumReturnType =
+      DenseSkewBinaryReturnType<OtherDerived,
+                                internal::scalar_sum_op<typename internal::traits<OtherDerived>::Scalar, Scalar>>;
+
+  template <typename OtherDerived>
+  using DenseSkewDifferenceReturnType = DenseSkewBinaryReturnType<
+      OtherDerived, internal::scalar_difference_op<typename internal::traits<OtherDerived>::Scalar, Scalar>>;
+
+  /** \returns the sum of a dense matrix \a lhs and the skew symmetric matrix \a rhs as a dense matrix.
+   *
+   * The LHS must be 3x3 at compile time (or Dynamic, which is checked at runtime by the conversion).
+   * Only the skew side is materialized via \c toDenseMatrix(); the LHS remains a lazy expression
+   * until the enclosing assignment. */
+  template <typename OtherDerived>
+  EIGEN_DEVICE_FUNC friend EIGEN_STRONG_INLINE DenseSkewSumReturnType<OtherDerived> operator+(
+      const MatrixBase<OtherDerived>& lhs, const SkewSymmetricBase& rhs) {
+    EIGEN_STATIC_ASSERT_SAME_MATRIX_SIZE(OtherDerived, DenseMatrixType);
+    return lhs.derived() + rhs.toDenseMatrix();
+  }
+
+  /** \returns the sum of the skew symmetric matrix \a lhs and a dense matrix \a rhs as a dense matrix.
+   *
+   * Sum is commutative, so this forwards to the \c dense+skew overload. */
+  template <typename OtherDerived>
+  EIGEN_DEVICE_FUNC friend EIGEN_STRONG_INLINE DenseSkewSumReturnType<OtherDerived> operator+(
+      const SkewSymmetricBase& lhs, const MatrixBase<OtherDerived>& rhs) {
+    return rhs + lhs;
+  }
+
+  /** \returns the difference of a dense matrix \a lhs and the skew symmetric matrix \a rhs as a dense matrix. */
+  template <typename OtherDerived>
+  EIGEN_DEVICE_FUNC friend EIGEN_STRONG_INLINE DenseSkewDifferenceReturnType<OtherDerived> operator-(
+      const MatrixBase<OtherDerived>& lhs, const SkewSymmetricBase& rhs) {
+    EIGEN_STATIC_ASSERT_SAME_MATRIX_SIZE(OtherDerived, DenseMatrixType);
+    return lhs.derived() - rhs.toDenseMatrix();
+  }
+
+  /** \returns the difference of the skew symmetric matrix \a lhs and a dense matrix \a rhs as a dense matrix. */
+  template <typename OtherDerived>
+  EIGEN_DEVICE_FUNC friend EIGEN_STRONG_INLINE DenseSkewDifferenceReturnType<OtherDerived> operator-(
+      const SkewSymmetricBase& lhs, const MatrixBase<OtherDerived>& rhs) {
+    EIGEN_STATIC_ASSERT_SAME_MATRIX_SIZE(OtherDerived, DenseMatrixType);
+    return lhs.toDenseMatrix() - rhs.derived();
+  }
 };
 
 /** \class SkewSymmetricMatrix3

test/skew_symmetric_matrix3.cpp

diff --git a/test/skew_symmetric_matrix3.cpp b/test/skew_symmetric_matrix3.cpp
index 48d490d..a7fcb47 100644
--- a/test/skew_symmetric_matrix3.cpp
+++ b/test/skew_symmetric_matrix3.cpp
@@ -81,6 +81,44 @@
   VERIFY_IS_APPROX((sk1 + sk1).vector(), 2 * v1);
   VERIFY((sk1 - sk1).vector().isZero());
   VERIFY((sk1 - sk1).toDenseMatrix().isZero());
+
+  // Mixed dense / skew-symmetric +/- operators (see issue #2925, MR !2308).
+  // Cases c1/c2/c3 below are exactly the three expressions from the issue.
+  {
+    const SquareMatrix A = SquareMatrix::Random();
+    const Vector b = Vector::Random();
+
+    const SquareMatrix skA = (A * b).asSkewSymmetric().toDenseMatrix();
+    const SquareMatrix skB = v1.asSkewSymmetric().toDenseMatrix();
+    const SquareMatrix skB_times_A = v1.asSkewSymmetric().toDenseMatrix() * A;
+    const SquareMatrix A_times_skB = A * v1.asSkewSymmetric().toDenseMatrix();
+
+    // c1: dense - skew
+    const SquareMatrix c1 = A - SkewSymmetricMatrix3<Scalar>(A * b);
+    VERIFY_IS_APPROX(c1, A - skA);
+
+    // c2: Product<dense,skew> - skew
+    const SquareMatrix c2 = A * SkewSymmetricMatrix3<Scalar>(v1) - SkewSymmetricMatrix3<Scalar>(A * b);
+    VERIFY_IS_APPROX(c2, A_times_skB - skA);
+
+    // c3: Product<skew,dense> - skew
+    const SquareMatrix c3 = SkewSymmetricMatrix3<Scalar>(v1) * A - SkewSymmetricMatrix3<Scalar>(A * b);
+    VERIFY_IS_APPROX(c3, skB_times_A - skA);
+
+    // M + S and S + M (commutativity check).
+    VERIFY_IS_APPROX(SquareMatrix(A + sk1), A + skB);
+    VERIFY_IS_APPROX(SquareMatrix(sk1 + A), skB + A);
+    VERIFY_IS_APPROX(SquareMatrix(A + sk1), SquareMatrix(sk1 + A));
+
+    // M - S and S - M (anti-commutativity sanity).
+    VERIFY_IS_APPROX(SquareMatrix(A - sk1), A - skB);
+    VERIFY_IS_APPROX(SquareMatrix(sk1 - A), skB - A);
+    VERIFY_IS_APPROX(SquareMatrix(A - sk1), -SquareMatrix(sk1 - A));
+
+    // Product expression on the LHS should still fuse (no explicit
+    // .toDenseMatrix() needed on the dense side).
+    VERIFY_IS_APPROX(SquareMatrix(A * A - sk1), (A * A).eval() - skB);
+  }
 }
 
 template <typename Scalar>