From 90706c8d4464eb24043dc41e34eed2f827bb92ae Mon Sep 17 00:00:00 2001 From: Sam Hocevar Date: Tue, 10 Dec 2013 09:32:45 +0000 Subject: [PATCH] lolremez: move some LolRemez matrix functions out of the engine. --- src/lol/base/assert.h | 12 +- src/lol/math/vector.h | 75 ------- tools/lolremez/Makefile.am | 3 +- tools/lolremez/lolremez.cpp | 10 +- tools/lolremez/lolremez.h | 376 ----------------------------------- tools/lolremez/solver.cpp | 377 ++++++++++++++++++++++++++++++++++++ tools/lolremez/solver.h | 57 ++++++ 7 files changed, 447 insertions(+), 463 deletions(-) delete mode 100644 tools/lolremez/lolremez.h create mode 100644 tools/lolremez/solver.cpp create mode 100644 tools/lolremez/solver.h diff --git a/src/lol/base/assert.h b/src/lol/base/assert.h index 38e79eef..bb0cde54 100644 --- a/src/lol/base/assert.h +++ b/src/lol/base/assert.h @@ -31,11 +31,11 @@ extern void DumpStack(); * on implementing __debugbreak() on POSIX systems. */ static inline void DebugAbort() { - DumpStack(); + lol::DumpStack(); #if defined _WIN32 __debugbreak(); #endif - Abort(); + lol::Abort(); } #define LOL_CALL(macro, args) macro args @@ -111,13 +111,13 @@ static inline void DebugAbort() */ #define LOL_ERROR_1(t) \ - Log::Error("assertion at %s:%d: %s\n", __FILE__, __LINE__, #t) + lol::Log::Error("assertion at %s:%d: %s\n", __FILE__, __LINE__, #t) #define LOL_ERROR_2(t, s) \ - Log::Error("assertion at %s:%d: %s\n", __FILE__, __LINE__, s) + lol::Log::Error("assertion at %s:%d: %s\n", __FILE__, __LINE__, s) #define LOL_ERROR_3(t, s, ...) \ - Log::Error("assertion at %s:%d: " s "\n", __FILE__, __LINE__, __VA_ARGS__) + lol::Log::Error("assertion at %s:%d: " s "\n", __FILE__, __LINE__, __VA_ARGS__) #if FINAL_RELEASE # define ASSERT(...) UNUSED(LOL_CALL(LOL_1ST, (__VA_ARGS__))) @@ -128,7 +128,7 @@ static inline void DebugAbort() LOL_CALL(LOL_CAT(LOL_ERROR_, LOL_CALL(LOL_COUNT_TO_3, \ (__VA_ARGS__))), \ (__VA_ARGS__)); \ - DebugAbort(); \ + lol::DebugAbort(); \ } #endif diff --git a/src/lol/math/vector.h b/src/lol/math/vector.h index 21d4f1eb..7afac1c3 100644 --- a/src/lol/math/vector.h +++ b/src/lol/math/vector.h @@ -1918,81 +1918,6 @@ template Mat2 inverse(Mat2 const &); template Mat3 inverse(Mat3 const &); template Mat4 inverse(Mat4 const &); -/* - * 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_MATH_VECTOR_H__ diff --git a/tools/lolremez/Makefile.am b/tools/lolremez/Makefile.am index bba1426d..e1ecd5fc 100644 --- a/tools/lolremez/Makefile.am +++ b/tools/lolremez/Makefile.am @@ -5,7 +5,8 @@ EXTRA_DIST += NEWS.txt lolremez.sln remez.vcxproj remez.vcxproj.filters noinst_PROGRAMS = lolremez -lolremez_SOURCES = lolremez.cpp +lolremez_SOURCES = \ + lolremez.cpp solver.cpp solver.h matrix.h lolremez_CPPFLAGS = $(AM_CPPFLAGS) lolremez_DEPENDENCIES = @LOL_DEPS@ diff --git a/tools/lolremez/lolremez.cpp b/tools/lolremez/lolremez.cpp index f4916a48..b6fc6b29 100644 --- a/tools/lolremez/lolremez.cpp +++ b/tools/lolremez/lolremez.cpp @@ -1,5 +1,5 @@ // -// Lol Engine - Sample math program: Chebyshev polynomials +// LolRemez - Remez algorithm implementation // // Copyright: (c) 2005-2013 Sam Hocevar // This program is free software; you can redistribute it and/or @@ -16,10 +16,9 @@ #include -#include "lolremez.h" +#include "solver.h" using lol::real; -using lol::RemezSolver; /* See the tutorial at http://lolengine.net/wiki/doc/maths/remez */ real f(real const &x) @@ -36,8 +35,9 @@ int main(int argc, char **argv) { UNUSED(argc, argv); - RemezSolver<4, real> solver; - solver.Run(-1, 1, f, g, 40); + RemezSolver solver; + solver.Run(4, 40, -1, 1, f, g); + return 0; } diff --git a/tools/lolremez/lolremez.h b/tools/lolremez/lolremez.h deleted file mode 100644 index 1675436d..00000000 --- a/tools/lolremez/lolremez.h +++ /dev/null @@ -1,376 +0,0 @@ -// -// Lol Engine -// -// 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://www.wtfpl.net/ for more details. -// - -// -// The RemezSolver class -// --------------------- -// - -#if !defined __LOL_MATH_REMEZ_H__ -#define __LOL_MATH_REMEZ_H__ - -#include - -#include "lol/math/vector.h" - -namespace lol -{ - -template class RemezSolver -{ -public: - typedef T RealFunc(T const &x); - - RemezSolver() - { - } - - void Run(T a, T b, RealFunc *func, RealFunc *weight, int decimals) - { - using std::printf; - - m_func = func; - m_weight = weight; - m_k1 = (b + a) / 2; - m_k2 = (b - a) / 2; - m_invk2 = re(m_k2); - m_invk1 = -m_k1 * m_invk2; - m_decimals = decimals; - m_epsilon = pow((T)10, (T)-(decimals + 2)); - - Init(); - - PrintPoly(); - - T error = -1; - - for (int n = 0; ; n++) - { - T newerror = FindExtrema(); - printf("Step %i error: ", n); - newerror.print(m_decimals); - printf("\n"); - - Step(); - - if (error >= (T)0 && fabs(newerror - error) < error * m_epsilon) - break; - error = newerror; - - PrintPoly(); - - FindZeroes(); - } - - PrintPoly(); - } - - inline void Run(T a, T b, RealFunc *func, int decimals) - { - Run(a, b, func, nullptr, decimals); - } - - T EvalCheby(T const &x) - { - T ret = 0.0, xn = 1.0; - - for (int i = 0; i < ORDER + 1; i++) - { - T mul = 0; - for (int j = 0; j < ORDER + 1; j++) - mul += coeff[j] * (T)Cheby(j, i); - ret += mul * xn; - xn *= x; - } - - return ret; - } - - void Init() - { - /* Pick up x_i where error will be 0 and compute f(x_i) */ - T fxn[ORDER + 1]; - for (int i = 0; i < ORDER + 1; i++) - { - zeroes[i] = (T)(2 * i - ORDER) / (T)(ORDER + 1); - fxn[i] = EvalFunc(zeroes[i]); - } - - /* We build a matrix of Chebishev evaluations: row i contains the - * evaluations of x_i for polynomial order n = 0, 1, ... */ - lol::Mat mat; - for (int i = 0; i < ORDER + 1; i++) - { - /* Compute the powers of x_i */ - T powers[ORDER + 1]; - powers[0] = 1.0; - for (int n = 1; n < ORDER + 1; n++) - powers[n] = powers[n - 1] * zeroes[i]; - - /* Compute the Chebishev evaluations at x_i */ - for (int n = 0; n < ORDER + 1; n++) - { - T sum = 0.0; - for (int k = 0; k < ORDER + 1; k++) - sum += (T)Cheby(n, k) * powers[k]; - mat.m[i][n] = sum; - } - } - - /* Solve the system */ - mat = mat.inv(); - - /* Compute interpolation coefficients */ - for (int j = 0; j < ORDER + 1; j++) - { - coeff[j] = 0; - for (int i = 0; i < ORDER + 1; i++) - coeff[j] += mat.m[j][i] * fxn[i]; - } - } - - void FindZeroes() - { - /* Find ORDER + 1 zeroes of the error function. No need to - * compute the relative error: its zeroes are at the same - * place as the absolute error! */ - for (int i = 0; i < ORDER + 1; i++) - { - struct { T value, error; } left, right, mid; - - left.value = control[i]; - left.error = EvalCheby(left.value) - EvalFunc(left.value); - right.value = control[i + 1]; - right.error = EvalCheby(right.value) - EvalFunc(right.value); - - static T limit = ldexp((T)1, -500); - static T zero = (T)0; - while (fabs(left.value - right.value) > limit) - { - mid.value = (left.value + right.value) / 2; - mid.error = EvalCheby(mid.value) - EvalFunc(mid.value); - - if ((left.error <= zero && mid.error <= zero) - || (left.error >= zero && mid.error >= zero)) - left = mid; - else - right = mid; - } - - zeroes[i] = mid.value; - } - } - - real FindExtrema() - { - using std::printf; - - /* 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. */ - T final = 0; - - for (int i = 0; i < ORDER + 2; i++) - { - T a = -1, b = 1; - if (i > 0) - a = zeroes[i - 1]; - if (i < ORDER + 1) - b = zeroes[i]; - - for (int round = 0; ; round++) - { - T maxerror = 0, maxweight = 0; - int best = -1; - - T c = a, delta = (b - a) / 4; - for (int k = 0; k <= 4; k++) - { - if (round == 0 || (k & 1)) - { - T error = fabs(EvalCheby(c) - EvalFunc(c)); - T weight = fabs(Weight(c)); - /* if error/weight >= maxerror/maxweight */ - if (error * maxweight >= maxerror * weight) - { - maxerror = error; - maxweight = weight; - best = k; - } - } - c += delta; - } - - switch (best) - { - case 0: - b = a + delta * 2; - break; - case 4: - a = b - delta * 2; - break; - default: - b = a + delta * (best + 1); - a = a + delta * (best - 1); - break; - } - - if (delta < m_epsilon) - { - T e = fabs(maxerror / maxweight); - if (e > final) - final = e; - control[i] = (a + b) / 2; - break; - } - } - } - - return final; - } - - void Step() - { - /* Pick up x_i where error will be 0 and compute f(x_i) */ - T fxn[ORDER + 2]; - for (int i = 0; i < ORDER + 2; i++) - fxn[i] = EvalFunc(control[i]); - - /* We build a matrix of Chebishev evaluations: row i contains the - * evaluations of x_i for polynomial order n = 0, 1, ... */ - lol::Mat mat; - for (int i = 0; i < ORDER + 2; i++) - { - /* Compute the powers of x_i */ - T powers[ORDER + 1]; - powers[0] = 1.0; - for (int n = 1; n < ORDER + 1; n++) - powers[n] = powers[n - 1] * control[i]; - - /* Compute the Chebishev evaluations at x_i */ - for (int n = 0; n < ORDER + 1; n++) - { - T sum = 0.0; - for (int k = 0; k < ORDER + 1; k++) - sum += (T)Cheby(n, k) * powers[k]; - mat.m[i][n] = sum; - } - if (i & 1) - mat.m[i][ORDER + 1] = fabs(Weight(control[i])); - else - mat.m[i][ORDER + 1] = -fabs(Weight(control[i])); - } - - /* Solve the system */ - mat = mat.inv(); - - /* Compute interpolation coefficients */ - for (int j = 0; j < ORDER + 1; j++) - { - coeff[j] = 0; - for (int i = 0; i < ORDER + 2; i++) - coeff[j] += mat.m[j][i] * fxn[i]; - } - - /* Compute the error */ - T error = 0; - for (int i = 0; i < ORDER + 2; i++) - error += mat.m[ORDER + 1][i] * fxn[i]; - } - - int Cheby(int n, int k) - { - if (k > n || k < 0) - return 0; - if (n <= 1) - return (n ^ k ^ 1) & 1; - return 2 * Cheby(n - 1, k - 1) - Cheby(n - 2, k); - } - - int Comb(int n, int k) - { - if (k == 0 || k == n) - return 1; - return Comb(n - 1, k - 1) + Comb(n - 1, k); - } - - void PrintPoly() - { - using std::printf; - - /* Transform Chebyshev polynomial weights into powers of X^i - * in the [-1..1] range. */ - 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] * (T)Cheby(j, i); - } - - /* Transform a polynomial in the [-1..1] range into a polynomial - * in the [a..b] range. */ - 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 : (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] += (T)Comb(j, i) * k1p[j - i] * bn[j]; - an[i] *= k2p[i]; - } - - printf("Polynomial estimate: "); - for (int j = 0; j < ORDER + 1; j++) - { - if (j) - printf(" + x**%i * ", j); - an[j].print(m_decimals); - } - printf("\n\n"); - } - - T EvalFunc(T const &x) - { - return m_func(x * m_k2 + m_k1); - } - - T Weight(T const &x) - { - if (m_weight) - return m_weight(x * m_k2 + m_k1); - return 1; - } - - /* ORDER + 1 Chebyshev coefficients and 1 error value */ - T coeff[ORDER + 2]; - /* ORDER + 1 zeroes of the error function */ - T zeroes[ORDER + 1]; - /* ORDER + 2 control points */ - T control[ORDER + 2]; - -private: - RealFunc *m_func, *m_weight; - T m_k1, m_k2, m_invk1, m_invk2, m_epsilon; - int m_decimals; -}; - -} /* namespace lol */ - -#endif /* __LOL_MATH_REMEZ_H__ */ - diff --git a/tools/lolremez/solver.cpp b/tools/lolremez/solver.cpp new file mode 100644 index 00000000..fe0ec999 --- /dev/null +++ b/tools/lolremez/solver.cpp @@ -0,0 +1,377 @@ +// +// LolRemez - Remez algorithm implementation +// +// Copyright: (c) 2005-2013 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://www.wtfpl.net/ for more details. +// + +#if defined HAVE_CONFIG_H +# include "config.h" +#endif + +#include "core.h" + +#include + +#include "matrix.h" +#include "solver.h" + +using lol::real; + +/* Some forward declarations first. */ +template<> void RemezSolver::Init(); +template<> void RemezSolver::FindZeroes(); +template<> real RemezSolver::FindExtrema(); +template<> void RemezSolver::Step(); +template<> int RemezSolver::Cheby(int n, int k); +template<> void RemezSolver::PrintPoly(); +template<> real RemezSolver::EvalFunc(real const &x); +template<> real RemezSolver::Weight(real const &x); + + +template<> +void RemezSolver::Run(int order, int decimals, real a, real b, + RemezSolver::RealFunc *func, + RemezSolver::RealFunc *weight) +{ + using std::printf; + + m_order = order; + m_func = func; + m_weight = weight; + m_k1 = (b + a) / 2; + m_k2 = (b - a) / 2; + m_invk2 = re(m_k2); + m_invk1 = -m_k1 * m_invk2; + m_decimals = decimals; + m_epsilon = pow((real)10, (real)-(decimals + 2)); + + Init(); + + PrintPoly(); + + real error = -1; + + for (int n = 0; ; n++) + { + real newerror = FindExtrema(); + printf("Step %i error: ", n); + newerror.print(m_decimals); + printf("\n"); + + Step(); + + if (error >= (real)0 && fabs(newerror - error) < error * m_epsilon) + break; + error = newerror; + + PrintPoly(); + + FindZeroes(); + } + + PrintPoly(); +} + +template<> +real RemezSolver::EvalCheby(real const &x) +{ + real ret = 0.0, xn = 1.0; + + for (int i = 0; i < m_order + 1; i++) + { + real mul = 0; + for (int j = 0; j < m_order + 1; j++) + mul += m_coeff[j] * (real)Cheby(j, i); + ret += mul * xn; + xn *= x; + } + + return ret; +} + +template<> +void RemezSolver::Init() +{ + /* m_order + 1 Chebyshev coefficients, plus 1 error value */ + m_coeff.Resize(m_order + 2); + + /* m_order + 1 zeroes of the error function */ + m_zeroes.Resize(m_order + 1); + + /* m_order + 2 control points */ + m_control.Resize(m_order + 2); + + /* Pick up x_i where error will be 0 and compute f(x_i) */ + real fxn[m_order + 1]; + for (int i = 0; i < m_order + 1; i++) + { + m_zeroes[i] = (real)(2 * i - m_order) / (real)(m_order + 1); + fxn[i] = EvalFunc(m_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(m_order + 1, m_order + 1); + for (int i = 0; i < m_order + 1; i++) + { + /* Compute the powers of x_i */ + real powers[m_order + 1]; + powers[0] = 1.0; + for (int n = 1; n < m_order + 1; n++) + powers[n] = powers[n - 1] * m_zeroes[i]; + + /* Compute the Chebishev evaluations at x_i */ + for (int n = 0; n < m_order + 1; n++) + { + real sum = 0.0; + for (int k = 0; k < m_order + 1; k++) + sum += (real)Cheby(n, k) * powers[k]; + mat.m(i, n) = sum; + } + } + + /* Solve the system */ + mat = mat.inv(); + + /* Compute interpolation coefficients */ + for (int j = 0; j < m_order + 1; j++) + { + m_coeff[j] = 0; + for (int i = 0; i < m_order + 1; i++) + m_coeff[j] += mat.m(j, i) * fxn[i]; + } +} + +template<> +void RemezSolver::FindZeroes() +{ + /* Find m_order + 1 zeroes of the error function. No need to + * compute the relative error: its zeroes are at the same + * place as the absolute error! */ + for (int i = 0; i < m_order + 1; i++) + { + struct { real value, error; } left, right, mid; + + left.value = m_control[i]; + left.error = EvalCheby(left.value) - EvalFunc(left.value); + right.value = m_control[i + 1]; + right.error = EvalCheby(right.value) - EvalFunc(right.value); + + static real limit = ldexp((real)1, -500); + static real zero = (real)0; + while (fabs(left.value - right.value) > limit) + { + mid.value = (left.value + right.value) / 2; + mid.error = EvalCheby(mid.value) - EvalFunc(mid.value); + + if ((left.error <= zero && mid.error <= zero) + || (left.error >= zero && mid.error >= zero)) + left = mid; + else + right = mid; + } + + m_zeroes[i] = mid.value; + } +} + +template<> +real RemezSolver::FindExtrema() +{ + using std::printf; + + /* Find m_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; + + for (int i = 0; i < m_order + 2; i++) + { + real a = -1, b = 1; + if (i > 0) + a = m_zeroes[i - 1]; + if (i < m_order + 1) + b = m_zeroes[i]; + + for (int round = 0; ; round++) + { + real maxerror = 0, maxweight = 0; + int best = -1; + + real c = a, delta = (b - a) / 4; + for (int k = 0; k <= 4; k++) + { + if (round == 0 || (k & 1)) + { + real error = fabs(EvalCheby(c) - EvalFunc(c)); + real weight = fabs(Weight(c)); + /* if error/weight >= maxerror/maxweight */ + if (error * maxweight >= maxerror * weight) + { + maxerror = error; + maxweight = weight; + best = k; + } + } + c += delta; + } + + switch (best) + { + case 0: + b = a + delta * 2; + break; + case 4: + a = b - delta * 2; + break; + default: + b = a + delta * (best + 1); + a = a + delta * (best - 1); + break; + } + + if (delta < m_epsilon) + { + real e = fabs(maxerror / maxweight); + if (e > final) + final = e; + m_control[i] = (a + b) / 2; + break; + } + } + } + + return final; +} + +template<> +void RemezSolver::Step() +{ + /* Pick up x_i where error will be 0 and compute f(x_i) */ + real fxn[m_order + 2]; + for (int i = 0; i < m_order + 2; i++) + fxn[i] = EvalFunc(m_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(m_order + 2, m_order + 2); + for (int i = 0; i < m_order + 2; i++) + { + /* Compute the powers of x_i */ + real powers[m_order + 1]; + powers[0] = 1.0; + for (int n = 1; n < m_order + 1; n++) + powers[n] = powers[n - 1] * m_control[i]; + + /* Compute the Chebishev evaluations at x_i */ + for (int n = 0; n < m_order + 1; n++) + { + real sum = 0.0; + for (int k = 0; k < m_order + 1; k++) + sum += (real)Cheby(n, k) * powers[k]; + mat.m(i, n) = sum; + } + if (i & 1) + mat.m(i, m_order + 1) = fabs(Weight(m_control[i])); + else + mat.m(i, m_order + 1) = -fabs(Weight(m_control[i])); + } + + /* Solve the system */ + mat = mat.inv(); + + /* Compute interpolation coefficients */ + for (int j = 0; j < m_order + 1; j++) + { + m_coeff[j] = 0; + for (int i = 0; i < m_order + 2; i++) + m_coeff[j] += mat.m(j, i) * fxn[i]; + } + + /* Compute the error */ + real error = 0; + for (int i = 0; i < m_order + 2; i++) + error += mat.m(m_order + 1, i) * fxn[i]; +} + +template<> +int RemezSolver::Cheby(int n, int k) +{ + if (k > n || k < 0) + return 0; + if (n <= 1) + return (n ^ k ^ 1) & 1; + return 2 * Cheby(n - 1, k - 1) - Cheby(n - 2, k); +} + +template<> +int RemezSolver::Comb(int n, int k) +{ + if (k == 0 || k == n) + return 1; + return Comb(n - 1, k - 1) + Comb(n - 1, k); +} + +template<> +void RemezSolver::PrintPoly() +{ + using std::printf; + + /* Transform Chebyshev polynomial weights into powers of X^i + * in the [-1..1] range. */ + real bn[m_order + 1]; + + for (int i = 0; i < m_order + 1; i++) + { + bn[i] = 0; + for (int j = 0; j < m_order + 1; j++) + bn[i] += m_coeff[j] * (real)Cheby(j, i); + } + + /* Transform a polynomial in the [-1..1] range into a polynomial + * in the [a..b] range. */ + real k1p[m_order + 1], k2p[m_order + 1]; + real an[m_order + 1]; + + for (int i = 0; i < m_order + 1; i++) + { + k1p[i] = i ? k1p[i - 1] * m_invk1 : (real)1; + k2p[i] = i ? k2p[i - 1] * m_invk2 : (real)1; + } + + for (int i = 0; i < m_order + 1; i++) + { + an[i] = 0; + for (int j = i; j < m_order + 1; j++) + an[i] += (real)Comb(j, i) * k1p[j - i] * bn[j]; + an[i] *= k2p[i]; + } + + printf("Polynomial estimate: "); + for (int j = 0; j < m_order + 1; j++) + { + if (j) + printf(" + x**%i * ", j); + an[j].print(m_decimals); + } + printf("\n\n"); +} + +template<> +real RemezSolver::EvalFunc(real const &x) +{ + return m_func(x * m_k2 + m_k1); +} + +template<> +real RemezSolver::Weight(real const &x) +{ + if (m_weight) + return m_weight(x * m_k2 + m_k1); + return 1; +} + diff --git a/tools/lolremez/solver.h b/tools/lolremez/solver.h new file mode 100644 index 00000000..3e8127a5 --- /dev/null +++ b/tools/lolremez/solver.h @@ -0,0 +1,57 @@ +// +// LolRemez - Remez algorithm implementation +// +// Copyright: (c) 2010-2013 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://www.wtfpl.net/ for more details. +// + +#pragma once + +// +// The RemezSolver class +// --------------------- +// + +#include + +template class RemezSolver +{ +public: + typedef NUMERIC_T RealFunc(NUMERIC_T const &x); + + inline RemezSolver() + : m_order(0) + { + } + + void Run(int order, int decimals, NUMERIC_T a, NUMERIC_T b, + RealFunc *func, RealFunc *weight = nullptr); + + NUMERIC_T EvalCheby(NUMERIC_T const &x); + void Init(); + void FindZeroes(); + NUMERIC_T FindExtrema(); + void Step(); + + int Cheby(int n, int k); + int Comb(int n, int k); + + void PrintPoly(); + NUMERIC_T EvalFunc(NUMERIC_T const &x); + NUMERIC_T Weight(NUMERIC_T const &x); + +private: + int m_order; + + lol::Array m_coeff; + lol::Array m_zeroes; + lol::Array m_control; + + RealFunc *m_func, *m_weight; + NUMERIC_T m_k1, m_k2, m_invk1, m_invk2, m_epsilon; + int m_decimals; +}; +