Sollya User's Manual

Name:

remez computes the minimax of a function on an interval.

Library names:

sollya_obj_t sollya_lib_remez(sollya_obj_t, sollya_obj_t, sollya_obj_t, ...) sollya_obj_t sollya_lib_v_remez(sollya_obj_t, sollya_obj_t, sollya_obj_t, va_list)

Usage:

remez(f, n, range, w, quality, bounds) : (function, integer, range, function, constant, range) -> function remez(f, L, range, w, quality, bounds) : (function, list, range, function, constant, range) -> function

Parameters:

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].

Description:

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.
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.
Examples of use:
[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).

Example 1:

   > p = remez(exp(x),5,[0;1]);
   > degree(p);
   5
   > dirtyinfnorm(p-exp(x),[0;1]);
   1.1295698151096148707171193829266077607222634589363e-6

Example 2:

   > p = remez(1,[|0,2,4,6,8|],[0,Pi/4],1/cos(x));
   > canonical=on!;
   > p;
   0.99999999994393732180959690352543887130348096061124 + -0.49999999571556857768772053063721544670949467222259 * x^2 + 4.1666613233473633009941059480570275870113220089059e-2 * x^4 + -1.3886529147145693651355523880319714051047635695061e-3 * x^6 + 2.4372679177224179934800328511009205218114284220126e-5 * x^8

Example 3:

   > p1 = remez(exp(x),5,[0;1],default,1e-5);
   > 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

Example 4:

   > L = [|exp(x), sin(x), cos(x)-1, sin(x^3)|];
   > 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

Example 5:

   > f = sin(x);
   > I = [-3b-5;-1b-1074];
   > time(popt = remez(1, [|1, 3, 4, 5, 7, 8, 9|], I, 1/f));
   0.116157125999999999999999999999999999997386771862303
   > time(p1 = remez(1, [|1, 3, 4, 5, 7, 8, 9|], I, 1/f, default, [0, 1b-73]));
   8.6851247000000000000000000000000000000254830624298e-2
   > time(p2 = remez(1, [|1, 3, 4, 5, 7, 8, 9|], I, 1/f, default, [3b-72, +@Inf@]));
   9.6260799999999999999999999999999999996109922438892e-2
   > 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

See also: dirtyinfnorm, infnorm, fpminimax, guessdegree, taylorform, taylor

Go back to the list of commands