Supply option for enabling clang and gcc compiler profiling

public/gal/algebra.hpp

@@ -22,9 +22,6 @@ using width_t = uint32_t;
 // order of "0" here, by convention, refers to an infinite order.
 struct ind
 {
-    // NOTE: the tag of ~0 - 1 is RESERVED by the library for the dual unit. This ensures that the
-    // dual unit (which has a high chance of extinguishing the monomial it's a factor of) comes
-    // first in the factor ordering.
     width_t id    = ~0u;
     int degree    = 0;
     uint8_t order = 0;
@@ -34,12 +31,12 @@ struct ind
 // If all indeterminates are identified (ID != ~0ull), the ordering becomes a total order.
 [[nodiscard]] constexpr bool operator==(ind lhs, ind rhs) noexcept
 {
-    return lhs.id == rhs.id && lhs.degree == rhs.degree && lhs.id != ~0u && rhs.id != ~0u;
+    return lhs.id == rhs.id && lhs.degree == rhs.degree;
 }
 
 [[nodiscard]] constexpr bool operator!=(ind lhs, ind rhs) noexcept
 {
-    return lhs.id != rhs.id || lhs.degree != rhs.degree || lhs.id == ~0u || rhs.id == ~0u;
+    return lhs.id != rhs.id || lhs.degree != rhs.degree;
 }
 
 [[nodiscard]] constexpr bool operator<(ind lhs, ind rhs) noexcept
@@ -428,6 +425,12 @@ struct const_term_it
     }
 };
 
+template <width_t IndMax, width_t MonMax, width_t TermMax>
+struct dense_mv
+{
+    std::array<std::pair<term, std::array<std::pair<mon, std::array<ind, IndMax>>, MonMax>>, TermMax> data;
+};
+
 // Multivector representation, intended to be a compile-time representation
 // A := Algebra
 // I := Indeterminate capacity
@@ -502,7 +505,7 @@ struct mv
     }
 
     template <width_t I2, width_t M2, width_t T2>
-    [[nodiscard]] constexpr auto shrink() const noexcept
+    [[nodiscard]] constexpr auto resize() const noexcept
     {
         if constexpr (I == I2 && M == M2 && T == T2)
         {
@@ -528,6 +531,59 @@ struct mv
         }
     }
 
+    // Reshape the multivector to be fully dense in the arrays to simplify compilation of the final reduction
+    template <width_t IndMax, width_t MonMax, width_t TermMax>
+    [[nodiscard]] constexpr auto regularize() const noexcept
+    {
+        dense_mv<IndMax, MonMax, TermMax> out;
+
+        for (width_t i = 0; i != size.term; ++i)
+        {
+            auto const& term = terms[i];
+            auto& out_term   = out.data[i];
+            out_term.first   = term;
+
+            for (width_t j = term.mon_offset; j != term.mon_offset + term.count; ++j)
+            {
+                auto const& mon = mons[j];
+                auto& out_mon   = out_term.second[j - term.mon_offset];
+                out_mon.first   = mon;
+
+                for (width_t k = mon.ind_offset; k != mon.ind_offset + mon.count; ++k)
+                {
+                    out_mon.second[k - mon.ind_offset] = inds[k];
+                }
+            }
+        }
+
+        return out;
+    }
+
+    [[nodiscard]] constexpr mv_size extent() const noexcept
+    {
+        // Determine the maximum number of indeterminates across all monomials and the maximum number of monomials
+        // across all terms
+
+        width_t mon_count = 0;
+        width_t ind_count = 0;
+        for (auto term_it = cbegin(); term_it != cend(); ++term_it)
+        {
+            if (term_it->count > mon_count)
+            {
+                mon_count = term_it->count;
+            }
+
+            for (auto mon_it = term_it.cbegin(); mon_it != term_it.cend(); ++mon_it)
+            {
+                if (mon_it->count > ind_count)
+                {
+                    ind_count = mon_it->count;
+                }
+            }
+        }
+        return {ind_count, mon_count, size.term};
+    }
+
     [[nodiscard]] constexpr const_term_it cbegin() const noexcept
     {
         return {inds.data(), mons.data(), terms.data()};

public/gal/engine.hpp

@@ -78,7 +78,7 @@ namespace detail
 
     // Data := flattened array of inputs
     template <auto const& ie, typename F, typename A, typename Data, size_t... I>
-    [[nodiscard]] static auto compute_entity(Data const& data, std::index_sequence<I...>) noexcept
+    [[nodiscard]] constexpr static auto compute_entity_typed(Data const& data, std::index_sequence<I...>) noexcept
     {
         using entity_t = entity<A, F, ie.terms[I].element...>;
         // NOTE: this should be flattened to arrays of constant size ideally
@@ -114,6 +114,45 @@ namespace detail
             tree);
     }
 
+    template <auto const& dmv, typename F, typename A, typename Data, size_t... I>
+    [[nodiscard]] constexpr static auto compute_entity(Data const& data, std::index_sequence<I...>) noexcept
+    {
+        using entity_t = entity<A, F, dmv.data[I].first.element...>;
+        // NOTE: this should be flattened to arrays of constant size ideally
+        return std::apply(
+            [&data](auto&&... t) {
+                return entity_t{std::apply(
+                    [&](auto&&... m) {
+                        if constexpr (sizeof...(m) == 0)
+                        {
+                            return F{0};
+                        }
+                        else
+                        {
+                            return ((m.first.q.num == 0
+                                         ? 0
+                                         : (static_cast<F>(m.first.q)
+                                            * std::apply(
+                                                [&](auto&&... i) {
+                                                    if constexpr (sizeof...(i) == 0)
+                                                    {
+                                                        return F{1};
+                                                    }
+                                                    else
+                                                    {
+                                                        return ((i.id == ~0u ? 1 : ::gal::pow(*data[i.id], i.degree))
+                                                                * ...);
+                                                    }
+                                                },
+                                                m.second)))
+                                    + ...);
+                        }
+                    },
+                    t.second)...};
+            },
+            dmv.data);
+    }
+
     // The indeterminate value is either a pointer to an entity's value or an evaluated expression
     template <typename T>
     struct ind_value
@@ -123,7 +162,7 @@ namespace detail
             T const* pointer;
             T value;
         };
-        
+
         bool is_value;
 
         [[nodiscard]] constexpr T operator*() const noexcept
@@ -156,18 +195,37 @@ namespace detail
     }
 
     template <typename A, typename V, typename T, typename D>
-    [[nodiscard]] static auto finalize_entity(D const& data)
+    [[nodiscard]] static auto finalize_entity_typed(D const& data)
     {
         constexpr static auto reified = reify<T>();
         if constexpr (detail::uses_null_basis<A>)
         {
             constexpr static auto null_conversion = detail::to_null_basis(reified);
-            return compute_entity<null_conversion, V, A>(
-                data, std::make_index_sequence<null_conversion.size.term>());
+            return compute_entity_typed<null_conversion, V, A>(data, std::make_index_sequence<null_conversion.size.term>());
+        }
+        else
+        {
+            return compute_entity_typed<reified, V, A>(data, std::make_index_sequence<reified.size.term>());
+        }
+    }
+
+    template <typename A, typename V, typename T, typename D>
+    [[nodiscard]] static auto finalize_entity(D const& data)
+    {
+        constexpr auto reified = reify<T>();
+        if constexpr (detail::uses_null_basis<A>)
+        {
+            constexpr auto null_conversion = detail::to_null_basis(reified);
+            constexpr auto extent          = null_conversion.extent();
+            constexpr static auto regularized
+                = null_conversion.template regularize<extent.ind, extent.mon, extent.term>();
+            return compute_entity<regularized, V, A>(data, std::make_index_sequence<null_conversion.size.term>());
         }
         else
         {
-            return compute_entity<reified, V, A>(data, std::make_index_sequence<reified.size.term>());
+            constexpr auto extent             = reified.extent();
+            constexpr static auto regularized = reified.template regularize<extent.ind, extent.mon, extent.term>();
+            return compute_entity<regularized, V, A>(data, std::make_index_sequence<reified.size.term>());
         }
     }
 } // namespace detail
@@ -205,7 +263,7 @@ template <typename L, typename... Data>
     // an evaluated property
     if constexpr (detail::is_tuple_v<ie_result_t>)
     {
-        if constexpr (std::tuple_size_v<ie_result_t> > 0)
+        if constexpr (std::tuple_size_v<ie_result_t>> 0)
         {
             using value_t   = typename std::tuple_element_t<0, ie_result_t>::value_t;
             using algebra_t = typename std::tuple_element_t<0, ie_result_t>::algebra_t;
@@ -216,7 +274,7 @@ template <typename L, typename... Data>
             return std::apply(
                 [&](auto&&... args) {
                     return std::make_tuple(
-                        detail::finalize_entity<algebra_t, value_t, std::decay_t<decltype(args)>>(data)...);
+                        detail::finalize_entity_typed<algebra_t, value_t, std::decay_t<decltype(args)>>(data)...);
                 },
                 ie_result_t{});
         }
@@ -228,7 +286,7 @@ template <typename L, typename... Data>
 
         std::array<detail::ind_value<value_t>, (Data::ind_count() + ...)> data{};
         detail::fill(data.data(), input...);
-        return detail::finalize_entity<algebra_t, value_t, ie_result_t>(data);
+        return detail::finalize_entity_typed<algebra_t, value_t, ie_result_t>(data);
     }
 }
 } // namespace gal

public/gal/expression.hpp

@@ -247,7 +247,7 @@ template <typename exp_t>
         if constexpr (detail::uses_null_basis<typename exp_t::algebra_t>)
         {
             constexpr auto out = detail::to_natural_basis(exp_t::lhs);
-            return out.template shrink<out.size.ind, out.size.mon, out.size.term>();
+            return out.template resize<out.size.ind, out.size.mon, out.size.term>();
         }
         else
         {
@@ -282,7 +282,7 @@ template <typename exp_t>
     else if constexpr (exp_t::op == expr_op::extract)
     {
         constexpr auto out = detail::extract(reify<typename exp_t::lhs_t>(), exp_t::elements);
-        return out.template shrink<out.size.ind, out.size.mon, out.size.term>();
+        return out.template resize<out.size.ind, out.size.mon, out.size.term>();
     }
     else if constexpr (exp_t::op == expr_op::select)
     {
@@ -296,18 +296,18 @@ template <typename exp_t>
         if constexpr (exp_t::op == expr_op::sum)
         {
             constexpr auto out = detail::sum(lhs, rhs);
-            return out.template shrink<out.size.ind, out.size.mon, out.size.term>();
+            return out.template resize<out.size.ind, out.size.mon, out.size.term>();
         }
         else if constexpr (exp_t::op == expr_op::difference)
         {
             constexpr auto n_rhs = detail::negate(rhs);
             constexpr auto out   = detail::sum(lhs, n_rhs);
-            return out.template shrink<out.size.ind, out.size.mon, out.size.term>();
+            return out.template resize<out.size.ind, out.size.mon, out.size.term>();
         }
         else if constexpr (exp_t::op == expr_op::geometric)
         {
             constexpr auto out = detail::product(typename decltype(lhs)::algebra_t::geometric{}, lhs, rhs);
-            return out.template shrink<out.size.ind, out.size.mon, out.size.term>();
+            return out.template resize<out.size.ind, out.size.mon, out.size.term>();
         }
         else if constexpr (exp_t::op == expr_op::sandwich)
         {
@@ -316,31 +316,31 @@ template <typename exp_t>
             constexpr auto temp  = detail::product(typename decltype(lhs)::algebra_t::geometric{}, lhs, rhs_reverse);
             constexpr auto temp2 = detail::product(typename decltype(lhs)::algebra_t::geometric{},
                                                    rhs,
-                                                   temp.template shrink<temp.size.ind, temp.size.mon, temp.size.term>());
-            return temp2.template shrink<temp2.size.ind, temp2.size.mon, temp2.size.term>();
+                                                   temp.template resize<temp.size.ind, temp.size.mon, temp.size.term>());
+            return temp2.template resize<temp2.size.ind, temp2.size.mon, temp2.size.term>();
         }
         else if constexpr (exp_t::op == expr_op::exterior)
         {
             constexpr auto out = detail::product(typename decltype(lhs)::algebra_t::exterior{}, lhs, rhs);
-            return out.template shrink<out.size.ind, out.size.mon, out.size.term>();
+            return out.template resize<out.size.ind, out.size.mon, out.size.term>();
         }
         else if constexpr (exp_t::op == expr_op::regressive)
         {
             // TODO: check if both lhs and rhs are dual
             constexpr auto lhs_dual = detail::poincare_dual(lhs);
             constexpr auto rhs_dual = detail::poincare_dual(rhs);
             constexpr auto out = detail::product(typename decltype(lhs)::algebra_t::exterior{}, lhs_dual, rhs_dual);
-            return detail::poincare_dual(out.template shrink<out.size.ind, out.size.mon, out.size.term>());
+            return detail::poincare_dual(out.template resize<out.size.ind, out.size.mon, out.size.term>());
         }
         else if constexpr (exp_t::op == expr_op::contract)
         {
             constexpr auto out = detail::product(typename decltype(lhs)::algebra_t::contract{}, lhs, rhs);
-            return out.template shrink<out.size.ind, out.size.mon, out.size.term>();
+            return out.template resize<out.size.ind, out.size.mon, out.size.term>();
         }
         else if constexpr (exp_t::op == expr_op::symmetric_inner)
         {
             constexpr auto out = detail::product(typename decltype(lhs)::algebra_t::symmetric_inner{}, lhs, rhs);
-            return out.template shrink<out.size.ind, out.size.mon, out.size.term>();
+            return out.template resize<out.size.ind, out.size.mon, out.size.term>();
         }
     }
 }

public/gal/numeric.hpp

@@ -8,20 +8,19 @@
 namespace gal
 {
 // right-to-left binary exponentiation
-template <typename T, int exponent>
-[[nodiscard]] constexpr T pow(T s, std::integral_constant<int, exponent>) noexcept
+template <typename T>
+[[nodiscard]] constexpr T pow(T s, int e) noexcept
 {
-    if constexpr (exponent < 0)
+    if (e < 0)
     {
-        return T{1} / pow(s, std::integral_constant<int, -exponent>{});
+        return T{1} / pow(s, -e);
     }
-    else if constexpr (exponent == 1)
+    else if (e == 1)
     {
         return s;
     }
-    else 
+    else
     {
-        int e = exponent;
         T temp1{1};
         T temp2{s};
         while (e > 1)
@@ -98,7 +97,7 @@ template <typename T>
 // overflows
 struct rat
 {
-    int num = 1;
+    int num = 0;
     int den = 1;
 
     [[nodiscard]] constexpr bool is_zero() const noexcept

public/gal/pga.hpp

@@ -152,6 +152,7 @@ namespace pga
 
         template <uint8_t... E>
         constexpr point(entity<pga_algebra, T, E...> in) noexcept
+            : data{}
         {
             auto input = in.template select<0b111, 0b1011, 0b1101, 0b1110>();
             auto w_inv = T{1} / input[3];
@@ -231,6 +232,7 @@ namespace pga
 
         template <uint8_t... E>
         constexpr vector(entity<pga_algebra, T, E...> in) noexcept
+            : data{}
         {
             auto input = in.template select<0b111, 0b1011, 0b1101>();
             z = -input[0];

test/CMakeLists.txt

@@ -16,7 +16,7 @@ add_executable(gal_test
     test_cga.cpp
     test_ega.cpp
     test_pga.cpp)
-    #test_ik.cpp)
+    # test_ik.cpp)
 
 target_link_libraries(gal_test PRIVATE gal doctest)
 target_compile_definitions(gal_test PRIVATE

test/test_ega.cpp

@@ -3,6 +3,7 @@
 #include <doctest/doctest.h>
 #include <gal/ega.hpp>
 #include <gal/engine.hpp>
+#include <gal/format.hpp>
 
 #include <cstdio>
 
@@ -118,6 +119,10 @@ TEST_CASE("rotors")
             r1,
             r2,
             v1);
+        auto r = gal::evaluate<rotor<float>, rotor<float>, vector<float>>{}.debug([](auto r1, auto r2, auto v1) {
+            auto r = r1 * r2;
+            return v1 % r;
+        });
         CHECK_EQ(rotated.x, doctest::Approx(0.0));
         CHECK_EQ(rotated.y, doctest::Approx(1.0));
         CHECK_EQ(rotated.z, doctest::Approx(0.0));