# Changeset 1014Tweet

Ignore:
Timestamp:
Oct 10, 2011, 3:33:15 AM (9 years ago)
Message:

test: various improvements to the Remez exchange solver.

Location:
trunk/test/math
Files:
2 edited

Unmodified
Added
Removed
• ## trunk/test/math/remez-solver.h

 r1010 } 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; for (int n = 0; n < steps; n++) { FindError(); FindExtrema(); Step(); } FindError(); FindExtrema(); Step(); void FindZeroes() { 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); while (fabs(a - b) > (real)1e-140) /* 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 { real 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; while (fabs(left.value - right.value) > limit) { real c = (a + b) * (real)0.5; real ec = ChebyEval(c) - Value(c); if ((ea < (real)0 && ec < (real)0) || (ea > (real)0 && ec > (real)0)) { a = c; ea = ec; } 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)) left = mid; else { b = c; eb = ec; } right = mid; } zeroes[i] = a; } } void FindError() { zeroes[i] = mid.value; } } 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 (;;) { 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; } } 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; } } 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])); } } 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"); } } real Error(real const &x) { return m_error(x * m_k2 + m_k1); real Weight(real const &x) { return m_weight(x * m_k2 + m_k1); } private: RealFunc *m_func, *m_error; RealFunc *m_func, *m_weight; real m_k1, m_k2, m_invk1, m_invk2; };
• ## trunk/test/math/remez.cpp

 r1012 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;
Note: See TracChangeset for help on using the changeset viewer.