summaryrefslogtreecommitdiff
path: root/lib
diff options
context:
space:
mode:
authorS. Gilles <sgilles@math.umd.edu>2019-05-11 07:55:37 -0400
committerS. Gilles <sgilles@math.umd.edu>2019-05-11 07:55:37 -0400
commita00e33e21435b8a4961d6e0f910e09f9c17f65ba (patch)
tree6d9a9b3bffbe8d43f72b7645d5ca1a63c8d767aa /lib
parent9dd2a40974da19018eeb33338cfcc486f4eeca8d (diff)
downloadmc-a00e33e21435b8a4961d6e0f910e09f9c17f65ba.tar.gz
Fix Remez algorithm.
Diffstat (limited to 'lib')
-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.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");