math: add a factory for Chebyshev polynomials.

src/lol/math/polynomial.h

@@ -15,36 +15,72 @@
 // --------------------
 //
 
+#include <functional>
+
 namespace lol
 {
 
 template<typename T>
 struct polynomial
 {
-    inline polynomial()
-    {
-    }
+    /* The zero polynomial */
+    explicit inline polynomial() {}
 
+    /* Create a polynomial from a list of coefficients */
     explicit polynomial(std::initializer_list<T> const &init)
     {
         for (auto a : init)
-            factors.Push(a);
+            m_coefficients.Push(a);
 
         reduce_degree();
     }
 
+    /* Factory for Chebyshev polynomials */
+    static polynomial<T> chebyshev(int n)
+    {
+        /* Use T0(x) = 1, T1(x) = x, Tn(x) = 2 x Tn-1(x) - Tn-2(x) */
+        std::function<int64_t (int, int)> 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<T> 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<T> operator +() const
@@ -64,8 +100,8 @@ struct polynomial
     polynomial<T> operator -() const
     {
         polynomial<T> 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<T> &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<T> factors;
+    /* The polynomial coefficients */
+    array<T> m_coefficients;
 };
 
+template<typename T>
+polynomial<T> operator *(T const &k, polynomial<T> const &p)
+{
+    /* We assume k * p is the same as p * k */
+    return p * k;
+}
+
 } /* namespace lol */
 

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<float> t0 = polynomial<float>::chebyshev(0);
+        polynomial<float> t1 = polynomial<float>::chebyshev(1);
+        polynomial<float> t2 = polynomial<float>::chebyshev(2);
+        polynomial<float> t3 = polynomial<float>::chebyshev(3);
+        polynomial<float> t4 = polynomial<float>::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 */