summaryrefslogtreecommitdiff log msg author committer range
diff options
 context: 12345678910152025303540 space: includeignore mode: unifiedssdiffstat only
author committer S. Gilles 2019-05-11 07:55:37 -0400 S. Gilles 2019-05-11 07:55:37 -0400 a00e33e21435b8a4961d6e0f910e09f9c17f65ba (patch) 6d9a9b3bffbe8d43f72b7645d5ca1a63c8d767aa 9dd2a40974da19018eeb33338cfcc486f4eeca8d (diff) mc-a00e33e21435b8a4961d6e0f910e09f9c17f65ba.tar.gz
Fix Remez algorithm.
-rw-r--r--lib/math/ancillary/generate-minimax-by-Remez.gp25
1 files changed, 13 insertions, 12 deletions
 diff --git a/lib/math/ancillary/generate-minimax-by-Remez.gp b/lib/math/ancillary/generate-minimax-by-Remez.gpindex a388f4e..2863120 100644--- a/lib/math/ancillary/generate-minimax-by-Remez.gp+++ b/lib/math/ancillary/generate-minimax-by-Remez.gp@@ -1,8 +1,6 @@ /* Attempts to find a minimax polynomial of degree n for f via Remez- algorithm. The Remez algorithm appears to be slightly shot, but- the initial step of approximating by Chebyshev nodes works, and- is usually good enough.+ algorithm. */ { chebyshev_node(a, b, k, n) = my(p, m, c);@@ -62,19 +60,21 @@ xnext = a + k*(b - a)/gran; ynext = err(f, n, v, xnext); - if(ycur > yprev && ycur > ynext, listput(l, [xcur, abs(ycur)]),);+ if(ycur > yprev && ycur > ynext, listput(l, [xcur, ycur]),); ); vecsort(l, 2); if(length(l) < n + 2 || l[1][2] < 2^(-2000), return("q"),); len = length(l); - lnext = listcreate();+ X = vector(n + 2);+ for(j = 1, n + 2,+ lnext = listcreate(); thisa = l[j][1] - (b-a)/gran; thisb = l[j][1] + (b-a)/gran; - xprev = thisa - (thisb - a)/gran;+ xprev = thisa - (thisb - thisa)/gran; yprev = err(f, n, v, xprev); xcur = thisa;@@ -93,14 +93,15 @@ xnext = thisa + k*(thisb - thisa)/gran; ynext = abs(f(xnext) - p_eval(n, v, xnext)); - if(ycur > yprev && ycur > ynext, listput(lnext, xcur),);+ if(ycur > yprev && ycur > ynext, listput(lnext, [xcur, ycur]),); );++ vecsort(lnext, 2);+ if(length(lnext) < 1, return("q"),);+ X[j] = lnext[1][1];+ listkill(lnext); );- vecsort(lnext, cmp); listkill(l);- X = vector(n + 2);- for (k = 1, min(n + 2, length(lnext)), X[k] = lnext[k]);- listkill(lnext); vecsort(X); return(X); }@@ -158,7 +159,7 @@ print("\n"); print("Minmimaxing x*cot(x), degree 8, on [-Pi/(4 * 256), Pi/(4 * 256)]:"); find_minimax(cotx, 8, -Pi/1024, Pi/1024) print("\n");-print("(Take the first v, and remember to divide by x)\n");+print("(Remember to divide by x)\n"); print("\n"); print("---\n"); print("\n");