From 0668d0d5a6f2db28a68c1ebdbf124365b6416db4 Mon Sep 17 00:00:00 2001 From: Sam Hocevar Date: Mon, 8 Dec 2014 08:05:27 +0000 Subject: [PATCH] math: add a factory for Chebyshev polynomials. --- src/lol/math/polynomial.h | 78 ++++++++++++++++++++++++++++++--------- src/t/math/polynomial.cpp | 35 ++++++++++++++++++ 2 files changed, 96 insertions(+), 17 deletions(-) diff --git a/src/lol/math/polynomial.h b/src/lol/math/polynomial.h index bba56b74..48f182a4 100644 --- a/src/lol/math/polynomial.h +++ b/src/lol/math/polynomial.h @@ -15,36 +15,72 @@ // -------------------- // +#include + namespace lol { template struct polynomial { - inline polynomial() - { - } + /* The zero polynomial */ + explicit inline polynomial() {} + /* Create a polynomial from a list of coefficients */ explicit polynomial(std::initializer_list const &init) { for (auto a : init) - factors.Push(a); + m_coefficients.Push(a); reduce_degree(); } + /* Factory for Chebyshev polynomials */ + static polynomial chebyshev(int n) + { + /* Use T0(x) = 1, T1(x) = x, Tn(x) = 2 x Tn-1(x) - Tn-2(x) */ + std::function coeff = [&](int i, int j) + { + if (i > j || i < 0 || ((j ^ i) & 1)) + return (int64_t)0; + if (j < 2) + return (int64_t)1; + return 2 * coeff(i - 1, j - 1) - coeff(i, j - 2); + }; + + polynomial ret; + for (int k = 0; k <= n; ++k) + ret.m_coefficients.Push(T(coeff(k, n))); + return ret; + } + /* We define the zero polynomial (with no coefficients) as having * degree -1 on purpose. */ inline int degree() const { - return (int)factors.Count() - 1; + return (int)m_coefficients.Count() - 1; + } + + /* Set one of the polynomial’s coefficients */ + void set(ptrdiff_t n, T const &a) + { + ASSERT(n >= 0); + + if (n > degree() && !a) + return; + + while (n > degree()) + m_coefficients.Push(T(0)); + + m_coefficients[n] = a; + reduce_degree(); } T eval(T x) const { T ret = this->operator[](degree()); for (int i = degree() - 1; i >= 0; --i) - ret = ret * x + factors[i]; + ret = ret * x + m_coefficients[i]; return ret; } @@ -53,7 +89,7 @@ struct polynomial if (n < 0 || n > degree()) return T(0); - return factors[n]; + return m_coefficients[n]; } polynomial operator +() const @@ -64,8 +100,8 @@ struct polynomial polynomial operator -() const { polynomial ret; - for (auto a : factors) - ret.factors.Push(-a); + for (auto a : m_coefficients) + ret.m_coefficients.Push(-a); return ret; } @@ -74,10 +110,10 @@ struct polynomial int min_degree = lol::min(p.degree(), degree()); for (int i = 0; i <= min_degree; ++i) - factors[i] += p[i]; + m_coefficients[i] += p[i]; for (int i = min_degree + 1; i <= p.degree(); ++i) - factors.Push(p[i]); + m_coefficients.Push(p[i]); reduce_degree(); return *this; @@ -100,7 +136,7 @@ struct polynomial polynomial &operator *=(T const &k) { - for (auto &a : factors) + for (auto &a : m_coefficients) a *= k; reduce_degree(); return *this; @@ -125,11 +161,11 @@ struct polynomial { int n = p.degree() + q.degree(); for (int i = 0; i <= n; ++i) - ret.factors.Push(T(0)); + ret.m_coefficients.Push(T(0)); for (int i = 0; i <= p.degree(); ++i) for (int j = 0; j <= q.degree(); ++j) - ret.factors[i + j] += p[i] * q[j]; + ret.m_coefficients[i + j] += p[i] * q[j]; ret.reduce_degree(); } @@ -140,12 +176,20 @@ struct polynomial private: void reduce_degree() { - while (factors.Count() && factors.Last() == T(0)) - factors.Pop(); + while (m_coefficients.Count() && m_coefficients.Last() == T(0)) + m_coefficients.Pop(); } - array factors; + /* The polynomial coefficients */ + array m_coefficients; }; +template +polynomial operator *(T const &k, polynomial const &p) +{ + /* We assume k * p is the same as p * k */ + return p * k; +} + } /* namespace lol */ diff --git a/src/t/math/polynomial.cpp b/src/t/math/polynomial.cpp index ca69b36b..72de936b 100644 --- a/src/t/math/polynomial.cpp +++ b/src/t/math/polynomial.cpp @@ -159,6 +159,41 @@ lolunit_declare_fixture(PolynomialTest) lolunit_assert_equal(r[2], 22.f); lolunit_assert_equal(r[3], 15.f); } + + lolunit_declare_test(Chebyshev) + { + polynomial t0 = polynomial::chebyshev(0); + polynomial t1 = polynomial::chebyshev(1); + polynomial t2 = polynomial::chebyshev(2); + polynomial t3 = polynomial::chebyshev(3); + polynomial t4 = polynomial::chebyshev(4); + + /* Taken from the sequence at http://oeis.org/A028297 */ + lolunit_assert_equal(t0.degree(), 0); + lolunit_assert_equal(t0[0], 1.f); + + lolunit_assert_equal(t1.degree(), 1); + lolunit_assert_equal(t1[0], 0.f); + lolunit_assert_equal(t1[1], 1.f); + + lolunit_assert_equal(t2.degree(), 2); + lolunit_assert_equal(t2[0], -1.f); + lolunit_assert_equal(t2[1], 0.f); + lolunit_assert_equal(t2[2], 2.f); + + lolunit_assert_equal(t3.degree(), 3); + lolunit_assert_equal(t3[0], 0.f); + lolunit_assert_equal(t3[1], -3.f); + lolunit_assert_equal(t3[2], 0.f); + lolunit_assert_equal(t3[3], 4.f); + + lolunit_assert_equal(t4.degree(), 4); + lolunit_assert_equal(t4[0], 1.f); + lolunit_assert_equal(t4[1], 0.f); + lolunit_assert_equal(t4[2], -8.f); + lolunit_assert_equal(t4[3], 0.f); + lolunit_assert_equal(t4[4], 8.f); + } }; } /* namespace lol */