math: move the Remez algorithm implementation to the core.

src/Makefile.am

@@ -2,7 +2,7 @@
 noinst_LIBRARIES = liblol.a
 
 liblol_a_SOURCES = \
-    core.h matrix.cpp matrix.h tiler.cpp tiler.h dict.cpp dict.h \
+    core.h matrix.cpp real.cpp tiler.cpp tiler.h dict.cpp dict.h \
     audio.cpp audio.h scene.cpp scene.h font.cpp font.h layer.cpp layer.h \
     map.cpp map.h entity.cpp entity.h ticker.cpp ticker.h lolgl.h \
     tileset.cpp tileset.h forge.cpp forge.h video.cpp video.h log.cpp log.h \
@@ -11,9 +11,9 @@ liblol_a_SOURCES = \
     text.cpp text.h emitter.cpp emitter.h numeric.h hash.cpp hash.h \
     worldentity.cpp worldentity.h gradient.cpp gradient.h half.cpp half.h \
     platform.cpp platform.h sprite.cpp sprite.h trig.cpp trig.h \
-    real.cpp real.h \
     \
     lol/unit.h \
+    lol/math/matrix.h lol/math/real.h lol/math/remez.h \
     \
     application/application.cpp application/application.h \
     eglapp.cpp eglapp.h \

src/core.h

@@ -64,8 +64,8 @@ static inline int isnan(float f)
 // Base types
 #include "trig.h"
 #include "half.h"
-#include "real.h"
-#include "matrix.h"
+#include "lol/math/real.h"
+#include "lol/math/matrix.h"
 #include "numeric.h"
 #include "timer.h"
 #include "thread/thread.h"

src/eglapp.h

@@ -16,7 +16,7 @@
 #if !defined __LOL_EGLAPP_H__
 #define __LOL_EGLAPP_H__
 
-#include "matrix.h"
+#include "lol/math/matrix.h"
 
 namespace lol
 {

src/image/image.h

@@ -16,7 +16,7 @@
 #if !defined __LOL_IMAGE_H__
 #define __LOL_IMAGE_H__
 
-#include "matrix.h"
+#include "lol/math/matrix.h"
 
 namespace lol
 {

src/input.h

@@ -16,7 +16,7 @@
 #if !defined __LOL_INPUT_H__
 #define __LOL_INPUT_H__
 
-#include "matrix.h"
+#include "lol/math/matrix.h"
 
 namespace lol
 {

src/matrix.h → src/lol/math/matrix.h

@@ -13,8 +13,8 @@
 // ------------------
 //
 
-#if !defined __LOL_MATRIX_H__
-#define __LOL_MATRIX_H__
+#if !defined __LOL_MATH_MATRIX_H__
+#define __LOL_MATH_MATRIX_H__
 
 #include <cmath>
 #if !defined __ANDROID__
@@ -24,6 +24,9 @@
 namespace lol
 {
 
+class half;
+class real;
+
 #define VECTOR_TYPES(tname, suffix) \
     template <typename T> struct tname; \
     typedef tname<half> f16##suffix; \
@@ -587,7 +590,82 @@ template <typename T> struct Mat4
     Vec4<T> v[4];
 };
 
+/*
+ * Arbitrarily-sized square matrices; for now this only supports
+ * naive inversion and is used for the Remez inversion method.
+ */
+
+template<int N, typename T> struct Mat
+{
+    inline Mat<N, T>() {}
+
+    Mat(T x)
+    {
+        for (int j = 0; j < N; j++)
+            for (int i = 0; i < N; i++)
+                if (i == j)
+                    m[i][j] = x;
+                else
+                    m[i][j] = 0;
+    }
+
+    /* Naive matrix inversion */
+    Mat<N, T> inv() const
+    {
+        Mat a = *this, b((T)1);
+
+        /* Inversion method: iterate through all columns and make sure
+         * all the terms are 1 on the diagonal and 0 everywhere else */
+        for (int i = 0; i < N; i++)
+        {
+            /* If the expected coefficient is zero, add one of
+             * the other lines. The first we meet will do. */
+            if (!a.m[i][i])
+            {
+                for (int j = i + 1; j < N; j++)
+                {
+                    if (!a.m[i][j])
+                        continue;
+                    /* Add row j to row i */
+                    for (int n = 0; n < N; n++)
+                    {
+                        a.m[n][i] += a.m[n][j];
+                        b.m[n][i] += b.m[n][j];
+                    }
+                    break;
+                }
+            }
+
+            /* Now we know the diagonal term is non-zero. Get its inverse
+             * and use that to nullify all other terms in the column */
+            T x = (T)1 / a.m[i][i];
+            for (int j = 0; j < N; j++)
+            {
+                if (j == i)
+                    continue;
+                T mul = x * a.m[i][j];
+                for (int n = 0; n < N; n++)
+                {
+                    a.m[n][j] -= mul * a.m[n][i];
+                    b.m[n][j] -= mul * b.m[n][i];
+                }
+            }
+
+            /* Finally, ensure the diagonal term is 1 */
+            for (int n = 0; n < N; n++)
+            {
+                a.m[n][i] *= x;
+                b.m[n][i] *= x;
+            }
+        }
+
+        return b;
+    }
+
+    T m[N][N];
+};
+
 } /* namespace lol */
 
-#endif // __LOL_MATRIX_H__
+#endif // __LOL_MATH_MATRIX_H__
 

src/real.h → src/lol/math/real.h

@@ -13,8 +13,8 @@
 // --------------
 //
 
-#if !defined __LOL_REAL_H__
-#define __LOL_REAL_H__
+#if !defined __LOL_MATH_REAL_H__
+#define __LOL_MATH_REAL_H__
 
 #include <stdint.h>
 
@@ -147,5 +147,5 @@ private:
 
 } /* namespace lol */
 
-#endif // __LOL_REAL_H__
+#endif // __LOL_MATH_REAL_H__
 

test/math/remez-solver.h → src/lol/math/remez.h

@@ -1,26 +1,34 @@
 //
-// Lol Engine - Sample math program: Chebyshev polynomials
+// Lol Engine
 //
-// Copyright: (c) 2005-2011 Sam Hocevar <sam@hocevar.net>
+// Copyright: (c) 2010-2011 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://sam.zoy.org/projects/COPYING.WTFPL for more details.
 //
 
-#if !defined __REMEZ_SOLVER_H__
-#define __REMEZ_SOLVER_H__
+//
+// The Remez class
+// ---------------
+//
+
+#if !defined __LOL_MATH_REMEZ_H__
+#define __LOL_MATH_REMEZ_H__
 
-template<int ORDER> class RemezSolver
+namespace lol
+{
+
+template<int ORDER, typename T> class RemezSolver
 {
 public:
-    typedef real RealFunc(real const &x);
+    typedef T RealFunc(T const &x);
 
     RemezSolver()
     {
     }
 
-    void Run(real a, real b, RealFunc *func, RealFunc *weight, int steps)
+    void Run(T a, T b, RealFunc *func, RealFunc *weight, int steps)
     {
         m_func = func;
         m_weight = weight;
@@ -49,15 +57,15 @@ public:
         PrintPoly();
     }
 
-    real ChebyEval(real const &x)
+    T ChebyEval(T const &x)
     {
-        real ret = 0.0, xn = 1.0;
+        T ret = 0.0, xn = 1.0;
 
         for (int i = 0; i < ORDER + 1; i++)
         {
-            real mul = 0;
+            T mul = 0;
             for (int j = 0; j < ORDER + 1; j++)
-                mul += coeff[j] * (real)Cheby(j, i);
+                mul += coeff[j] * (T)Cheby(j, i);
             ret += mul * xn;
             xn *= x;
         }
@@ -68,20 +76,20 @@ public:
     void Init()
     {
         /* Pick up x_i where error will be 0 and compute f(x_i) */
-        real fxn[ORDER + 1];
+        T fxn[ORDER + 1];
         for (int i = 0; i < ORDER + 1; i++)
         {
-            zeroes[i] = (real)(2 * i - ORDER) / (real)(ORDER + 1);
+            zeroes[i] = (T)(2 * i - ORDER) / (T)(ORDER + 1);
             fxn[i] = Value(zeroes[i]);
         }
 
         /* We build a matrix of Chebishev evaluations: row i contains the
          * evaluations of x_i for polynomial order n = 0, 1, ... */
-        Matrix<ORDER + 1> mat;
+        lol::Mat<ORDER + 1, T> mat;
         for (int i = 0; i < ORDER + 1; i++)
         {
             /* Compute the powers of x_i */
-            real powers[ORDER + 1];
+            T powers[ORDER + 1];
             powers[0] = 1.0;
             for (int n = 1; n < ORDER + 1; n++)
                  powers[n] = powers[n - 1] * zeroes[i];
@@ -89,9 +97,9 @@ public:
             /* Compute the Chebishev evaluations at x_i */
             for (int n = 0; n < ORDER + 1; n++)
             {
-                real sum = 0.0;
+                T sum = 0.0;
                 for (int k = 0; k < ORDER + 1; k++)
-                    sum += (real)Cheby(n, k) * powers[k];
+                    sum += (T)Cheby(n, k) * powers[k];
                 mat.m[i][n] = sum;
             }
         }
@@ -115,21 +123,22 @@ public:
          * place as the absolute error! */
         for (int i = 0; i < ORDER + 1; i++)
         {
-            struct { real value, error; } left, right, mid;
+            struct { T value, error; } left, right, mid;
 
             left.value = control[i];
             left.error = ChebyEval(left.value) - Value(left.value);
             right.value = control[i + 1];
             right.error = ChebyEval(right.value) - Value(right.value);
 
-            static real limit = real::R_1 >> 500;
+            static T limit = (T)1 >> 500;
+            static T zero = (T)0;
             while (fabs(left.value - right.value) > limit)
             {
                 mid.value = (left.value + right.value) >> 1;
                 mid.error = ChebyEval(mid.value) - Value(mid.value);
 
-                if ((left.error < real::R_0 && mid.error < real::R_0)
-                     || (left.error > real::R_0 && mid.error > real::R_0))
+                if ((left.error < zero && mid.error < zero)
+                     || (left.error > zero && mid.error > zero))
                     left = mid;
                 else
                     right = mid;
@@ -146,11 +155,11 @@ public:
         /* Find ORDER + 2 extrema of the error function. We need to
          * compute the relative error, since its extrema are at slightly
          * different locations than the absolute error’s. */
-        real final = 0;
+        T final = 0;
 
         for (int i = 0; i < ORDER + 2; i++)
         {
-            real a = -1, b = 1;
+            T a = -1, b = 1;
             if (i > 0)
                 a = zeroes[i - 1];
             if (i < ORDER + 1)
@@ -158,14 +167,14 @@ public:
 
             for (;;)
             {
-                real c = a, delta = (b - a) >> 3;
-                real maxerror = 0;
-                real maxweight = 0;
+                T c = a, delta = (b - a) >> 3;
+                T maxerror = 0;
+                T maxweight = 0;
                 int best = -1;
                 for (int k = 1; k <= 7; k++)
                 {
-                    real error = ChebyEval(c) - Value(c);
-                    real weight = Weight(c);
+                    T error = ChebyEval(c) - Value(c);
+                    T weight = Weight(c);
                     if (fabs(error * maxweight) >= fabs(maxerror * weight))
                     {
                         maxerror = error;
@@ -175,12 +184,12 @@ public:
                     c += delta;
                 }
 
-                b = a + (real)(best + 1) * delta;
-                a = a + (real)(best - 1) * delta;
+                b = a + (T)(best + 1) * delta;
+                a = a + (T)(best - 1) * delta;
 
-                if (b - a < (real)1e-18)
+                if (b - a < (T)1e-18)
                 {
-                    real e = maxerror / maxweight;
+                    T e = maxerror / maxweight;
                     if (e > final)
                         final = e;
                     control[i] = (a + b) >> 1;
@@ -196,17 +205,17 @@ public:
     void Step()
     {
         /* Pick up x_i where error will be 0 and compute f(x_i) */
-        real fxn[ORDER + 2];
+        T fxn[ORDER + 2];
         for (int i = 0; i < ORDER + 2; i++)
             fxn[i] = Value(control[i]);
 
         /* We build a matrix of Chebishev evaluations: row i contains the
          * evaluations of x_i for polynomial order n = 0, 1, ... */
-        Matrix<ORDER + 2> mat;
+        lol::Mat<ORDER + 2, T> mat;
         for (int i = 0; i < ORDER + 2; i++)
         {
             /* Compute the powers of x_i */
-            real powers[ORDER + 1];
+            T powers[ORDER + 1];
             powers[0] = 1.0;
             for (int n = 1; n < ORDER + 1; n++)
                  powers[n] = powers[n - 1] * control[i];
@@ -214,9 +223,9 @@ public:
             /* Compute the Chebishev evaluations at x_i */
             for (int n = 0; n < ORDER + 1; n++)
             {
-                real sum = 0.0;
+                T sum = 0.0;
                 for (int k = 0; k < ORDER + 1; k++)
-                    sum += (real)Cheby(n, k) * powers[k];
+                    sum += (T)Cheby(n, k) * powers[k];
                 mat.m[i][n] = sum;
             }
             if (i & 1)
@@ -237,7 +246,7 @@ public:
         }
 
         /* Compute the error */
-        real error = 0;
+        T error = 0;
         for (int i = 0; i < ORDER + 2; i++)
             error += mat.m[ORDER + 1][i] * fxn[i];
     }
@@ -264,31 +273,31 @@ public:
 
         /* Transform Chebyshev polynomial weights into powers of X^i
          * in the [-1..1] range. */
-        real bn[ORDER + 1];
+        T bn[ORDER + 1];
 
         for (int i = 0; i < ORDER + 1; i++)
         {
             bn[i] = 0;
             for (int j = 0; j < ORDER + 1; j++)
-                bn[i] += coeff[j] * (real)Cheby(j, i);
+                bn[i] += coeff[j] * (T)Cheby(j, i);
         }
 
         /* Transform a polynomial in the [-1..1] range into a polynomial
          * in the [a..b] range. */
-        real k1p[ORDER + 1], k2p[ORDER + 1];
-        real an[ORDER + 1];
+        T k1p[ORDER + 1], k2p[ORDER + 1];
+        T an[ORDER + 1];
 
         for (int i = 0; i < ORDER + 1; i++)
         {
-            k1p[i] = i ? k1p[i - 1] * m_invk1 : real::R_1;
-            k2p[i] = i ? k2p[i - 1] * m_invk2 : real::R_1;
+            k1p[i] = i ? k1p[i - 1] * m_invk1 : (T)1;
+            k2p[i] = i ? k2p[i - 1] * m_invk2 : (T)1;
         }
 
         for (int i = 0; i < ORDER + 1; i++)
         {
             an[i] = 0;
             for (int j = i; j < ORDER + 1; j++)
-                an[i] += (real)Comb(j, i) * k1p[j - i] * bn[j];
+                an[i] += (T)Comb(j, i) * k1p[j - i] * bn[j];
             an[i] *= k2p[i];
         }
 
@@ -297,33 +306,35 @@ public:
         {
             if (j)
                 printf("+");
-            printf("x^%i*", j);
+            printf("x**%i*", j);
             an[j].print(40);
         }
         printf("\n");
     }
 
-    real Value(real const &x)
+    T Value(T const &x)
     {
         return m_func(x * m_k2 + m_k1);
     }
 
-    real Weight(real const &x)
+    T Weight(T const &x)
     {
         return m_weight(x * m_k2 + m_k1);
     }
 
     /* ORDER + 1 Chebyshev coefficients and 1 error value */
-    real coeff[ORDER + 2];
+    T coeff[ORDER + 2];
     /* ORDER + 1 zeroes of the error function */
-    real zeroes[ORDER + 1];
+    T zeroes[ORDER + 1];
     /* ORDER + 2 control points */
-    real control[ORDER + 2];
+    T control[ORDER + 2];
 
 private:
     RealFunc *m_func, *m_weight;
-    real m_k1, m_k2, m_invk1, m_invk2;
+    T m_k1, m_k2, m_invk1, m_invk2;
 };
 
-#endif /* __REMEZ_SOLVER_H__ */
+} /* namespace lol */
+
+#endif /* __LOL_MATH_REMEZ_H__ */
 

src/platform/nacl/naclapp.h

@@ -16,7 +16,7 @@
 #if !defined __LOL_NACLAPP_H__
 #define __LOL_NACLAPP_H__
 
-#include "matrix.h"
+#include "lol/math/matrix.h"
 
 namespace lol
 {

src/platform/ps3/ps3app.h

@@ -16,7 +16,7 @@
 #if !defined __LOL_PS3APP_H__
 #define __LOL_PS3APP_H__
 
-#include "matrix.h"
+#include "lol/math/matrix.h"
 
 namespace lol
 {

src/platform/sdl/sdlapp.h

@@ -16,7 +16,7 @@
 #if !defined __LOL_SDLAPP_H__
 #define __LOL_SDLAPP_H__
 
-#include "matrix.h"
+#include "lol/math/matrix.h"
 
 namespace lol
 {

test/math/Makefile.am

@@ -22,7 +22,7 @@ poly_CPPFLAGS = @LOL_CFLAGS@ @PIPI_CFLAGS@
 poly_LDFLAGS = $(top_builddir)/src/liblol.a @LOL_LIBS@ @PIPI_LIBS@
 poly_DEPENDENCIES = $(top_builddir)/src/liblol.a
 
-remez_SOURCES = remez.cpp remez-matrix.h remez-solver.h
+remez_SOURCES = remez.cpp
 remez_CPPFLAGS = @LOL_CFLAGS@ @PIPI_CFLAGS@
 remez_LDFLAGS = $(top_builddir)/src/liblol.a @LOL_LIBS@ @PIPI_LIBS@
 remez_DEPENDENCIES = $(top_builddir)/src/liblol.a

test/math/remez-matrix.h

@@ -1,97 +0,0 @@
-//
-// Lol Engine - Sample math program: Chebyshev polynomials
-//
-// Copyright: (c) 2005-2011 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://sam.zoy.org/projects/COPYING.WTFPL for more details.
-//
-
-#if !defined __REMEZ_MATRIX_H__
-#define __REMEZ_MATRIX_H__
-
-template<int N> struct Matrix
-{
-    inline Matrix() {}
-
-    Matrix(real x)
-    {
-        for (int j = 0; j < N; j++)
-            for (int i = 0; i < N; i++)
-                if (i == j)
-                    m[i][j] = x;
-                else
-                    m[i][j] = 0;
-    }
-
-    /* Naive matrix inversion */
-    Matrix<N> inv() const
-    {
-        Matrix a = *this, b((real)1.0);
-
-        /* Inversion method: iterate through all columns and make sure
-         * all the terms are 1 on the diagonal and 0 everywhere else */
-        for (int i = 0; i < N; i++)
-        {
-            /* If the expected coefficient is zero, add one of
-             * the other lines. The first we meet will do. */
-            if ((double)a.m[i][i] == 0.0)
-            {
-                for (int j = i + 1; j < N; j++)
-                {
-                    if ((double)a.m[i][j] == 0.0)
-                        continue;
-                    /* Add row j to row i */
-                    for (int n = 0; n < N; n++)
-                    {
-                        a.m[n][i] += a.m[n][j];
-                        b.m[n][i] += b.m[n][j];
-                    }
-                    break;
-                }
-            }
-
-            /* Now we know the diagonal term is non-zero. Get its inverse
-             * and use that to nullify all other terms in the column */
-            real x = (real)1.0 / a.m[i][i];
-            for (int j = 0; j < N; j++)
-            {
-                if (j == i)
-                    continue;
-                real mul = x * a.m[i][j];
-                for (int n = 0; n < N; n++)
-                {
-                    a.m[n][j] -= mul * a.m[n][i];
-                    b.m[n][j] -= mul * b.m[n][i];
-                }
-            }
-
-            /* Finally, ensure the diagonal term is 1 */
-            for (int n = 0; n < N; n++)
-            {
-                a.m[n][i] *= x;
-                b.m[n][i] *= x;
-            }
-        }
-
-        return b;
-    }
-
-    void print() const
-    {
-        using std::printf;
-
-        for (int j = 0; j < N; j++)
-        {
-            for (int i = 0; i < N; i++)
-                printf("%9.5f ", (double)m[j][i]);
-            printf("\n");
-        }
-    }
-
-    real m[N][N];
-};
-
-#endif /* __REMEZ_MATRIX_H__ */
-

test/math/remez.cpp

@@ -12,19 +12,19 @@
 #   include "config.h"
 #endif
 
-#include <cstring>
 #include <cstdio>
+#include <cstdlib>
 
 #if USE_SDL && defined __APPLE__
 #   include <SDL_main.h>
 #endif
 
-#include "core.h"
+#include "lol/math/real.h"
+#include "lol/math/matrix.h"
+#include "lol/math/remez.h"
 
 using lol::real;
-
-#include "remez-matrix.h"
-#include "remez-solver.h"
+using lol::RemezSolver;
 
 /* The function we want to approximate */
 real myfun(real const &y)
@@ -41,7 +41,7 @@ real myerr(real const &y)
 
 int main(int argc, char **argv)
 {
-    RemezSolver<2> solver;
+    RemezSolver<2, real> solver;
     solver.Run(real::R_1 >> 400, real::R_PI_2 * real::R_PI_2, myfun, myerr, 40);
 
     return EXIT_SUCCESS;