diff --git a/src/analysis/Polynomial.hpp b/src/analysis/Polynomial.hpp
index 3d534dbff3a140872a89e2c73ff5ff26c379de72..e43c4535221682d6b2e5c96b5492ad6d96f780c7 100644
--- a/src/analysis/Polynomial.hpp
+++ b/src/analysis/Polynomial.hpp
@@ -308,13 +308,13 @@ class Polynomial
const size_t Mr = P2.realDegree();
R.coefficients() = P1.coefficients();
Q.coefficients() = zero;
- for (ssize_t k = Nr - Mr; k >= 0; --k) {
- Q.coefficients()[k] = R.coefficients()[Mr + k] / P2.coefficients()[Mr];
- for (ssize_t j = Mr + k; j >= k; --j) {
- R.coefficients()[j] -= Q.coefficients()[k] * P2.coefficients()[j - k];
+ for (size_t k = Nr - Mr + 1; k > 0; --k) {
+ Q.coefficients()[k - 1] = R.coefficients()[Mr + k - 1] / P2.coefficients()[Mr];
+ for (size_t j = Mr + k - 1; j > (k - 1); --j) {
+ R.coefficients()[j] -= Q.coefficients()[k - 1] * P2.coefficients()[j - k];
}
}
- for (size_t j = Mr; j <= Nr; ++j) {
+ for (size_t j = Mr; j < Nr + 1; ++j) {
R.coefficients()[j] = 0;
}
}
@@ -402,7 +402,7 @@ class Polynomial
PUGS_INLINE
constexpr friend Polynomial<N>
- lagrangeToCanonical(const TinyVector<N + 1> lagrange_coefs, const TinyVector<N + 1, Polynomial<N>>& basis)
+ lagrangeToCanonical(const TinyVector<N + 1> lagrange_coefs, const std::array<Polynomial<N>, N + 1>& basis)
{
Polynomial<N> P(zero);
// lagrangeBasis({0, 0, 0}, basis);
@@ -426,12 +426,11 @@ template <size_t N>
PUGS_INLINE constexpr Polynomial<N> lagrangePolynomial(const TinyVector<N + 1> zeros, const size_t k);
template <size_t N>
-PUGS_INLINE constexpr TinyVector<N, Polynomial<N - 1>>
+PUGS_INLINE constexpr std::array<Polynomial<N - 1>, N>
lagrangeBasis(const TinyVector<N>& zeros)
{
static_assert(N >= 1, "invalid degree");
- TinyVector<N, Polynomial<N - 1>> basis;
- Polynomial<N - 1> lj;
+ std::array<Polynomial<N - 1>, N> basis;
for (size_t j = 0; j < N; ++j) {
basis[j] = lagrangePolynomial<N - 1>(zeros, j);
}
@@ -464,7 +463,7 @@ PUGS_INLINE constexpr auto
derivative(const Polynomial<N>& P)
{
if constexpr (N == 0) {
- return Polynomial<0>(0);
+ return Polynomial<0>(zero);
} else {
TinyVector<N> coefs;
for (size_t i = 0; i < N; ++i) {
@@ -488,7 +487,7 @@ lagrangePolynomial(const TinyVector<N + 1> zeros, const size_t k)
if (k == i)
continue;
double factor = 1. / (zeros[k] - zeros[i]);
- Polynomial<1> P({-zeros[i] * factor, factor});
+ Polynomial<1> P(TinyVector<2>{-zeros[i] * factor, factor});
lk *= P;
}
return lk;
diff --git a/src/analysis/PolynomialBasis.hpp b/src/analysis/PolynomialBasis.hpp
index 5bcd8af5ce7b5a062af6fd37ace63beadf8b418e..cf64978b11c21252dfa6b6e98382806e71703306 100644
--- a/src/analysis/PolynomialBasis.hpp
+++ b/src/analysis/PolynomialBasis.hpp
@@ -18,7 +18,7 @@ class PolynomialBasis
{
private:
static_assert((N >= 0), "Number of elements in the basis must be non-negative");
- TinyVector<N + 1, Polynomial<N>> m_elements;
+ std::array<Polynomial<N>, N + 1> m_elements;
BasisType m_basis_type;
PUGS_INLINE
constexpr PolynomialBasis<N>&
@@ -65,7 +65,7 @@ class PolynomialBasis
if (i == j)
continue;
double adim = 1. / (zeros[i] - zeros[j]);
- elements()[i] *= Polynomial<1>{{-zeros[j] * adim, adim}};
+ elements()[i] *= Polynomial<1>{TinyVector<2>{-zeros[j] * adim, adim}};
}
}
return *this;
@@ -113,14 +113,14 @@ class PolynomialBasis
}
PUGS_INLINE
- constexpr const TinyVector<N + 1, Polynomial<N>>&
+ constexpr const std::array<Polynomial<N>, N + 1>&
elements() const
{
return m_elements;
}
PUGS_INLINE
- constexpr TinyVector<N + 1, Polynomial<N>>&
+ constexpr std::array<Polynomial<N>, N + 1>&
elements()
{
return m_elements;
diff --git a/tests/test_Polynomial.cpp b/tests/test_Polynomial.cpp
index 74f8e295e0b14b5d77c141b3d7d00fed28b8dfc4..232ffcca02f445314d14b3418bc14924c5d19db8 100644
--- a/tests/test_Polynomial.cpp
+++ b/tests/test_Polynomial.cpp
@@ -22,75 +22,75 @@ TEST_CASE("Polynomial", "[analysis]")
{
SECTION("construction")
{
- REQUIRE_NOTHROW(Polynomial<2>({2, 3, 4}));
+ REQUIRE_NOTHROW(Polynomial<2>{TinyVector<3>{2, 3, 4}});
}
SECTION("degree")
{
- Polynomial<2> P({2, 3, 4});
+ Polynomial<2> P(TinyVector<3>{2, 3, 4});
REQUIRE(P.degree() == 2);
}
SECTION("equality")
{
- Polynomial<2> P({2, 3, 4});
- Polynomial<2> Q({2, 3, 4});
- Polynomial<2> S({2, 3, 5});
+ Polynomial<2> P(TinyVector<3>{2, 3, 4});
+ Polynomial<2> Q(TinyVector<3>{2, 3, 4});
+ Polynomial<2> S(TinyVector<3>{2, 3, 5});
REQUIRE(P == Q);
REQUIRE(P != S);
}
SECTION("addition")
{
- 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});
+ 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});
REQUIRE(S == (P + Q));
REQUIRE((T + P) == U);
}
SECTION("opposed")
{
- Polynomial<2> P({2, 3, 4});
+ Polynomial<2> P(TinyVector<3>{2, 3, 4});
Polynomial<2> Q = -P;
- REQUIRE(Q == Polynomial<2>({-2, -3, -4}));
+ REQUIRE(Q == Polynomial<2>(TinyVector<3>{-2, -3, -4}));
}
SECTION("difference")
{
- Polynomial<2> P({2, 3, 4});
- Polynomial<2> Q({3, 4, 5});
- Polynomial<2> D({-1, -1, -1});
+ Polynomial<2> P(TinyVector<3>{2, 3, 4});
+ Polynomial<2> Q(TinyVector<3>{3, 4, 5});
+ Polynomial<2> D(TinyVector<3>{-1, -1, -1});
REQUIRE(D == (P - Q));
- Polynomial<3> R({2, 3, 4, 1});
+ Polynomial<3> R(TinyVector<4>{2, 3, 4, 1});
REQUIRE(D == (P - Q));
- REQUIRE((P - R) == Polynomial<3>({0, 0, 0, -1}));
+ REQUIRE((P - R) == Polynomial<3>(TinyVector<4>{0, 0, 0, -1}));
R -= P;
- REQUIRE(R == Polynomial<3>({0, 0, 0, 1}));
+ REQUIRE(R == Polynomial<3>(TinyVector<4>{0, 0, 0, 1}));
}
SECTION("product_by_scalar")
{
- Polynomial<2> P({2, 3, 4});
- Polynomial<2> M({6, 9, 12});
+ Polynomial<2> P(TinyVector<3>{2, 3, 4});
+ Polynomial<2> M(TinyVector<3>{6, 9, 12});
REQUIRE(M == (P * 3));
REQUIRE(M == (3 * P));
}
SECTION("product")
{
- Polynomial<2> P({2, 3, 4});
- Polynomial<3> Q({1, 2, -1, 1});
+ Polynomial<2> P(TinyVector<3>{2, 3, 4});
+ Polynomial<3> Q(TinyVector<4>{1, 2, -1, 1});
Polynomial<4> R;
Polynomial<5> S;
R = P;
S = P;
S *= Q;
- REQUIRE(Polynomial<5>({2, 7, 8, 7, -1, 4}) == (P * Q));
- REQUIRE(Polynomial<5>({2, 7, 8, 7, -1, 4}) == S);
+ 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_THROWS_AS(R *= Q, AssertError);
}
SECTION("divide")
{
- Polynomial<2> P({1, 0, 1});
- Polynomial<1> Q({0, 1});
- Polynomial<1> Q1({0, 1});
+ Polynomial<2> P(TinyVector<3>{1, 0, 1});
+ Polynomial<1> Q(TinyVector<2>{0, 1});
+ Polynomial<1> Q1(TinyVector<2>{0, 1});
Polynomial<2> R;
Polynomial<2> S;
@@ -99,26 +99,12 @@ TEST_CASE("Polynomial", "[analysis]")
REQUIRE(Q1.realDegree() == 1);
divide(P, Q1, R, S);
- REQUIRE(Polynomial<2>({1, 0, 0}) == S);
- REQUIRE(Polynomial<2>({0, 1, 0}) == R);
- }
- SECTION("divide2")
- {
- Polynomial<7> P({1, -3, 0.6, -2.6, 7, 1.8, 4, -2});
- Polynomial<3> Q({-1, 1, 0, 2});
-
- Polynomial<7> R;
- Polynomial<7> S;
- REQUIRE(P.realDegree() == 7);
- REQUIRE(Q.realDegree() == 3);
-
- divide(P, Q, R, S);
- REQUIRE(Polynomial<7>({-1, 2, 1.4, 2, -1, 0, 0, 0}) == R);
- REQUIRE(Polynomial<7>({0, 0, 0, 0, 0, 0, 0, 0}) == S);
+ REQUIRE(Polynomial<2>(TinyVector<3>{1, 0, 0}) == S);
+ REQUIRE(Polynomial<2>(TinyVector<3>{0, 1, 0}) == R);
}
SECTION("evaluation")
{
- Polynomial<2> P({2, -3, 4});
+ Polynomial<2> P(TinyVector<3>{2, -3, 4});
double result = P.evaluate(3);
REQUIRE(result == 29);
result = P(3);
@@ -126,16 +112,16 @@ TEST_CASE("Polynomial", "[analysis]")
}
SECTION("primitive")
{
- Polynomial<2> P({2, -3, 4});
+ Polynomial<2> P(TinyVector<3>{2, -3, 4});
TinyVector<4> coefs = zero;
Polynomial<3> Q(coefs);
Q = primitive(P);
- Polynomial<3> R({0, 2, -3. / 2, 4. / 3});
+ Polynomial<3> R(TinyVector<4>{0, 2, -3. / 2, 4. / 3});
REQUIRE(Q == R);
}
SECTION("integrate")
{
- Polynomial<2> P({2, -3, 3});
+ Polynomial<2> P(TinyVector<3>{2, -3, 3});
double xinf = -1;
double xsup = 1;
double result = integrate(P, xinf, xsup);
@@ -145,66 +131,66 @@ TEST_CASE("Polynomial", "[analysis]")
}
SECTION("derivative")
{
- Polynomial<2> P({2, -3, 3});
+ Polynomial<2> P(TinyVector<3>{2, -3, 3});
Polynomial<1> Q = derivative(P);
- REQUIRE(Q == Polynomial<1>({-3, 6}));
+ REQUIRE(Q == Polynomial<1>(TinyVector<2>{-3, 6}));
- Polynomial<0> P2(3);
+ Polynomial<0> P2(TinyVector<1>{3});
- Polynomial<0> R(0);
+ Polynomial<0> R(TinyVector<1>{0});
REQUIRE(derivative(P2) == R);
}
SECTION("affectation")
{
- Polynomial<2> Q({2, -3, 3});
- Polynomial<4> R({2, -3, 3, 0, 0});
- Polynomial<4> P({0, 1, 2, 3, 3});
+ 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});
P = Q;
REQUIRE(P == R);
}
SECTION("affectation addition")
{
- Polynomial<2> Q({2, -3, 3});
- Polynomial<4> R({2, -2, 5, 3, 3});
- Polynomial<4> P({0, 1, 2, 3, 3});
+ 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});
P += Q;
REQUIRE(P == R);
}
SECTION("power")
{
- Polynomial<2> P({2, -3, 3});
- Polynomial<4> R({4, -12, 21, -18, 9});
- Polynomial<1> Q({0, 2});
+ 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> S = Q.pow<2>(2);
REQUIRE(P.pow<2>(2) == R);
- REQUIRE(S == Polynomial<2>({0, 0, 4}));
+ REQUIRE(S == Polynomial<2>(TinyVector<3>{0, 0, 4}));
}
SECTION("composition")
{
- 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}));
+ 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}));
}
SECTION("Lagrange polynomial")
{
- Polynomial<1> S({0.5, -0.5});
+ Polynomial<1> S(TinyVector<2>{0.5, -0.5});
Polynomial<1> Q;
- Q = lagrangePolynomial<1>({-1, 1}, 0);
+ Q = lagrangePolynomial<1>(TinyVector<2>{-1, 1}, 0);
REQUIRE(S == Q);
- Polynomial<2> P({0, -0.5, 0.5});
+ Polynomial<2> P(TinyVector<3>{0, -0.5, 0.5});
Polynomial<2> R;
- R = lagrangePolynomial<2>({-1, 0, 1}, 0);
+ R = lagrangePolynomial<2>(TinyVector<3>{-1, 0, 1}, 0);
REQUIRE(R == P);
- const TinyVector<3, Polynomial<2>> basis = lagrangeBasis(TinyVector<3>{-1, 0, 1});
- REQUIRE(lagrangeToCanonical({1, 0, 1}, basis) == Polynomial<2>({0, 0, 1}));
+ const std::array<Polynomial<2>, 3> basis = lagrangeBasis(TinyVector<3>{-1, 0, 1});
+ REQUIRE(lagrangeToCanonical(TinyVector<3>{1, 0, 1}, basis) == Polynomial<2>(TinyVector<3>{0, 0, 1}));
}
}
diff --git a/tests/test_PolynomialBasis.cpp b/tests/test_PolynomialBasis.cpp
index b2c072e331e2e9bef939bf058bcc90550d5660f7..3490498fdff4f108e518664254bf432c4df8e7d1 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>{{0, 1, 0}});
- REQUIRE(B.elements()[2] == Polynomial<2>{{0, 0, 1}});
+ REQUIRE(B.elements()[1] == Polynomial<2>{TinyVector<3>{0, 1, 0}});
+ REQUIRE(B.elements()[2] == Polynomial<2>{TinyVector<3>{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>{{1, -2, 1}});
+ REQUIRE(C.elements()[2] == Polynomial<2>{TinyVector<3>{1, -2, 1}});
}
}
diff --git a/tests/test_TinyMatrix.cpp b/tests/test_TinyMatrix.cpp
index 4d5667dbff133ca432e115c8af2dad1b7665caf5..39bb8a9c1353b8336d60db9b4f77afb4a04c4710 100644
--- a/tests/test_TinyMatrix.cpp
+++ b/tests/test_TinyMatrix.cpp
@@ -225,8 +225,8 @@ TEST_CASE("TinyMatrix", "[algebra]")
REQUIRE(I(2, 2) == Catch::Approx(1).epsilon(1E-14));
}
{
- const TinyMatrix<4, double> A4(2, 3, 1, 4, -1, 5, -2, 3, 4, 1, 2, 3, 4, -1, -1, 0);
- const TinyMatrix<4, double> I = inverse(A4) * A4;
+ const TinyMatrix<4> A4(2, 3, 1, 4, -1, 5, -2, 3, 4, 1, 2, 3, 4, -1, -1, 0);
+ const TinyMatrix<4> I = inverse(A4) * A4;
REQUIRE(I(0, 0) == Catch::Approx(1).epsilon(1E-14));
REQUIRE(I(0, 1) == Catch::Approx(0).margin(1E-14));