Implement a template loop for compose

src/analysis/Polynomial.hpp

@@ -199,6 +199,36 @@ class Polynomial
     }
     return P;
   }
+  template <size_t M, size_t I>
+  PUGS_INLINE constexpr Polynomial<M * N>
+  power(const Polynomial<M>& Q) const
+  {
+    return coefficients()[I] * Q.template pow<N>(I);
+  }
+
+  template <size_t M, size_t... I>
+  PUGS_INLINE constexpr Polynomial<M * N>
+  compose2(const Polynomial<M>& Q, std::index_sequence<I...>) const
+  {
+    Polynomial<M * N> P;
+    P.coefficients() = zero;
+    P                = (power<M, I>(Q) + ...);
+    return P;
+  }
+
+  template <size_t M>
+  PUGS_INLINE constexpr Polynomial<M * N>
+  compose2(const Polynomial<M>& Q) const
+  {
+    Polynomial<M * N> P;
+    P.coefficients()    = zero;
+    using IndexSequence = std::make_index_sequence<N + 1>;
+    return compose2<M>(Q, IndexSequence{});
+    //    for (size_t i = 0; i <= N; ++i) {
+    //      P += Q.template pow<N>(i) * coefficients()[i];
+    //    }
+    return P;
+  }
 
   template <size_t M>
   PUGS_INLINE constexpr Polynomial<M * N>
@@ -408,8 +438,8 @@ class Polynomial
   }
 
   PUGS_INLINE
-  constexpr Polynomial<N>
-  lagrangeToCanonical(const TinyVector<N + 1> lagrange_coefs, TinyVector<N + 1, Polynomial<N>>& basis)
+  constexpr friend Polynomial<N>
+  lagrangeToCanonical(const TinyVector<N + 1> lagrange_coefs, const TinyVector<N + 1, Polynomial<N>>& basis)
   {
     Polynomial<N> P(zero);
     //    lagrangeBasis({0, 0, 0}, basis);
@@ -428,4 +458,18 @@ class Polynomial
   constexpr Polynomial() noexcept = default;
   ~Polynomial()                   = default;
 };
+
+template <size_t N>
+PUGS_INLINE constexpr TinyVector<N, Polynomial<N - 1>>
+lagrangeBasis(const TinyVector<N>& zeros)
+{
+  static_assert(N >= 1, "invalid degree");
+  TinyVector<N, Polynomial<N - 1>> basis;
+  Polynomial<N - 1> lj;
+  for (size_t j = 0; j < N; ++j) {
+    lagrangePolynomial(zeros, j, basis[j]);
+  }
+  return basis;
+}
+
 #endif   // POLYNOMIAL_HPP

tests/test_Polynomial.cpp

@@ -174,6 +174,7 @@ TEST_CASE("Polynomial", "[analysis]")
     Polynomial<2> R({2, -1, 3});
     Polynomial<2> S({1, 2, 2});
     REQUIRE(P.compose(Q) == Polynomial<2>({2, -6, 12}));
+    REQUIRE(P.compose2(Q) == Polynomial<2>({2, -6, 12}));
     REQUIRE(R(S) == Polynomial<4>({4, 10, 22, 24, 12}));
   }
 
@@ -187,12 +188,11 @@ TEST_CASE("Polynomial", "[analysis]")
     Polynomial<2> R;
     lagrangePolynomial({-1, 0, 1}, 0, R);
     REQUIRE(R == P);
-    TinyVector<3, Polynomial<2>> basis;
-    //    lagrangeBasis<2>({-1,0,1},basis);
-    basis[0] = Polynomial<2>({0, -0.5, 0.5});
-    basis[1] = Polynomial<2>({1, 0, -1});
-    basis[2] = Polynomial<2>({0, 0.5, 0.5});
-    // Q = lagrangeToCanonical({1,0,1},basis);
-    // REQUIRE(lagrangeToCanonical<2>({1, 0, 1}, basis) == Polynomial<2>({0, 0, 1}));
+    const TinyVector<3, Polynomial<2>> basis = lagrangeBasis(TinyVector<3>{-1, 0, 1});
+    //    lagrangeBasis({-1, 0, 1}, basis);
+    // basis[0] = Polynomial<2>({0, -0.5, 0.5});
+    // basis[1] = Polynomial<2>({1, 0, -1});
+    // basis[2] = Polynomial<2>({0, 0.5, 0.5});
+    REQUIRE(lagrangeToCanonical({1, 0, 1}, basis) == Polynomial<2>({0, 0, 1}));
   }
 }