From 1d0797efa5ac4c5487b9db871d24a8fbe3266d6d Mon Sep 17 00:00:00 2001 From: Sam Hocevar Date: Mon, 10 Oct 2011 01:33:15 +0000 Subject: [PATCH] test: various improvements to the Remez exchange solver. --- test/math/remez-solver.h | 96 ++++++++++++++++++++++------------------ test/math/remez.cpp | 18 +++----- 2 files changed, 58 insertions(+), 56 deletions(-) diff --git a/test/math/remez-solver.h b/test/math/remez-solver.h index c926bff9..4d3f7974 100644 --- a/test/math/remez-solver.h +++ b/test/math/remez-solver.h @@ -20,10 +20,10 @@ public: { } - void Run(real a, real b, RealFunc *func, RealFunc *error, int steps) + void Run(real a, real b, RealFunc *func, RealFunc *weight, int steps) { m_func = func; - m_error = error; + m_weight = weight; m_k1 = (b + a) >> 1; m_k2 = (b - a) >> 1; m_invk2 = re(m_k2); @@ -35,7 +35,7 @@ public: for (int n = 0; n < steps; n++) { - FindError(); + FindExtrema(); Step(); PrintPoly(); @@ -43,7 +43,7 @@ public: FindZeroes(); } - FindError(); + FindExtrema(); Step(); PrintPoly(); @@ -110,37 +110,40 @@ public: 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++) { - real a = control[i]; - real ea = ChebyEval(a) - Value(a); - real b = control[i + 1]; - real eb = ChebyEval(b) - Value(b); + struct { real value, error; } left, right, mid; - while (fabs(a - b) > (real)1e-140) + 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; + while (fabs(left.value - right.value) > limit) { - real c = (a + b) * (real)0.5; - real ec = ChebyEval(c) - Value(c); + mid.value = (left.value + right.value) >> 1; + mid.error = ChebyEval(mid.value) - Value(mid.value); - if ((ea < (real)0 && ec < (real)0) - || (ea > (real)0 && ec > (real)0)) - { - a = c; - ea = ec; - } + if ((left.error < real::R_0 && mid.error < real::R_0) + || (left.error > real::R_0 && mid.error > real::R_0)) + left = mid; else - { - b = c; - eb = ec; - } + right = mid; } - zeroes[i] = a; + zeroes[i] = mid.value; } } - void FindError() + void FindExtrema() { + /* 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; for (int i = 0; i < ORDER + 2; i++) @@ -154,34 +157,33 @@ public: printf("Error for [%g..%g]: ", (double)a, (double)b); for (;;) { - real c = a, delta = (b - a) / (real)10.0; + real c = a, delta = (b - a) >> 3; real maxerror = 0; + real maxweight = 0; int best = -1; - for (int k = 0; k <= 10; k++) + for (int k = 1; k <= 7; k++) { - real e = fabs(ChebyEval(c) - Value(c)); - if (e > maxerror) + real error = ChebyEval(c) - Value(c); + real weight = Weight(c); + if (fabs(error * maxweight) >= fabs(maxerror * weight)) { - maxerror = e; + maxerror = error; + maxweight = weight; best = k; } c += delta; } - if (best == 0) - best = 1; - if (best == 10) - best = 9; - b = a + (real)(best + 1) * delta; a = a + (real)(best - 1) * delta; - if (b - a < (real)1e-15) + if (b - a < (real)1e-18) { - if (maxerror > final) - final = maxerror; - control[i] = (a + b) * (real)0.5; - printf("%g (at %g)\n", (double)maxerror, (double)control[i]); + real e = maxerror / maxweight; + if (e > final) + final = e; + control[i] = (a + b) >> 1; + printf("%g (at %g)\n", (double)e, (double)control[i]); break; } } @@ -218,9 +220,9 @@ public: mat.m[i][n] = sum; } if (i & 1) - mat.m[i][ORDER + 1] = fabs(Error(control[i])); + mat.m[i][ORDER + 1] = fabs(Weight(control[i])); else - mat.m[i][ORDER + 1] = -fabs(Error(control[i])); + mat.m[i][ORDER + 1] = -fabs(Weight(control[i])); } /* Solve the system */ @@ -288,8 +290,14 @@ public: an[i] *= k2p[i]; } + printf("Final polynomial:\n"); for (int j = 0; j < ORDER + 1; j++) - printf("%s%18.16gx^%i", j && (an[j] >= real::R_0) ? "+" : "", (double)an[j], j); + { + if (j) + printf("+"); + printf("x^%i*", j); + an[j].print(40); + } printf("\n"); } @@ -298,9 +306,9 @@ public: return m_func(x * m_k2 + m_k1); } - real Error(real const &x) + real Weight(real const &x) { - return m_error(x * m_k2 + m_k1); + return m_weight(x * m_k2 + m_k1); } /* ORDER + 1 Chebyshev coefficients and 1 error value */ @@ -311,7 +319,7 @@ public: real control[ORDER + 2]; private: - RealFunc *m_func, *m_error; + RealFunc *m_func, *m_weight; real m_k1, m_k2, m_invk1, m_invk2; }; diff --git a/test/math/remez.cpp b/test/math/remez.cpp index 9469d916..1f0e0832 100644 --- a/test/math/remez.cpp +++ b/test/math/remez.cpp @@ -26,26 +26,20 @@ using namespace std; /* The function we want to approximate */ static real myfun(real const &x) { - static real const a0 = real::R_1; - static real const a1 = real(-11184811) >> 26; - static real const b1 = real(-1) / real(6); - static real const b2 = real(1) / real(120); real y = sqrt(x); - return sin(y) / y; + return (sin(y) - y) / (x * y); } static real myerr(real const &x) { - return real::R_1; + real y = sqrt(x); + return re(x * y); } -int main(int argc, char **argv) +int main(void) { - RemezSolver<6> solver; - //solver.Run(0, 1, myfun, myfun, 15); - solver.Run(exp((real)-100), real::R_PI * real::R_PI >> 2, myfun, myfun, 15); - //solver.Run(-1, 1, myfun, myfun, 15); - //solver.Run(0, real::R_PI * real::R_PI >> 4, myfun, myfun, 15); + RemezSolver<4> solver; + solver.Run(real::R_1 >> 400, real::R_PI_2 * real::R_PI_2, myfun, myerr, 15); return EXIT_SUCCESS; }