math: add a simple polynomial template class.

src/Makefile.am

@@ -44,7 +44,7 @@ liblolcore_headers = \
     lol/math/functions.h lol/math/vector.h lol/math/half.h lol/math/real.h \
     lol/math/geometry.h lol/math/interp.h lol/math/rand.h lol/math/array2d.h \
     lol/math/array3d.h lol/math/constants.h lol/math/matrix.h lol/math/ops.h \
-    lol/math/transform.h \
+    lol/math/transform.h lol/math/polynomial.h lol/math/simplex_interpolator.h \
     \
     lol/algorithm/all.h \
     lol/algorithm/sort.h lol/algorithm/portal.h lol/algorithm/aabb_tree.h \

src/lol/math/all.h

@@ -22,3 +22,5 @@
 #include <lol/math/geometry.h>
 #include <lol/math/interp.h>
 #include <lol/math/rand.h>
+#include <lol/math/polynomial.h>
+#include <lol/math/simplex_interpolator.h>

src/lol/math/polynomial.h

@@ -0,0 +1,151 @@
+//
+// Lol Engine
+//
+// Copyright: (c) 2010-2014 Sam Hocevar <sam@hocevar.net>
+//   This program is free software; you can redistribute it and/or
+//   modify it under the terms of the Do What The Fuck You Want To
+//   Public License, Version 2, as published by Sam Hocevar. See
+//   http://www.wtfpl.net/ for more details.
+//
+
+#pragma once
+
+//
+// The polynomial class
+// --------------------
+//
+
+namespace lol
+{
+
+template<typename T>
+struct polynomial
+{
+    inline polynomial()
+    {
+    }
+
+    explicit polynomial(std::initializer_list<T> const &init)
+    {
+        for (auto a : init)
+            factors.Push(a);
+
+        reduce_degree();
+    }
+
+    /* We define the zero polynomial (with no coefficients) as having
+     * degree -1 on purpose. */
+    inline int degree() const
+    {
+        return (int)factors.Count() - 1;
+    }
+
+    T eval(T x) const
+    {
+        T ret = this->operator[](degree());
+        for (int i = degree() - 1; i >= 0; --i)
+            ret = ret * x + factors[i];
+        return ret;
+    }
+
+    inline T operator[](ptrdiff_t n) const
+    {
+        if (n < 0 || n > degree())
+            return T(0);
+
+        return factors[n];
+    }
+
+    polynomial<T> operator +() const
+    {
+        return *this;
+    }
+
+    polynomial<T> operator -() const
+    {
+        polynomial<T> ret;
+        for (auto a : factors)
+            ret.factors.Push(-a);
+        return ret;
+    }
+
+    polynomial<T> &operator +=(polynomial<T> const &p)
+    {
+        int min_degree = lol::min(p.degree(), degree());
+
+        for (int i = 0; i <= min_degree; ++i)
+            factors[i] += p[i];
+
+        for (int i = min_degree + 1; i <= p.degree(); ++i)
+            factors.Push(p[i]);
+
+        reduce_degree();
+        return *this;
+    }
+
+    polynomial<T> operator +(polynomial<T> const &p) const
+    {
+        return polynomial<T>(*this) += p;
+    }
+
+    polynomial<T> &operator -=(polynomial<T> const &p)
+    {
+        return *this += -p;
+    }
+
+    polynomial<T> operator -(polynomial<T> const &p) const
+    {
+        return polynomial<T>(*this) += -p;
+    }
+
+    polynomial<T> &operator *=(T const &k)
+    {
+        for (auto &a : factors)
+            a *= k;
+        reduce_degree();
+        return *this;
+    }
+
+    polynomial<T> operator *(T const &k) const
+    {
+        return polynomial<T>(*this) *= k;
+    }
+
+    polynomial<T> &operator *=(polynomial<T> const &p)
+    {
+        return *this = *this * p;
+    }
+
+    polynomial<T> operator *(polynomial<T> const &p) const
+    {
+        polynomial<T> ret;
+        polynomial<T> const &q = *this;
+
+        if (p.degree() >= 0 && q.degree() >= 0)
+        {
+            int n = p.degree() + q.degree();
+            for (int i = 0; i <= n; ++i)
+                ret.factors.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.reduce_degree();
+        }
+
+        return ret;
+    }
+
+private:
+    void reduce_degree()
+    {
+        while (factors.Count() && factors.Last() == T(0))
+            factors.Pop();
+    }
+
+    array<T> factors;
+};
+
+} /* namespace lol */
+

src/lolcore.vcxproj

@@ -333,9 +333,11 @@
     <ClInclude Include="lol\math\interp.h" />
     <ClInclude Include="lol\math\matrix.h" />
     <ClInclude Include="lol\math\ops.h" />
+    <ClInclude Include="lol\math\polynomial.h" />
     <ClInclude Include="lol\math\rand.h" />
     <ClInclude Include="lol\math\real.h" />
     <ClInclude Include="lol\math\remez.h" />
+    <ClInclude Include="lol\math\simplex_interpolator.h" />
     <ClInclude Include="lol\math\transform.h" />
     <ClInclude Include="lol\math\vector.h" />
     <ClInclude Include="lol\public.h" />

src/lolcore.vcxproj.filters

@@ -453,6 +453,9 @@
     <ClInclude Include="lol\math\ops.h">
       <Filter>lol\math</Filter>
     </ClInclude>
+    <ClInclude Include="lol\math\polynomial.h">
+      <Filter>lol\math</Filter>
+    </ClInclude>
     <ClInclude Include="lol\math\rand.h">
       <Filter>lol\math</Filter>
     </ClInclude>
@@ -462,6 +465,9 @@
     <ClInclude Include="lol\math\remez.h">
       <Filter>lol\math</Filter>
     </ClInclude>
+    <ClInclude Include="lol\math\simplex_interpolator.h">
+      <Filter>lol\math</Filter>
+    </ClInclude>
     <ClInclude Include="lol\math\transform.h">
       <Filter>lol\math</Filter>
     </ClInclude>

src/t/Makefile.am

@@ -20,7 +20,8 @@ test_math_SOURCES = test-common.cpp \
     math/array2d.cpp math/array3d.cpp math/arraynd.cpp math/box.cpp \
     math/cmplx.cpp math/half.cpp math/interp.cpp math/matrix.cpp \
     math/quat.cpp math/rand.cpp math/real.cpp math/rotation.cpp \
-    math/trig.cpp math/vector.cpp math/simplex_interpolator.cpp
+    math/trig.cpp math/vector.cpp math/simplex_interpolator.cpp \
+    math/polynomial.cpp
 test_math_CPPFLAGS = $(AM_CPPFLAGS) -I$(top_srcdir)/tools/lolunit
 test_math_DEPENDENCIES = @LOL_DEPS@
 

src/t/math/polynomial.cpp

@@ -0,0 +1,165 @@
+//
+// Lol Engine
+//
+// Copyright: (c) 2010-2014 Sam Hocevar <sam@hocevar.net>
+//   This program is free software; you can redistribute it and/or
+//   modify it under the terms of the Do What The Fuck You Want To
+//   Public License, Version 2, as published by Sam Hocevar. See
+//   http://www.wtfpl.net/ for more details.
+//
+
+#include <lol/engine-internal.h>
+
+#include <lolunit.h>
+
+namespace lol
+{
+
+lolunit_declare_fixture(PolynomialTest)
+{
+    lolunit_declare_test(Declaration)
+    {
+        polynomial<float> p;
+        polynomial<real> q;
+    }
+
+    lolunit_declare_test(Init)
+    {
+        polynomial<float> p { };
+        lolunit_assert_equal(p[0], 0.f);
+        lolunit_assert_equal(p[1], 0.f);
+        lolunit_assert_equal(p[2], 0.f);
+        lolunit_assert_equal(p.degree(), -1);
+
+        polynomial<float> q { 1.f };
+        lolunit_assert_equal(q[0], 1.f);
+        lolunit_assert_equal(q[1], 0.f);
+        lolunit_assert_equal(q[2], 0.f);
+        lolunit_assert_equal(q.degree(), 0);
+
+        polynomial<float> r { 1.f, 2.f };
+        lolunit_assert_equal(r[0], 1.f);
+        lolunit_assert_equal(r[1], 2.f);
+        lolunit_assert_equal(r[2], 0.f);
+        lolunit_assert_equal(r.degree(), 1);
+
+        polynomial<float> s { 0.f };
+        lolunit_assert_equal(s[0], 0.f);
+        lolunit_assert_equal(s[1], 0.f);
+        lolunit_assert_equal(s[2], 0.f);
+        lolunit_assert_equal(s.degree(), -1);
+    }
+
+    lolunit_declare_test(Eval)
+    {
+        /* Special null polynomial */
+        polynomial<float> p;
+
+        float a = p.eval(42.f);
+        lolunit_assert_equal(a, 0.f);
+    }
+
+    lolunit_declare_test(Eval0)
+    {
+        /* Constant polynomial p(x) = 1 */
+        polynomial<float> p { 1.f };
+
+        float a = p.eval(42.f);
+        lolunit_assert_equal(a, 1.f);
+    }
+
+    lolunit_declare_test(Eval1)
+    {
+        /* p(x) = 1 + 2x */
+        polynomial<float> p { 1.f, 2.f };
+
+        float a = p.eval(0.f);
+        lolunit_assert_equal(a, 1.f);
+
+        float b = p.eval(1.f);
+        lolunit_assert_equal(b, 3.f);
+
+        float c = p.eval(2.f);
+        lolunit_assert_equal(c, 5.f);
+    }
+
+    lolunit_declare_test(Eval2)
+    {
+        /* p(x) = 1 + 2x + 3x² */
+        polynomial<float> p { 1.f, 2.f, 3.f };
+
+        float a = p.eval(0.f);
+        lolunit_assert_equal(a, 1.f);
+
+        float b = p.eval(1.f);
+        lolunit_assert_equal(b, 6.f);
+
+        float c = p.eval(2.f);
+        lolunit_assert_equal(c, 17.f);
+    }
+
+    lolunit_declare_test(UnaryPlusMinus)
+    {
+        /* p(x) = 1 + 2x + 3x² */
+        polynomial<float> p { 1.f, 2.f, 3.f };
+        polynomial<float> q = +p;
+        polynomial<float> r = -p;
+
+        lolunit_assert_equal(q[0], 1.f);
+        lolunit_assert_equal(q[1], 2.f);
+        lolunit_assert_equal(q[2], 3.f);
+
+        lolunit_assert_equal(r[0], -1.f);
+        lolunit_assert_equal(r[1], -2.f);
+        lolunit_assert_equal(r[2], -3.f);
+    }
+
+    lolunit_declare_test(Addition)
+    {
+        /* p(x) = 1 + 2x + 3x² */
+        /* q(x) = 4 + 5x */
+        polynomial<float> p { 1.f, 2.f, 3.f };
+        polynomial<float> q { 4.f, 5.f };
+
+        /* r(x) = 5 + 7x + 3x² */
+        polynomial<float> r = p + q;
+        lolunit_assert_equal(r.degree(), 2);
+        lolunit_assert_equal(r[0], 5.f);
+        lolunit_assert_equal(r[1], 7.f);
+        lolunit_assert_equal(r[2], 3.f);
+    }
+
+    lolunit_declare_test(Subtraction)
+    {
+        /* p(x) = 1 + 2x + 3x² */
+        /* q(x) = 4 + 5x */
+        polynomial<float> p { 1.f, 2.f, 3.f };
+        polynomial<float> q { 4.f, 5.f };
+
+        /* r(x) = -3 + -3x + 3x² */
+        polynomial<float> r = p - q;
+        lolunit_assert_equal(r.degree(), 2);
+        lolunit_assert_equal(r[0], -3.f);
+        lolunit_assert_equal(r[1], -3.f);
+        lolunit_assert_equal(r[2], 3.f);
+    }
+
+    lolunit_declare_test(Multiplication)
+    {
+        /* p(x) = 1 + 2x + 3x² */
+        /* q(x) = 4 + 5x */
+        polynomial<float> p { 1.f, 2.f, 3.f };
+        polynomial<float> q { 4.f, 5.f };
+
+        /* r(x) = 4 + 13x + 22x² + 15x³ */
+        polynomial<float> r = p * q;
+        lolunit_assert_equal(r.degree(), 3);
+        lolunit_assert_equal(r[0], 4.f);
+        lolunit_assert_equal(r[1], 13.f);
+        lolunit_assert_equal(r[2], 22.f);
+        lolunit_assert_equal(r[3], 15.f);
+    }
+};
+
+} /* namespace lol */
+

src/t/test-math.vcxproj

@@ -52,10 +52,12 @@
     <ClCompile Include="math\half.cpp" />
     <ClCompile Include="math\interp.cpp" />
     <ClCompile Include="math\matrix.cpp" />
+    <ClCompile Include="math\polynomial.cpp" />
     <ClCompile Include="math\quat.cpp" />
     <ClCompile Include="math\rand.cpp" />
     <ClCompile Include="math\real.cpp" />
     <ClCompile Include="math\rotation.cpp" />
+    <ClCompile Include="math\simplex_interpolation.cpp" />
     <ClCompile Include="math\trig.cpp" />
     <ClCompile Include="math\vector.cpp" />
   </ItemGroup>