diff --git a/src/analysis/Polynomial.hpp b/src/analysis/Polynomial.hpp index e43c4535221682d6b2e5c96b5492ad6d96f780c7..9014679bfbe861078cb35802a839f5cf1f2972e9 100644 --- a/src/analysis/Polynomial.hpp +++ b/src/analysis/Polynomial.hpp @@ -9,7 +9,7 @@ class Polynomial private: using Coefficients = TinyVector<N + 1, double>; Coefficients m_coefficients; - static_assert((N >= 0), "Polynomial degree must be non-negative"); + static_assert(N >= 0, "Polynomial degree must be non-negative"); public: PUGS_INLINE @@ -18,6 +18,7 @@ class Polynomial { return N; } + PUGS_INLINE constexpr size_t realDegree() const @@ -29,37 +30,32 @@ class Polynomial } return 0; } + PUGS_INLINE - constexpr const TinyVector<N + 1>& + constexpr const Coefficients& coefficients() const { return m_coefficients; } PUGS_INLINE - constexpr TinyVector<N + 1>& + constexpr Coefficients& coefficients() { return m_coefficients; } - template <size_t M> PUGS_INLINE constexpr bool - operator==(const Polynomial<M>& q) const + operator==(const Polynomial& q) const { - if (M != N) - return false; - if (m_coefficients != q.coefficients()) { - return false; - } - return true; + return m_coefficients == q.m_coefficients; } PUGS_INLINE constexpr bool - operator!=(const Polynomial<N>& q) const + operator!=(const Polynomial& q) const { - return not this->operator==(q); + return not(*this == q); } template <size_t M> @@ -81,12 +77,10 @@ class Polynomial return P; } - PUGS_INLINE constexpr Polynomial<N> + PUGS_INLINE constexpr Polynomial operator-() const { - Polynomial<N> P; - P.coefficients() = -coefficients(); - return P; + return Polynomial{-m_coefficients}; } template <size_t M> @@ -109,7 +103,7 @@ class Polynomial } template <size_t M> - PUGS_INLINE constexpr Polynomial<N>& + PUGS_INLINE constexpr Polynomial& operator=(const Polynomial<M>& Q) { coefficients() = zero; @@ -124,7 +118,7 @@ class Polynomial } template <size_t M> - PUGS_INLINE constexpr Polynomial<N>& + PUGS_INLINE constexpr Polynomial& operator+=(const Polynomial<M>& Q) { static_assert(N >= M, "Polynomial degree to small in affectation addition"); @@ -135,7 +129,7 @@ class Polynomial } template <size_t M> - PUGS_INLINE constexpr Polynomial<N>& + PUGS_INLINE constexpr Polynomial& operator-=(const Polynomial<M>& Q) { static_assert(N >= M, "Polynomial degree to small in affectation addition"); @@ -160,14 +154,14 @@ class Polynomial } template <size_t M> - PUGS_INLINE constexpr Polynomial<N>& + PUGS_INLINE constexpr Polynomial& operator*=(const Polynomial<M>& Q) { static_assert(N >= M, "Degree to small in affectation product"); for (size_t i = N - M + 1; i <= N; ++i) { Assert(coefficients()[i] == 0, "Degree of affectation product greater than the degree of input polynomial"); } - Polynomial<N> P(zero); + Polynomial P(zero); for (size_t i = 0; i <= N - M; ++i) { for (size_t j = 0; j <= M; ++j) { P.coefficients()[i + j] += coefficients()[i] * Q.coefficients()[j]; @@ -178,12 +172,10 @@ class Polynomial } PUGS_INLINE - constexpr Polynomial<N> + constexpr Polynomial operator*(const double& lambda) const { - TinyVector<N + 1> mult_coefs = lambda * m_coefficients; - Polynomial<N> M(mult_coefs); - return M; + return Polynomial(lambda * m_coefficients); } template <size_t M> @@ -268,36 +260,22 @@ class Polynomial } PUGS_INLINE - constexpr friend Polynomial<N> - operator*(const double& lambda, const Polynomial<N> P) + constexpr friend Polynomial + operator*(const double& lambda, const Polynomial& P) { return P * lambda; } - // evaluation using Horner's method https://en.wikipedia.org/wiki/Horner's_method - PUGS_INLINE - constexpr double - evaluate(const double& x) const - { - TinyVector<N + 1> coefs = this->coefficients(); - double bcoef = coefs[N]; - for (size_t i = N; i > 0; --i) { - bcoef *= x; - bcoef += coefs[i - 1]; - } - return bcoef; - } PUGS_INLINE constexpr double - operator()(const double x) const + operator()(double x) const { - TinyVector<N + 1> coefs = this->coefficients(); - double bcoef = coefs[N]; + double p_x = m_coefficients[N]; for (size_t i = N; i > 0; --i) { - bcoef *= x; - bcoef += coefs[i - 1]; + p_x *= x; + p_x += m_coefficients[i - 1]; } - return bcoef; + return p_x; } template <size_t M> @@ -321,7 +299,7 @@ class Polynomial PUGS_INLINE constexpr friend Polynomial<N + 1> - primitive(const Polynomial<N>& P) + primitive(const Polynomial& P) { TinyVector<N + 2> coefs; for (size_t i = 0; i < N + 1; ++i) { @@ -333,7 +311,7 @@ class Polynomial PUGS_INLINE constexpr friend std::ostream& - operator<<(std::ostream& os, const Polynomial<N>& P) + operator<<(std::ostream& os, const Polynomial& P) { // os << "P(x) = "; bool all_coef_zero = true; @@ -412,16 +390,36 @@ class Polynomial return P; } - PUGS_INLINE constexpr Polynomial(const TinyVector<N + 1>& coefficients) noexcept : m_coefficients{coefficients} {} + PUGS_INLINE constexpr Polynomial& operator=(const Polynomial&) = default; + PUGS_INLINE constexpr Polynomial& operator=(Polynomial&&) = default; + + PUGS_INLINE + constexpr Polynomial(const TinyVector<N + 1>& coefficients) noexcept : m_coefficients(coefficients) {} PUGS_INLINE - constexpr Polynomial(TinyVector<N + 1>&& coefficients) noexcept : m_coefficients{coefficients} {} + constexpr Polynomial(const Polynomial&) noexcept = default; + + PUGS_INLINE + constexpr Polynomial(Polynomial&&) noexcept = default; + + template <typename... T> + explicit PUGS_INLINE constexpr Polynomial(T&&... coefficients) noexcept : m_coefficients(coefficients...) + { + // static_assert(sizeof...(T) == N + 1, "invalid number of parameters"); + // static_assert((std::is_convertible_v<T, double> and ...), "arguments must be convertible to double"); + } PUGS_INLINE constexpr Polynomial() noexcept = default; - ~Polynomial() = default; + + ~Polynomial() = default; }; +// template <size_t N> +// template <> +// PUGS_INLINE constexpr Polynomial(TinyVector<N + 1>&& coefficients) noexcept : m_coefficients{coefficients} +// {} + template <size_t N> PUGS_INLINE constexpr Polynomial<N> lagrangePolynomial(const TinyVector<N + 1> zeros, const size_t k); diff --git a/src/analysis/PolynomialBasis.hpp b/src/analysis/PolynomialBasis.hpp index cf64978b11c21252dfa6b6e98382806e71703306..a0b07ecafd27bbdc24dc41b4f210a2861cac46c1 100644 --- a/src/analysis/PolynomialBasis.hpp +++ b/src/analysis/PolynomialBasis.hpp @@ -38,8 +38,8 @@ class PolynomialBasis { TinyVector<N + 1> coefficients(zero); elements()[0] = Polynomial<N>(coefficients); - elements()[0] += Polynomial<0>(TinyVector<1>{1}); - Polynomial<1> unit(TinyVector<2>{-shift, 1}); + elements()[0] += Polynomial<0>(1); + Polynomial<1> unit(-shift, 1); for (size_t i = 1; i <= N; i++) { elements()[i] = elements()[i - 1] * unit; } @@ -51,7 +51,7 @@ class PolynomialBasis _buildLagrangeBasis(const TinyVector<N + 1>& zeros) { if constexpr (N == 0) { - elements()[0] = Polynomial<0>(TinyVector<1>{1}); + elements()[0] = Polynomial<0>(1); return *this; } else { for (size_t i = 0; i < N; ++i) { @@ -65,7 +65,7 @@ class PolynomialBasis if (i == j) continue; double adim = 1. / (zeros[i] - zeros[j]); - elements()[i] *= Polynomial<1>{TinyVector<2>{-zeros[j] * adim, adim}}; + elements()[i] *= Polynomial<1>{-zeros[j] * adim, adim}; } } return *this; @@ -128,7 +128,7 @@ class PolynomialBasis PUGS_INLINE constexpr PolynomialBasis<N>& - build(BasisType basis_type, const double& shift = 0, const TinyVector<N + 1>& zeros = TinyVector<N + 1>(zero)) + build(BasisType basis_type, const double& shift = 0, const TinyVector<N + 1>& zeros = zero) { type() = basis_type; switch (basis_type) { diff --git a/tests/test_Polynomial.cpp b/tests/test_Polynomial.cpp index 232ffcca02f445314d14b3418bc14924c5d19db8..e9f92dbea15b854e902c76544cdf50a137fccbb3 100644 --- a/tests/test_Polynomial.cpp +++ b/tests/test_Polynomial.cpp @@ -22,75 +22,83 @@ TEST_CASE("Polynomial", "[analysis]") { SECTION("construction") { - REQUIRE_NOTHROW(Polynomial<2>{TinyVector<3>{2, 3, 4}}); + REQUIRE_NOTHROW(Polynomial<2>(2, 3, 4)); } + SECTION("degree") { - Polynomial<2> P(TinyVector<3>{2, 3, 4}); + Polynomial<2> P(2, 3, 4); REQUIRE(P.degree() == 2); } + SECTION("equality") { - Polynomial<2> P(TinyVector<3>{2, 3, 4}); - Polynomial<2> Q(TinyVector<3>{2, 3, 4}); - Polynomial<2> S(TinyVector<3>{2, 3, 5}); + Polynomial<2> P(2, 3, 4); + Polynomial<2> Q(2, 3, 4); + Polynomial<2> S(2, 3, 5); REQUIRE(P == Q); REQUIRE(P != S); } + SECTION("addition") { - Polynomial<2> P(TinyVector<3>{2, 3, 4}); - Polynomial<2> Q(TinyVector<3>{-1, -3, 2}); - Polynomial<2> S(TinyVector<3>{1, 0, 6}); - Polynomial<3> T(TinyVector<4>{0, 3, 1, -2}); - Polynomial<3> U(TinyVector<4>{2, 6, 5, -2}); + Polynomial<2> P(2, 3, 4); + Polynomial<2> Q(-1, -3, 2); + Polynomial<2> S(1, 0, 6); + Polynomial<3> T(0, 3, 1, -2); + Polynomial<3> U(2, 6, 5, -2); REQUIRE(S == (P + Q)); REQUIRE((T + P) == U); } + SECTION("opposed") { - Polynomial<2> P(TinyVector<3>{2, 3, 4}); + Polynomial<2> P(2, 3, 4); Polynomial<2> Q = -P; - REQUIRE(Q == Polynomial<2>(TinyVector<3>{-2, -3, -4})); + REQUIRE(Q == Polynomial<2>(-2, -3, -4)); } + SECTION("difference") { - Polynomial<2> P(TinyVector<3>{2, 3, 4}); - Polynomial<2> Q(TinyVector<3>{3, 4, 5}); - Polynomial<2> D(TinyVector<3>{-1, -1, -1}); + Polynomial<2> P(2, 3, 4); + Polynomial<2> Q(3, 4, 5); + Polynomial<2> D(-1, -1, -1); REQUIRE(D == (P - Q)); - Polynomial<3> R(TinyVector<4>{2, 3, 4, 1}); + Polynomial<3> R(2, 3, 4, 1); REQUIRE(D == (P - Q)); - REQUIRE((P - R) == Polynomial<3>(TinyVector<4>{0, 0, 0, -1})); + REQUIRE((P - R) == Polynomial<3>{0, 0, 0, -1}); R -= P; - REQUIRE(R == Polynomial<3>(TinyVector<4>{0, 0, 0, 1})); + REQUIRE(R == Polynomial<3>(0, 0, 0, 1)); } + SECTION("product_by_scalar") { - Polynomial<2> P(TinyVector<3>{2, 3, 4}); - Polynomial<2> M(TinyVector<3>{6, 9, 12}); + Polynomial<2> P(2, 3, 4); + Polynomial<2> M(6, 9, 12); REQUIRE(M == (P * 3)); REQUIRE(M == (3 * P)); } + SECTION("product") { - Polynomial<2> P(TinyVector<3>{2, 3, 4}); - Polynomial<3> Q(TinyVector<4>{1, 2, -1, 1}); + Polynomial<2> P(2, 3, 4); + Polynomial<3> Q(1, 2, -1, 1); Polynomial<4> R; Polynomial<5> S; R = P; S = P; S *= Q; - REQUIRE(Polynomial<5>(TinyVector<6>{2, 7, 8, 7, -1, 4}) == (P * Q)); - REQUIRE(Polynomial<5>(TinyVector<6>{2, 7, 8, 7, -1, 4}) == S); + REQUIRE(Polynomial<5>(2, 7, 8, 7, -1, 4) == (P * Q)); + REQUIRE(Polynomial<5>(2, 7, 8, 7, -1, 4) == S); // REQUIRE_THROWS_AS(R *= Q, AssertError); } + SECTION("divide") { - Polynomial<2> P(TinyVector<3>{1, 0, 1}); - Polynomial<1> Q(TinyVector<2>{0, 1}); - Polynomial<1> Q1(TinyVector<2>{0, 1}); + Polynomial<2> P(1, 0, 1); + Polynomial<1> Q(0, 1); + Polynomial<1> Q1(0, 1); Polynomial<2> R; Polynomial<2> S; @@ -99,29 +107,29 @@ TEST_CASE("Polynomial", "[analysis]") REQUIRE(Q1.realDegree() == 1); divide(P, Q1, R, S); - REQUIRE(Polynomial<2>(TinyVector<3>{1, 0, 0}) == S); - REQUIRE(Polynomial<2>(TinyVector<3>{0, 1, 0}) == R); + REQUIRE(Polynomial<2>{1, 0, 0} == S); + REQUIRE(Polynomial<2>{0, 1, 0} == R); } + SECTION("evaluation") { - Polynomial<2> P(TinyVector<3>{2, -3, 4}); - double result = P.evaluate(3); - REQUIRE(result == 29); - result = P(3); - REQUIRE(result == 29); + Polynomial<2> P(2, -3, 4); + REQUIRE(P(3) == 29); } + SECTION("primitive") { - Polynomial<2> P(TinyVector<3>{2, -3, 4}); + Polynomial<2> P(2, -3, 4); TinyVector<4> coefs = zero; Polynomial<3> Q(coefs); Q = primitive(P); - Polynomial<3> R(TinyVector<4>{0, 2, -3. / 2, 4. / 3}); + Polynomial<3> R(0, 2, -3. / 2, 4. / 3); REQUIRE(Q == R); } + SECTION("integrate") { - Polynomial<2> P(TinyVector<3>{2, -3, 3}); + Polynomial<2> P(2, -3, 3); double xinf = -1; double xsup = 1; double result = integrate(P, xinf, xsup); @@ -129,64 +137,65 @@ TEST_CASE("Polynomial", "[analysis]") result = symmetricIntegrate(P, 2); REQUIRE(result == 24); } + SECTION("derivative") { - Polynomial<2> P(TinyVector<3>{2, -3, 3}); + Polynomial<2> P(2, -3, 3); Polynomial<1> Q = derivative(P); - REQUIRE(Q == Polynomial<1>(TinyVector<2>{-3, 6})); + REQUIRE(Q == Polynomial<1>(-3, 6)); - Polynomial<0> P2(TinyVector<1>{3}); + Polynomial<0> P2(3); - Polynomial<0> R(TinyVector<1>{0}); + Polynomial<0> R(0); REQUIRE(derivative(P2) == R); } SECTION("affectation") { - Polynomial<2> Q(TinyVector<3>{2, -3, 3}); - Polynomial<4> R(TinyVector<5>{2, -3, 3, 0, 0}); - Polynomial<4> P(TinyVector<5>{0, 1, 2, 3, 3}); + Polynomial<2> Q(2, -3, 3); + Polynomial<4> R(2, -3, 3, 0, 0); + Polynomial<4> P(0, 1, 2, 3, 3); P = Q; REQUIRE(P == R); } SECTION("affectation addition") { - Polynomial<2> Q(TinyVector<3>{2, -3, 3}); - Polynomial<4> R(TinyVector<5>{2, -2, 5, 3, 3}); - Polynomial<4> P(TinyVector<5>{0, 1, 2, 3, 3}); + Polynomial<2> Q(2, -3, 3); + Polynomial<4> R(2, -2, 5, 3, 3); + Polynomial<4> P(0, 1, 2, 3, 3); P += Q; REQUIRE(P == R); } SECTION("power") { - Polynomial<2> P(TinyVector<3>{2, -3, 3}); - Polynomial<4> R(TinyVector<5>{4, -12, 21, -18, 9}); - Polynomial<1> Q(TinyVector<2>{0, 2}); + Polynomial<2> P(2, -3, 3); + Polynomial<4> R(4, -12, 21, -18, 9); + Polynomial<1> Q(0, 2); Polynomial<2> S = Q.pow<2>(2); REQUIRE(P.pow<2>(2) == R); - REQUIRE(S == Polynomial<2>(TinyVector<3>{0, 0, 4})); + REQUIRE(S == Polynomial<2>(0, 0, 4)); } SECTION("composition") { - Polynomial<2> P(TinyVector<3>{2, -3, 3}); - Polynomial<1> Q(TinyVector<2>{0, 2}); - Polynomial<2> R(TinyVector<3>{2, -1, 3}); - Polynomial<2> S(TinyVector<3>{1, 2, 2}); - REQUIRE(P.compose(Q) == Polynomial<2>(TinyVector<3>{2, -6, 12})); - REQUIRE(P.compose2(Q) == Polynomial<2>(TinyVector<3>{2, -6, 12})); - REQUIRE(R(S) == Polynomial<4>(TinyVector<5>{4, 10, 22, 24, 12})); + Polynomial<2> P(2, -3, 3); + Polynomial<1> Q(0, 2); + 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)); } SECTION("Lagrange polynomial") { - Polynomial<1> S(TinyVector<2>{0.5, -0.5}); + Polynomial<1> S(0.5, -0.5); Polynomial<1> Q; Q = lagrangePolynomial<1>(TinyVector<2>{-1, 1}, 0); REQUIRE(S == Q); - Polynomial<2> P(TinyVector<3>{0, -0.5, 0.5}); + Polynomial<2> P(0, -0.5, 0.5); Polynomial<2> R; R = lagrangePolynomial<2>(TinyVector<3>{-1, 0, 1}, 0); REQUIRE(R == P); diff --git a/tests/test_PolynomialBasis.cpp b/tests/test_PolynomialBasis.cpp index 3490498fdff4f108e518664254bf432c4df8e7d1..446a7278a36c5358b93069f5f0d5dc10ac00aa5a 100644 --- a/tests/test_PolynomialBasis.cpp +++ b/tests/test_PolynomialBasis.cpp @@ -34,12 +34,12 @@ TEST_CASE("PolynomialBasis", "[analysis]") PolynomialBasis<2> B; REQUIRE(B.displayType() == "undefined"); B.build(BasisType::canonical); - REQUIRE(B.elements()[1] == Polynomial<2>{TinyVector<3>{0, 1, 0}}); - REQUIRE(B.elements()[2] == Polynomial<2>{TinyVector<3>{0, 0, 1}}); + REQUIRE(B.elements()[1] == Polynomial<2>{0, 1, 0}); + REQUIRE(B.elements()[2] == Polynomial<2>{0, 0, 1}); PolynomialBasis<2> C; C.build(BasisType::taylor); REQUIRE(B.elements()[1] == C.elements()[1]); C.build(BasisType::taylor, 1); - REQUIRE(C.elements()[2] == Polynomial<2>{TinyVector<3>{1, -2, 1}}); + REQUIRE(C.elements()[2] == Polynomial<2>{1, -2, 1}); } }