wip interior product

Change-Id: I672378174643d5bbf580375b61bc9552ed0d62cf

geometry/BUILD.bazel

@@ -58,6 +58,7 @@ cc_library(
         "src/exterior_product.hpp",
         "src/geometric_product.hpp",
         "src/get.hpp",
+        "src/interior_product.hpp",
         "src/multivector_for.hpp",
         "src/regressive_product.hpp",
         "src/sum.hpp",

geometry/geometry.hpp

@@ -8,6 +8,7 @@
 #include "geometry/src/exterior_product.hpp"
 #include "geometry/src/geometric_product.hpp"
 #include "geometry/src/get.hpp"
+#include "geometry/src/interior_product.hpp"
 #include "geometry/src/multivector_for.hpp"
 #include "geometry/src/regressive_product.hpp"
 #include "geometry/src/sum.hpp"

geometry/src/exterior_product.hpp

@@ -20,12 +20,12 @@ inline constexpr class
       class D1 = std::decay_t<typename T1::eval_type>,
       class D2 = std::decay_t<typename T2::eval_type>,
       class Algebra = common_algebra_type_t<D1, D2>,
-      class L = detail::pending_dimensions_list_t<D1, D2>>
+      class L = detail::pending_dimensions_list_t<D1, D2>,
+      class R =
+          typename tmp::convert_to_sequence_t<L, Algebra::template reified_blade>::type>
   constexpr auto
   operator()(expression_template::op<expression_template::multiplies, T1, T2>)
-      const -> std::bool_constant<(
-          tmp::convert_to_sequence_t<L, Algebra::template reified_blade>::type::
-              grade < D1::grade + D2::grade)>
+      const -> std::bool_constant<R::grade != D1::grade + D2::grade>
   {
     return {};
   }

geometry/src/interior_product.hpp

@@ -0,0 +1,63 @@
+#pragma once
+
+#include "geometry/expression_template.hpp"
+#include "geometry/src/common_algebra_type.hpp"
+#include "geometry/src/detail/multivector_product.hpp"
+#include "geometry/src/detail/pending_dimensions_list.hpp"
+#include "geometry/type_metaprogramming.hpp"
+
+#include <type_traits>
+
+namespace geometry {
+namespace detail {
+
+inline constexpr class
+{
+public:
+  template <
+      class T1,
+      class T2,
+      class D1 = std::decay_t<typename T1::eval_type>,
+      class D2 = std::decay_t<typename T2::eval_type>,
+      class Algebra = common_algebra_type_t<D1, D2>,
+      class L = detail::pending_dimensions_list_t<D1, D2>,
+      class R =
+          typename tmp::convert_to_sequence_t<L, Algebra::template reified_blade>::type,
+      int k = tmp::convert_to_sequence_t<
+          L,
+          Algebra::template reified_blade_coefficient>::value>
+  constexpr auto
+  operator()(expression_template::op<expression_template::multiplies, T1, T2>)
+      const -> std::bool_constant<
+          (R::grade !=
+           std::max(D1::grade, D2::grade) - std::min(D1::grade, D2::grade)) or
+          (k == 0)>
+  {
+    return {};
+  }
+} non_grade_contracting{};
+
+}  // namespace detail
+
+/// interior product
+///
+/// @{
+
+inline constexpr auto interior_product =  //
+    detail::multivector_product_with(     //
+        detail::non_grade_contracting,
+        [](auto algebra) {
+          return typename decltype(algebra)::template blade<>{};
+        });
+
+template <class T1, class T2, class A = common_algebra_type_t<T1, T2>>
+[[nodiscard]]
+constexpr auto
+operator|(const T1& x, const T2& y)
+{
+  return interior_product(x, y);
+}
+
+/// @}
+
+}  // namespace geometry

test/BUILD.bazel

@@ -79,3 +79,14 @@ cc_test(
         "@skytest",
     ],
 )
+
+cc_test(
+    name = "algebra_interior_product_test",
+    size = "small",
+    srcs = ["algebra_interior_product_test.cpp"],
+    deps = [
+        ":test_same",
+        "//:geometry",
+        "@skytest",
+    ],
+)

test/algebra_interior_product_test.cpp

@@ -0,0 +1,44 @@
+#include "geometry/geometry.hpp"
+#include "skytest/skytest.hpp"
+
+#include "test/test_same.hpp"
+
+#include <cstddef>
+#include <tuple>
+#include <utility>
+
+auto main() -> int
+{
+  using namespace ::skytest::literals;
+  using ::geometry::get;
+  using ::geometry::test::same;
+  using ::skytest::eq;
+  using ::skytest::expect;
+  using ::skytest::ne;
+
+  "line orthogonal to point (2D)"_ctest = [] {
+    using algebra = ::geometry::algebra<double, 2>;
+
+    // https://bivector.net/PROJECTIVE_GEOMETRIC_ALGEBRA.pdf
+    // section 7.1
+
+    constexpr auto x_hat = algebra::e<2, 0>;
+    constexpr auto y_hat = algebra::e<0, 1>;
+    constexpr auto E0 = algebra::e<1, 2>;
+
+    const auto P = 0 * x_hat + 0 * y_hat + E0;
+
+    // a line of the form
+    // a * x + b * y + c = 0
+    // is represented as
+    // a * e1 + b * e2 + c * e0 = 0
+
+    //using e0 = algebra::blade<0>;
+    constexpr auto x = algebra::e<1>;
+    constexpr auto y = algebra::e<2>;
+
+    return expect(eq(x, y | P) and
+                  eq(y, x | P));
+  };
+
+}