diff options
author | S. Gilles <sgilles@math.umd.edu> | 2019-05-11 07:55:37 -0400 |
---|---|---|
committer | S. Gilles <sgilles@math.umd.edu> | 2019-05-11 07:55:37 -0400 |
commit | a00e33e21435b8a4961d6e0f910e09f9c17f65ba (patch) | |
tree | 6d9a9b3bffbe8d43f72b7645d5ca1a63c8d767aa | |
parent | 9dd2a40974da19018eeb33338cfcc486f4eeca8d (diff) | |
download | mc-a00e33e21435b8a4961d6e0f910e09f9c17f65ba.tar.gz |
Fix Remez algorithm.
-rw-r--r-- | lib/math/ancillary/generate-minimax-by-Remez.gp | 25 |
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.gp index 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"); |