define canonical form for blade types

Change-Id: If2396afbe3d57b84214445eb6728e660dc6c313a

rigid_geometric_algebra/BUILD.bazel

@@ -5,6 +5,8 @@ cc_library(
     srcs = [
         "algebra.hpp",
         "blade.hpp",
+        "detail/sorted_dimensions.hpp",
+        "detail/swaps_to_sorted_dimensions.hpp",
         "detail/unique_dimensions.hpp",
     ],
     hdrs = ["rigid_geometric_algebra.hpp"],

rigid_geometric_algebra/blade.hpp

@@ -1,9 +1,14 @@
 #pragma once
 
 #include "rigid_geometric_algebra/algebra.hpp"
+#include "rigid_geometric_algebra/detail/sorted_dimensions.hpp"
+#include "rigid_geometric_algebra/detail/swaps_to_sorted_dimensions.hpp"
 #include "rigid_geometric_algebra/detail/unique_dimensions.hpp"
 
 #include <cstddef>
+#include <functional>
+#include <type_traits>
+#include <utility>
 
 namespace rigid_geometric_algebra {
 
@@ -35,15 +40,70 @@ template <class A, std::size_t... Is>
 struct blade
 {
   static_assert(((Is < algebra_dimension_v<A>) and ...));
-  static_assert(detail::unique_dimensions(Is...));
+  static_assert(detail::unique_dimensions<Is...>());
 
+  /// number of factors
+  ///
   static constexpr auto grade = sizeof...(Is);
 
+  /// blade scalar type
+  ///
   using value_type = algebra_field_t<A>;
 
+  /// blade type with indices in canonical form
+  ///
+  using canonical_type =
+      decltype([]<std::size_t... Js>(
+                   std::index_sequence<Js...>) -> blade<A, Js...> {
+        return {};
+      }(detail::sorted_dimensions<Is...>{}));
+
+  /// blade scalar coefficient
+  ///
   value_type coefficient{};
 
+  /// negation
+  ///
+  template <class Self>
+  constexpr auto operator-(this Self&& self) -> blade
+  {
+    return {-std::forward<Self>(self).coefficient};
+  }
+
+  /// obtain the blade with indices in canonical form
+  ///
+  /// Returns same blade expressed in canonical form - i.e. with indices in
+  /// increasing order. The blade coefficient is negated if sorting the indices
+  /// requires an odd number of swaps.
+  ///
+  template <class Self>
+  constexpr auto canonical(this Self&& self) -> canonical_type
+  {
+    constexpr auto even = [](std::size_t value) { return value % 2UZ == 0UZ; };
+    using maybe_negate = std::conditional_t<
+        even(detail::swaps_to_sorted_dimensions<Is...>()),
+        std::identity,
+        std::negate<>>;
+
+    return {maybe_negate{}(std::forward<Self>(self).coefficient)};
+  }
+
+  /// equality comparison
+  ///
+  /// @{
+
   friend auto operator==(const blade&, const blade&) -> bool = default;
+
+  template <std::size_t... Js>
+    requires (std::is_same_v<typename blade::canonical_type,
+                             typename blade<A, Js...>::canonical_type>)
+  friend constexpr auto
+  operator==(const blade& lhs, const blade<A, Js...>& rhs) -> bool
+  {
+    return lhs.canonical() == rhs.canonical();
+  }
+
+  /// @}
 };
 
 }  // namespace rigid_geometric_algebra

rigid_geometric_algebra/detail/sorted_dimensions.hpp

@@ -0,0 +1,51 @@
+#pragma once
+
+#include "rigid_geometric_algebra/detail/unique_dimensions.hpp"
+
+#include <cstddef>
+#include <tuple>
+#include <type_traits>
+
+namespace rigid_geometric_algebra {
+namespace detail {
+
+inline constexpr class
+{
+  template <std::size_t N>
+  using constant = std::integral_constant<std::size_t, N>;
+
+  template <std::size_t N, class... Is>
+  static constexpr auto make_sequence_from(Is... is)
+  {
+    if constexpr (((N > Is::value) and ...)) {
+      return std::tuple{};
+    } else if constexpr (((N == Is::value) or ...)) {
+      return std::tuple_cat(
+          std::tuple<constant<N>>{}, make_sequence_from<N + 1>(is...));
+    } else {
+      return make_sequence_from<N + 1>(is...);
+    }
+  }
+
+public:
+  constexpr auto operator()() const noexcept { return std::index_sequence<>{}; }
+
+  template <class... Is>
+  constexpr auto operator()(Is... is) const noexcept
+  {
+    return []<class... Js>(
+               std::tuple<Js...>) -> std::index_sequence<Js::value...> {
+      return {};
+    }(make_sequence_from<0>(is...));
+  }
+
+} sort_dimensions_impl{};
+
+template <std::size_t... Is>
+  requires (unique_dimensions<Is...>())
+struct sorted_dimensions
+    : decltype(sort_dimensions_impl(
+          std::integral_constant<std::size_t, Is>{}...)){};
+
+}  // namespace detail
+}  // namespace rigid_geometric_algebra

rigid_geometric_algebra/detail/swaps_to_sorted_dimensions.hpp

@@ -0,0 +1,67 @@
+#pragma once
+
+#include "rigid_geometric_algebra/detail/sorted_dimensions.hpp"
+#include "rigid_geometric_algebra/detail/unique_dimensions.hpp"
+
+#include <algorithm>
+#include <array>
+#include <cassert>
+#include <cstddef>
+#include <iterator>
+#include <ranges>
+#include <type_traits>
+
+namespace rigid_geometric_algebra {
+namespace detail {
+
+template <std::size_t... Is>
+class swaps_to_sorted_dimensions_fn
+{
+  template <
+      std::ranges::random_access_range R1,
+      std::ranges::random_access_range R2>
+    requires std::ranges::view<R1> and std::ranges::view<R2>
+  static constexpr auto impl(R1 unsorted, R2 sorted) -> std::size_t
+  {
+    assert(unsorted.size() == sorted.size());
+
+    if (unsorted.empty()) {
+      return 0UZ;
+    }
+
+    if (unsorted.size() == 1) {
+      assert(unsorted.front() == sorted.front());
+      return 0Uz;
+    }
+
+    auto it = std::ranges::find(unsorted, sorted.front());
+    assert(it != unsorted.end());
+    std::iter_swap(it, unsorted.begin());
+
+    return std::size_t(it != unsorted.begin()) +
+           impl(unsorted | std::views::drop(1), sorted | std::views::drop(1));
+  }
+
+public:
+  constexpr auto operator()() const noexcept -> std::size_t
+  {
+    assert(unique_dimensions<Is...>());
+
+    constexpr auto into_array =
+        []<std::size_t... Js>(std::index_sequence<Js...>) {
+          return std::array<std::size_t, sizeof...(Js)>{Js...};
+        };
+
+    auto unsorted = into_array(std::index_sequence<Is...>{});
+    const auto sorted = into_array(sorted_dimensions<Is...>{});
+
+    return impl(std::ranges::subrange(unsorted), std::ranges::subrange(sorted));
+  }
+};
+
+template <std::size_t... Is>
+inline constexpr auto swaps_to_sorted_dimensions =
+    swaps_to_sorted_dimensions_fn<Is...>{};
+
+}  // namespace detail
+}  // namespace rigid_geometric_algebra

rigid_geometric_algebra/detail/unique_dimensions.hpp

@@ -2,22 +2,28 @@
 
 #include <algorithm>
 #include <array>
-#include <concepts>
 #include <cstddef>
 
 namespace rigid_geometric_algebra {
 namespace detail {
 
-template <std::same_as<std::size_t>... Ts>
-constexpr auto unique_dimensions(Ts... is)
+template <std::size_t... Is>
+class unique_dimensions_fn
 {
-  auto dims = std::array<std::size_t, sizeof...(is)>{is...};
+public:
+  constexpr auto operator()() const noexcept -> bool
+  {
+    auto dims = std::array<std::size_t, sizeof...(Is)>{Is...};
 
-  std::ranges::sort(dims);
+    std::ranges::sort(dims);
 
-  return dims.cend() ==
-         std::ranges::adjacent_find(dims, std::ranges::equal_to{});
-}
+    return dims.cend() ==
+           std::ranges::adjacent_find(dims, std::ranges::equal_to{});
+  }
+};
+
+template <std::size_t... Is>
+inline constexpr auto unique_dimensions = unique_dimensions_fn<Is...>{};
 
 }  // namespace detail
 }  // namespace rigid_geometric_algebra

test/blade_test.cpp

@@ -10,7 +10,7 @@ auto main() -> int
 
   using rga = ::rigid_geometric_algebra::algebra<double, 2>;
 
-  "blade constructible"_test = [] {
+  "constructible"_test = [] {
     const auto a = rga::blade<>{};
     const auto b = rga::blade<0>{0};
     const auto c = rga::blade<1>{1};
@@ -21,15 +21,15 @@ auto main() -> int
         eq(1, c.coefficient) and eq(2, d.coefficient));
   };
 
-  "blade comparable"_test = [] {
+  "comparable"_test = [] {
     const auto a = rga::blade<>{};
     const auto b = rga::blade<>{};
     const auto c = rga::blade<>{1};
 
     return expect(eq(a, b) and ne(a, c) and ne(b, c));
   };
 
-  "blade grade"_test = [] {
+  "grade"_test = [] {
     return expect(
         eq(0, rga::blade<>::grade) and      //
         eq(1, rga::blade<0>::grade) and     //
@@ -40,4 +40,31 @@ auto main() -> int
         eq(2, rga::blade<0, 2>::grade) and  //
         eq(3, rga::blade<0, 1, 2>::grade));
   };
+
+  "negatable"_test = [] {
+    return expect(eq(rga::blade<0, 1>{1}, -rga::blade<0, 1>{-1}));
+  };
+
+  "canonical form"_test = [] {
+    return expect(
+        eq(rga::blade<>{}, rga::blade<>{}.canonical()) and
+        eq(rga::blade<0>{}, rga::blade<0>{}.canonical()) and
+        eq(rga::blade<1>{}, rga::blade<1>{}.canonical()) and
+        eq(rga::blade<2>{}, rga::blade<2>{}.canonical()) and
+        eq(rga::blade<0, 1>{}, rga::blade<0, 1>{}.canonical()) and
+        eq(rga::blade<0, 1>{}, rga::blade<1, 0>{}.canonical()) and
+        eq(rga::blade<0, 2>{}, rga::blade<0, 2>{}.canonical()) and
+        eq(rga::blade<0, 2>{}, rga::blade<2, 0>{}.canonical()) and
+        eq(rga::blade<0, 1, 2>{}, rga::blade<0, 2, 1>{}.canonical()) and
+        eq(rga::blade<0, 1, 2>{}, rga::blade<1, 2, 0>{}.canonical()) and
+        eq(rga::blade<0, 1, 2>{}, rga::blade<0, 1, 2>{}.canonical()));
+  };
+
+  "sign conversion for canonical form"_test = [] {
+    return expect(
+        eq(rga::blade<0, 1>{1}, -rga::blade<1, 0>{1}) and
+        eq(rga::blade<0, 1, 2>{1}, rga::blade<1, 2, 0>{1}) and
+        eq(rga::blade<0, 1, 2>{1}, -rga::blade<0, 2, 1>{1}) and
+        eq(rga::blade<0, 1, 2>{1}, -rga::blade<2, 1, 0>{1}));
+  };
 }