ref: bfc9204963426bf0e8121b416458550b8f2616ad
dir: /lib/math/ancillary/generate-minimax-by-Remez.gp/
/*
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.
*/
{ chebyshev_node(a, b, k, n) = my(p, m, c);
p = (b + a)/2;
m = (b - a)/2;
c = cos(Pi * (2*k - 1)/(2*n));
return(p + m*c);
}
{ remez_step(f, n, a, b, x) = my(M, xx, bvec, k);
M = matrix(n + 2, n + 2);
bvec = vector(n + 2);
for (k = 1, n + 2,
xx = x[k];
for (j = 1, n + 1,
M[k,j] = xx^(j - 1);
);
M[k, n + 2] = (-1)^k;
bvec[k] = f(xx);
);
return(mattranspose(M^-1 * mattranspose(bvec)));
}
{ p_eval(n, v, x) = my(s, k);
s = 0;
for (k = 1, n + 1,
s = s + v[k]*x^(k - 1)
);
return(s);
}
{ err(f, n, v, x) =
return(abs(f(x) - p_eval(n, v, x)));
}
{ find_M(f, n, v, a, b, depth) = my(X, gran, l, lnext, len, xprev, xcur, xnext, yprev, ycur, ynext, thisa, thisb, k, j);
gran = 10000 * depth;
l = listcreate();
xprev = a - (b - a)/gran;
yprev = err(f, n, v, xprev);
xcur = a;
ycur = err(f, n, v, xprev);
xnext = a + (b - a)/gran;
ynext = err(f, n, v, xprev);
for (k = 2, gran,
xprev = xcur;
yprev = ycur;
xcur = xnext;
ycur = ynext;
xnext = a + k*(b - a)/gran;
ynext = err(f, n, v, xnext);
if(ycur > yprev && ycur > ynext, listput(l, [xcur, abs(ycur)]),);
);
vecsort(l, 2);
if(length(l) < n + 2 || l[1][2] < 2^(-2000), return("q"),);
len = length(l);
lnext = listcreate();
for(j = 1, n + 2,
thisa = l[j][1] - (b-a)/gran;
thisb = l[j][1] + (b-a)/gran;
xprev = thisa - (thisb - a)/gran;
yprev = err(f, n, v, xprev);
xcur = thisa;
ycur = err(f, n, v, xprev);
xnext = thisa + (thisb - thisa)/gran;
ynext = err(f, n, v, xprev);
for (k = 2, gran,
xprev = xcur;
yprev = ycur;
xcur = xnext;
ycur = ynext;
xnext = thisa + k*(thisb - thisa)/gran;
ynext = abs(f(xnext) - p_eval(n, v, xnext));
if(ycur > yprev && ycur > ynext, listput(lnext, xcur),);
);
);
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);
}
{ find_minimax(f, n, a, b) = my(c, k, j);
c = vector(n + 2);
for (k = 1, n + 2,
c[k] = chebyshev_node(a, b, k, n + 2);
);
for(j = 1, 100,
v = remez_step(f, n, a, b, c);
print("v = ", v);
c = find_M(f, n, v, a, b, j);
if(c == "q", return,);
print("c = ", c);
);
}
{ sinoverx(x) =
return(if(x == 0, 1, sin(x)/x));
}
{ tanoverx(x) =
return(if(x == 0, 1, tan(x)/x));
}
{ atanxoverx(x) =
return(if(x == 0, 1, atan(x)/x));
}
{ cotx(x) =
return(1/tanoverx(x));
}
print("\n");
print("Minimaxing sin(x) / x, degree 6, on [-Pi/(4 * 256), Pi/(4 * 256)]:");
find_minimax(sinoverx, 6, -Pi/1024, Pi/1024)
print("\n");
print("(You'll need to add a 0x0 at the beginning to make a degree 7...\n");
print("\n");
print("---\n");
print("\n");
print("Minimaxing cos(x), degree 7, on [-Pi/(4 * 256), Pi/(4 * 256)]:");
find_minimax(cos, 7, -Pi/1024, Pi/1024)
print("\n");
print("---\n");
print("\n");
print("Minmimaxing tan(x) / x, degree 6, on [-Pi/(4 * 256), Pi/(4 * 256)]:");
find_minimax(tanoverx, 6, -Pi/1024, Pi/1024)
print("\n");
print("(You'll need to add a 0x0 at the beginning to make a degree 7...\n");
print("\n");
print("---\n");
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("\n");
print("---\n");
print("\n");
print("Minmimaxing tan(x) / x, degree 10, on [0, 15.5/256]:");
find_minimax(tanoverx, 10, 0, 15.5/256)
print("\n");
print("(You'll need to add a 0x0 at the beginning to make a degree 11...\n");
print("\n");
print("---\n");
print("\n");
print("Minmimaxing atan(x) / x, degree 12, on [0, 15.5/256]:");
find_minimax(atanxoverx, 12, 0, 1/16)
print("\n");
print("(You'll need to add a 0x0 at the beginning to make a degree 13...\n");
print("\n");
print("---\n");
print("Remember that there's that extra, ugly E term at the end of the vector that you want to lop off.\n");