diff --git a/tools/lolremez/solver.cpp b/tools/lolremez/solver.cpp index e58038ee..d1dac5c7 100644 --- a/tools/lolremez/solver.cpp +++ b/tools/lolremez/solver.cpp @@ -12,9 +12,12 @@ # include "config.h" #endif +#include + #include #include +#include #include "matrix.h" #include "solver.h" @@ -68,34 +71,16 @@ void RemezSolver::Run(real a, real b, PrintPoly(); } -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; -} - 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) */ + /* Initial estimates for the x_i where the error will be zero and + * precompute f(x_i). */ array fxn; for (int i = 0; i < m_order + 1; i++) { @@ -106,33 +91,25 @@ void RemezSolver::Init() /* We build a matrix of Chebishev evaluations: row i contains the * evaluations of x_i for polynomial order n = 0, 1, ... */ LinearSystem system(m_order + 1); - for (int i = 0; i < m_order + 1; i++) + for (int n = 0; n < m_order + 1; n++) { - /* Compute the powers of x_i */ - array powers; - powers.Push(real(1.0)); - for (int n = 1; n < m_order + 1; n++) - powers.Push(powers.Last() * 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]; - system.m(i, n) = sum; - } + auto p = polynomial::chebyshev(n); + for (int i = 0; i < m_order + 1; i++) + system.m(i, n) = p.eval(m_zeroes[i]); } /* Solve the system */ system = system.inv(); - /* Compute interpolation coefficients */ - for (int j = 0; j < m_order + 1; j++) + /* Compute new Chebyshev estimate */ + m_estimate = polynomial(); + for (int n = 0; n < m_order + 1; n++) { - m_coeff[j] = 0; + real weight = 0; for (int i = 0; i < m_order + 1; i++) - m_coeff[j] += system.m(j, i) * fxn[i]; + weight += system.m(n, i) * fxn[i]; + + m_estimate += weight * polynomial::chebyshev(n); } } @@ -168,6 +145,7 @@ void RemezSolver::FindZeroes() } } +/* XXX: this is the real costly function */ real RemezSolver::FindExtrema() { using std::printf; @@ -177,6 +155,9 @@ real RemezSolver::FindExtrema() * different locations than the absolute error’s. */ real final = 0; + m_stats_cheby = m_stats_func = m_stats_weight = 0.f; + int evals = 0, rounds = 0; + for (int i = 0; i < m_order + 2; i++) { real a = -1, b = 1; @@ -185,7 +166,7 @@ real RemezSolver::FindExtrema() if (i < m_order + 1) b = m_zeroes[i]; - for (int round = 0; ; round++) + for (int r = 0; ; ++r, ++rounds) { real maxerror = 0, maxweight = 0; int best = -1; @@ -193,10 +174,11 @@ real RemezSolver::FindExtrema() real c = a, delta = (b - a) / 4; for (int k = 0; k <= 4; k++) { - if (round == 0 || (k & 1)) + if (r == 0 || (k & 1)) { + ++evals; real error = fabs(EvalCheby(c) - EvalFunc(c)); - real weight = fabs(Weight(c)); + real weight = fabs(EvalWeight(c)); /* if error/weight >= maxerror/maxweight */ if (error * maxweight >= maxerror * weight) { @@ -233,6 +215,10 @@ real RemezSolver::FindExtrema() } } + printf("Iterations: Rounds %d Evals %d\n", rounds, evals); + printf("Times: Cheby %f Func %f EvalWeight %f\n", + m_stats_cheby, m_stats_func, m_stats_weight); + return final; } @@ -246,37 +232,32 @@ void RemezSolver::Step() /* We build a matrix of Chebishev evaluations: row i contains the * evaluations of x_i for polynomial order n = 0, 1, ... */ LinearSystem system(m_order + 2); + for (int n = 0; n < m_order + 1; n++) + { + auto p = polynomial::chebyshev(n); + for (int i = 0; i < m_order + 2; i++) + system.m(i, n) = p.eval(m_control[i]); + } + + /* The last line of the system is the oscillating error */ for (int i = 0; i < m_order + 2; i++) { - /* Compute the powers of x_i */ - array powers; - powers.Push(real(1.0)); - for (int n = 1; n < m_order + 1; n++) - powers.Push(powers.Last() * 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]; - system.m(i, n) = sum; - } - if (i & 1) - system.m(i, m_order + 1) = fabs(Weight(m_control[i])); - else - system.m(i, m_order + 1) = -fabs(Weight(m_control[i])); + real error = fabs(EvalWeight(m_control[i])); + system.m(i, m_order + 1) = (i & 1) ? error : -error; } /* Solve the system */ system = system.inv(); - /* Compute interpolation coefficients */ - for (int j = 0; j < m_order + 1; j++) + /* Compute new polynomial estimate */ + m_estimate = polynomial(); + for (int n = 0; n < m_order + 1; n++) { - m_coeff[j] = 0; + real weight = 0; for (int i = 0; i < m_order + 2; i++) - m_coeff[j] += system.m(j, i) * fxn[i]; + weight += system.m(n, i) * fxn[i]; + + m_estimate += weight * polynomial::chebyshev(n); } /* Compute the error */ @@ -285,39 +266,11 @@ void RemezSolver::Step() error += system.m(m_order + 1, i) * fxn[i]; } -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); -} - -int RemezSolver::Comb(int n, int k) -{ - if (k == 0 || k == n) - return 1; - return Comb(n - 1, k - 1) + Comb(n - 1, k); -} - void RemezSolver::PrintPoly() { using std::printf; - /* Transform Chebyshev polynomial weights into powers of X^i - * in the [-1..1] range. */ - array bn; - - for (int i = 0; i < m_order + 1; i++) - { - real tmp = 0; - for (int j = 0; j < m_order + 1; j++) - tmp += m_coeff[j] * (real)Cheby(j, i); - bn.Push(tmp); - } - - /* Transform a polynomial in the [-1..1] range into a polynomial + /* Transform our polynomial in the [-1..1] range into a polynomial * in the [a..b] range. */ array k1p, k2p, an; @@ -327,11 +280,18 @@ void RemezSolver::PrintPoly() k2p.Push(i ? k2p[i - 1] * m_invk2 : (real)1); } + std::function comb = [&](int n, int k) + { + if (k == 0 || k == n) + return 1; + return comb(n - 1, k - 1) + comb(n - 1, k); + }; + for (int i = 0; i < m_order + 1; i++) { real tmp = 0; for (int j = i; j < m_order + 1; j++) - tmp += (real)Comb(j, i) * k1p[j - i] * bn[j]; + tmp += (real)comb(j, i) * k1p[j - i] * m_estimate[j]; an.Push(tmp * k2p[i]); } @@ -345,15 +305,27 @@ void RemezSolver::PrintPoly() printf("\n\n"); } +real RemezSolver::EvalCheby(real const &x) +{ + Timer t; + real ret = m_estimate.eval(x); + m_stats_cheby += t.Get(); + return ret; +} + real RemezSolver::EvalFunc(real const &x) { - return m_func(x * m_k2 + m_k1); + Timer t; + real ret = m_func(x * m_k2 + m_k1); + m_stats_func += t.Get(); + return ret; } -real RemezSolver::Weight(real const &x) +real RemezSolver::EvalWeight(real const &x) { - if (m_weight) - return m_weight(x * m_k2 + m_k1); - return 1; + Timer t; + real ret = m_weight ? m_weight(x * m_k2 + m_k1) : real(1); + m_stats_weight += t.Get(); + return ret; } diff --git a/tools/lolremez/solver.h b/tools/lolremez/solver.h index 45980814..3de110d5 100644 --- a/tools/lolremez/solver.h +++ b/tools/lolremez/solver.h @@ -28,28 +28,31 @@ public: RealFunc *func, RealFunc *weight = nullptr); private: - lol::real EvalCheby(lol::real const &x); void Init(); void FindZeroes(); lol::real FindExtrema(); void Step(); - - int Cheby(int n, int k); - int Comb(int n, int k); - void PrintPoly(); + + lol::real EvalCheby(lol::real const &x); lol::real EvalFunc(lol::real const &x); - lol::real Weight(lol::real const &x); + lol::real EvalWeight(lol::real const &x); private: int m_order; - lol::array m_coeff; + lol::polynomial m_estimate; + lol::array m_zeroes; lol::array m_control; RealFunc *m_func, *m_weight; lol::real m_k1, m_k2, m_invk1, m_invk2, m_epsilon; int m_decimals; + + /* Statistics */ + float m_stats_cheby; + float m_stats_func; + float m_stats_weight; };