I have implemented the Remez algorithm for Chebyshev rational approximation using a multiple precision numerics library (mpmath via SymPy). Wikipedia has suggested using $\max\{|z_i|\} - \min\{|z_i|\}$ as a stopping criterion. Here $z_i$ are the values of the error function $E(x) = f(x) - r_i(x)$ evaluated at the local extrema, that is, $z_i = E(x_i)$. Other sources have suggested stopping when the alternates $x_i$ stop changing, but Wikipedia's suggestion seems more reasonable to me, as the theory says that the perfect minimax rational approximation has equioscillating error, and $\max\{|z_i|\} - \min\{|z_i|\}$ is a measure of closeness to that.
Empirically, when I perform the computations with $N$ digits of precision, the convergence happens when $\max\{|z_i|\} - \min\{|z_i|\}$ goes under $10^N$ (that is, iterating past that doesn't produce any further improvement), so I've set that as my stopping condition. If anyone knows why this might be the case (or might not actually be the case) that would be great to know.
My question is, once the approximation has converged using $N$ digits of precision and $\max\{|z_i|\} - \min\{|z_i|\} < \epsilon$, how many digits of the coefficients in final approximation $r(x)$ are correct? Obviously my end goal is to get at least enough digits for some double floating point calculations.
My only goal post thus far has been a paper with published coefficients.
My question is twofold:
1. How can I tell how many digits of the coefficients in my computed approximation are accurate?
2. How can I put a lower bound on the number of accurate digits (preferably by increasing the number of digits of precision in the algorithm, $N$, or by changing the convergence conditions)?
I only care about approximating $e^{-x}$ on $[0, \infty)$, although general answers would be nice. I've been using the technique of the above paper (translate to $[-1, 1]$ via $z = c\frac{1 + x}{1 - x}$ with suitably chosen $c$ to spread out the errors, then translate back with $x = \frac{z - c}{z + c}$).
So far, I've thought of the following:
Suppose we have a polynomial approximation instead of a rational approximation (i.e., $E(x) = f(x) - p(x) < \epsilon$ is the error for the best minimax approximation). If we perturb a coefficient of $p(x)$ by $\delta$, and let $$E^* = f(x) - p^*(x),$$ where $$p^*(x) = a_0 + \cdots + (a_k + \delta)x^k + \cdots + a_nx^n,$$ then $$|E^*| = |E - \delta x^k| <= \epsilon + \delta$$ because $x\in [-1, 1].$ But I'm not sure if this helps me. I also don't know how to apply it to the rational case (how does the denominator affect things?), and finally, I don't know how the inverse transformation $\frac{z - c}{z + c}$ will affect things.
fpminimax
function of the Sollya tool, the second part an open research problem. $\endgroup$ – njuffa Dec 6 '16 at 6:17