What are the most accurate algorithms for evaluating a polynomial using floating point arithmetic? The internet seems to suggest that Horner's method is commonly used. In particular I have a cubic polynomial of the form f(x) = a x^3 + b x^2 + c x + d, which I am evaluating in the naive way with floating point arithmetic, and was curious whether some recursive approach or some refactoring would give a more accurate result. That is, a result closer to what I would get if I were to evaluate the polynomial exactly or with a metric ton of precision in Mathematica.
$\begingroup$
$\endgroup$
2
-
$\begingroup$ A key question, also, is: what is your input? Sometimes you get coefficients in different "bases", for instance you may want to evaluate the polynomial such that $p(0)=0$, $p(1)=5$, $p(2)=3$, $p(3)=5$ in $x=10$ (interpolation). $\endgroup$– Federico PoloniCommented Sep 1, 2020 at 14:52
-
$\begingroup$ Are you dealing with piecewise cubic splines? $\endgroup$– Brian BorchersCommented Sep 2, 2020 at 15:48
Add a comment
|
3 Answers
$\begingroup$
$\endgroup$
1
Horner is indeed the most stable way to evaluate a polynomial (and you get the bonus of evaluating its derivatives with not too much extra cost). Higham presents a nice error analysis of the algorithm (Accuracy and stability of numerical algorithms, 2nd edition, p.94). He also presents an algorithm that includes a running error bound so you have an idea on the difference between the real value and what Horner calculated (Algorithm 5.1 on p.95).
If your case is limited to cubic polynomials and you're worried about loop overhead, you could unroll the loops by hand. But I doubt you will gain much. Stick to the tried and tested.
-
3$\begingroup$ The accuracy of the Horner scheme can be improved by using compensated arithmetic. See: Philippe Langlois and Nicolas Louvet. "How to ensure a faithful polynomial evaluation with the compensated Horner algorithm." In 18th IEEE Symposium on Computer Arithmetic (ARITH'07), pp. 141-149. IEEE, 2007 (draft online); Stef Graillat, Philippe Langlois, and Nicolas Louvet. "Algorithms for Accurate, Validated and Fast Polynomial Evaluation", Japan J. Indust. Appl. Math., 26 (2009), 191–214 (online) $\endgroup$– njuffaCommented Jun 19, 2016 at 18:46
$\begingroup$
$\endgroup$
The stability of Horner can be improved by subtracting a constant from each $x$ in the Horner form.
$\begingroup$
$\endgroup$
As an alternative you may use the barycentric Lagrange basis which has the same evaluation complexity as the Newton basis. The barycentric formulation is also used in the MATLAB Chebfun implementation. For details see papers of L. N. Trefethen.
