- f is the function to be approximated
- n is the degree of the polynomial that must approximate f
- L is a list of integers or a list of functions and indicates the basis for the approximation of f
- range is the interval where the function must be approximated
- w (optional) is a weight function. Default is 1.
- quality (optional) is a parameter that controls the quality of the returned polynomial p, with respect to the exact minimax p*. Default is 1e-5.
- bounds (optional) is a parameter that allows the user to make the algorithm stop earlier, whenever a given accuracy is reached or a given accuracy is proved unreachable. Default is [0, +Inf].

- remez computes an approximation of the function f with respect to the weight function w on the interval range. More precisely, it searches p such that ||p*w-f|| is (almost) minimal among all p of a certain form. The norm is the infinity norm, e.g. ||g|| = max {|g(x)|, x in range}.
- If w=1 (the default case), it consists in searching the best polynomial approximation of f with respect to the absolute error. If f=1 and w is of the form 1/g, it consists in searching the best polynomial approximation of g with respect to the relative error.
- If n is given, p is searched among the polynomials with degree not greater than n. If L is given and is a list of integers, p is searched as a linear combination of monomials X^k where k belongs to L. In the case when L is a list of integers, it may contain ellipses but cannot be end-elliptic. If L is given and is a list of functions g_k, p is searched as a linear combination of the g_k. In that case L cannot contain ellipses. It is the user responsibility to check that the g_k are linearly independent over the interval range. Moreover, the functions w*g_k must be at least twice differentiable over range. If these conditions are not fulfilled, the algorithm might fail or even silently return a result as if it successfully found the minimax, though the returned p is not optimal.
- The polynomial is obtained by a convergent iteration called Remez'
algorithm (and an extension of this algorithm, due to Stiefel).
The algorithm computes a sequence p1, ..., pk, ...
such that ek = ||pk*w-f|| converges towards
the optimal value e. The algorithm is stopped when the relative error
between ek and e is less than quality.

For this reason, the returned polynomial usually has dyadic coefficients. However, one should not rely on that property. Especially, a noticeable exception to that statement is when we manage to detect that f/w exactly simplifies to a polynomial of the required degree, in which case we return it as it is without further rounding. - The optional argument bounds is an interval [satisfying_err, target_err]
with the following behavior:
- if, during the algorithm, we manage to prove that target_err is unreachable, we stop the algorithm returning the last computed polynomial.
- if, during the algorithm, we obtain a polynomial with an error smaller than satisfying_err, we stop the algorithm returning that polynomial.
- otherwise we loop until we find an optimal polynomial with the required quality, as usual.

[0, +Inf] (compute the optimal polynomial with the required quality)

[target_err] (stops as soon as a polynomial achieving target_err is obtained or as soon as such a polynomial is proved not to exist).

[0, target_err] (finds the optimal polynomial, but provided that its error is smaller than target_err).

[satisfying_err, +Inf] (stops as soon as a polynomial achieving satisfying_err is obtained. If such a polynomial does not exist, returns the optimal polynomial).

> degree(p);

5

> dirtyinfnorm(p-exp(x),[0;1]);

1.1295698151096148707171193829266077607222634589363e-6

> canonical=on!;

> p;

0.99999999994393732180959690352543887130348096061124 - 0.49999999571556857768772053063721544670949467222259 * x^2 + 4.1666613233473633009941059480570275870113220089059e-2 * x^4 - 1.3886529147145693651355523880319714051047635695061e-3 * x^6 + 2.4372679177224179934800328511009205218114284220126e-5 * x^8

> p2 = remez(exp(x),5,[0;1],default,1e-10);

> p3 = remez(exp(x),5,[0;1],default,1e-15);

> dirtyinfnorm(p1-exp(x),[0;1]);

1.1295698151096148707171193829266077607222634589363e-6

> dirtyinfnorm(p2-exp(x),[0;1]);

1.12956980227478675612619255125474525171079325793124e-6

> dirtyinfnorm(p3-exp(x),[0;1]);

1.12956980227478675612619255125474525171079325793124e-6

> g = (2^x-1)/x;

> p1 = remez(g, L, [-1/16;1/16]);

> p2 = remez(g, 3, [-1/16;1/16]);

> dirtyinfnorm(p1 - g, [-1/16;1/16]);

9.8841323829271038137685646777951687620288462194746e-8

> dirtyinfnorm(p2 - g, [-1/16;1/16]);

2.54337800593461418356437401152248866818783932027105e-9

> I = [-3b-5;-1b-1074];

> time(popt = remez(1, [|1, 3, 4, 5, 7, 8, 9|], I, 1/f));

0.185799695000000000000000000000000000002291314966024

> time(p1 = remez(1, [|1, 3, 4, 5, 7, 8, 9|], I, 1/f, default, [0, 1b-73]));

0.134645933000000000000000000000000000002249871361033

> time(p2 = remez(1, [|1, 3, 4, 5, 7, 8, 9|], I, 1/f, default, [3b-72, +@Inf@]));

0.150942799000000000000000000000000000010601473971874

> dirtyinfnorm(popt/f-1, I);

2.06750931454112835098093903810531156576504665659064e-22

> dirtyinfnorm(p1/f-1, I);

2.49711266837493110470637913808914046704452778960875e-22

> dirtyinfnorm(p2/f-1, I);

5.4567247553615435246376977231253834265248756996947e-22

> 1b-73;

1.05879118406787542383540312584955245256423950195312e-22

> 3b-72;

6.3527471044072525430124187550973147153854370117187e-22