diff --git a/src/analysis/Polynomial.hpp b/src/analysis/Polynomial.hpp index c768d2d6757545bdddfdd4ceb8797f2e8630c607..887fbf12ae2672d379016f5a0b3f1a4ecb0eb7a4 100644 --- a/src/analysis/Polynomial.hpp +++ b/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 diff --git a/tests/test_Polynomial.cpp b/tests/test_Polynomial.cpp index 1d227564f97810f40a855e46df436336e5c9bf34..fd1f07a2c24cb8113de96feb5049f94682e3f266 100644 --- a/tests/test_Polynomial.cpp +++ b/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})); } }