From 48bf48a4e4408c0fa4d2adffbc8f86b406895878 Mon Sep 17 00:00:00 2001 From: Sam Hocevar Date: Thu, 29 Dec 2011 18:36:48 +0000 Subject: [PATCH] math: move the Remez algorithm implementation to the core. --- src/Makefile.am | 4 +- src/core.h | 4 +- src/eglapp.h | 2 +- src/image/image.h | 2 +- src/input.h | 2 +- src/{ => lol/math}/matrix.h | 84 ++++++++++++- src/{ => lol/math}/real.h | 6 +- .../remez-solver.h => src/lol/math/remez.h | 117 ++++++++++-------- src/platform/nacl/naclapp.h | 2 +- src/platform/ps3/ps3app.h | 2 +- src/platform/sdl/sdlapp.h | 2 +- test/math/Makefile.am | 2 +- test/math/remez-matrix.h | 97 --------------- test/math/remez.cpp | 12 +- 14 files changed, 165 insertions(+), 173 deletions(-) rename src/{ => lol/math}/matrix.h (87%) rename src/{ => lol/math}/real.h (97%) rename test/math/remez-solver.h => src/lol/math/remez.h (75%) delete mode 100644 test/math/remez-matrix.h diff --git a/src/Makefile.am b/src/Makefile.am index 38789e6b..6a567616 100644 --- a/src/Makefile.am +++ b/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 \ diff --git a/src/core.h b/src/core.h index 4db4d973..0ce798b4 100644 --- a/src/core.h +++ b/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" diff --git a/src/eglapp.h b/src/eglapp.h index 7364d775..25c739fd 100644 --- a/src/eglapp.h +++ b/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 { diff --git a/src/image/image.h b/src/image/image.h index 458c4794..d4e85d50 100644 --- a/src/image/image.h +++ b/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 { diff --git a/src/input.h b/src/input.h index 39df90ca..5768e2fe 100644 --- a/src/input.h +++ b/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 { diff --git a/src/matrix.h b/src/lol/math/matrix.h similarity index 87% rename from src/matrix.h rename to src/lol/math/matrix.h index 8b2c94ab..fc520516 100644 --- a/src/matrix.h +++ b/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 #if !defined __ANDROID__ @@ -24,6 +24,9 @@ namespace lol { +class half; +class real; + #define VECTOR_TYPES(tname, suffix) \ template struct tname; \ typedef tname f16##suffix; \ @@ -587,7 +590,82 @@ template struct Mat4 Vec4 v[4]; }; +/* + * Arbitrarily-sized square matrices; for now this only supports + * naive inversion and is used for the Remez inversion method. + */ + +template struct Mat +{ + inline Mat() {} + + 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 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__ diff --git a/src/real.h b/src/lol/math/real.h similarity index 97% rename from src/real.h rename to src/lol/math/real.h index acaaecfe..488a6a6d 100644 --- a/src/real.h +++ b/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 @@ -147,5 +147,5 @@ private: } /* namespace lol */ -#endif // __LOL_REAL_H__ +#endif // __LOL_MATH_REAL_H__ diff --git a/test/math/remez-solver.h b/src/lol/math/remez.h similarity index 75% rename from test/math/remez-solver.h rename to src/lol/math/remez.h index d3d7cc08..4714151f 100644 --- a/test/math/remez-solver.h +++ b/src/lol/math/remez.h @@ -1,26 +1,34 @@ // -// Lol Engine - Sample math program: Chebyshev polynomials +// Lol Engine // -// Copyright: (c) 2005-2011 Sam Hocevar +// Copyright: (c) 2010-2011 Sam Hocevar // 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 class RemezSolver +namespace lol +{ + +template 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 mat; + lol::Mat 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 mat; + lol::Mat 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__ */ diff --git a/src/platform/nacl/naclapp.h b/src/platform/nacl/naclapp.h index d654465a..ddf260a0 100644 --- a/src/platform/nacl/naclapp.h +++ b/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 { diff --git a/src/platform/ps3/ps3app.h b/src/platform/ps3/ps3app.h index cd02d4b2..227d3b7e 100644 --- a/src/platform/ps3/ps3app.h +++ b/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 { diff --git a/src/platform/sdl/sdlapp.h b/src/platform/sdl/sdlapp.h index d7195e62..06232042 100644 --- a/src/platform/sdl/sdlapp.h +++ b/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 { diff --git a/test/math/Makefile.am b/test/math/Makefile.am index a245a8e1..3a16e700 100644 --- a/test/math/Makefile.am +++ b/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 diff --git a/test/math/remez-matrix.h b/test/math/remez-matrix.h deleted file mode 100644 index 1832ffd9..00000000 --- a/test/math/remez-matrix.h +++ /dev/null @@ -1,97 +0,0 @@ -// -// Lol Engine - Sample math program: Chebyshev polynomials -// -// Copyright: (c) 2005-2011 Sam Hocevar -// 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 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 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__ */ - diff --git a/test/math/remez.cpp b/test/math/remez.cpp index 1375a136..4f5dc276 100644 --- a/test/math/remez.cpp +++ b/test/math/remez.cpp @@ -12,19 +12,19 @@ # include "config.h" #endif -#include #include +#include #if USE_SDL && defined __APPLE__ # include #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;