Knuth–Eve algorithm: Difference between revisions

Revision as of 06:36, 26 July 2025

Revision as of 05:06, 27 July 2025

Line 2:

<!-- Important, do not remove anything above this line before article has been created. -->

The '''Knuth-Eve algorithm''' is an [[algorithm]] for [[polynomial evaluation]]. Specifically, it can evaluate an arbitrary [[polynomial]] with [[real number|real]] coefficients of one real variable. It [[Polynomial_evaluation#Evaluation_with_preprocessing|preprocesses]] the coefficients of the polynomial to reduce the number of multiplications required at [[Execution_(computing)#Runtime|runtime]].

The key ideas used in this algorithm were originally proposed by [[Donald Knuth]]. His procedure opportunistically exploits structure in the polynomial being evaluated.<ref name="knuth1962"/> Eve determined for which polynomials this structure exists, and they gave a simple method of [[#"Preconditioning"|"preconditioning"]] polynomials to endow them with that structure.<ref name="eve1964"/>

Line 12:

Consider an arbitrary polynomial <math>p \in \mathbb{R}[x]</math> of [[Degree of a polynomial|degree]] <math>n</math>. Assume that <math>n \geq 3</math>. Define <math>m</math> such that: if <math>n</math> is odd then <math>n = 2m+1</math>, and if <math>n</math> is even then <math>n = 2m+2</math>.

Unless otherwise stated, all variables represent either [[Real number|real numbers]] or univariate [[Polynomial|polynomials]] with real coefficients. All operations are done over <math>\mathbb{R}</math>.

Again, the goal is to create an algorithm that returns <math>p(x)</math> given any <math>x</math>. The algorithm is allowed to depend on the polynomial <math>p</math> itself.

Line 24:

<math display="block"> p(x) = q(x) \cdot (x^2 - \alpha) + (\beta x + \gamma), </math>

where <math>x^2 - \alpha</math> is the divisor. Picking a value for <math>\alpha</math> fixes both the quotient <math>q</math> and the coefficients in the remainder <math>\beta</math> and <math>\gamma</math>. Then, theThe key idea is to cleverly choose <math>\alpha</math> such that <math>\beta = 0</math>, so that

<math display="block"> p(x) = q(x) \cdot (x^2 - \alpha) + \gamma. </math>

Finally, weWe apply this procedure [[Recursion (computer science)|recursively]] to <math>q</math>, expressing

<math display="block"> p(x) = \left( \left( q_mq(x) \cdot (x^2 - \alpha_m) + \gamma_m \right) \cdots \right) \cdot (x^2 - \alpha_1) + \gamma_1. </math>

After <math>m</math> recursive calls, the quotient <math>q_mq</math> is either a [[Linear function (calculus)|linear]] or a [[Quadratic function|quadratic]] polynomial. In this base case, the polynomial can be evaluated with (say) [[Horner's method]].<ref name="knuth1962"/>

==== "Preconditioning" ====

For arbitrary <math>p</math>, it may not be possible to force <math>\beta = 0</math> at every step of the recursion.<ref name="knuth1962"/> Still, considerConsider the polynomials <math>p^e</math> and <math>p^o</math> with coefficients taken from the even and odd terms of <math>p</math> respectively, so that

<math display="block"> p(x) = p^e(x^2) + x \cdot p^o(x^2). </math>

Line 42:

If every [[Zero of a function#Polynomial roots|root]] of <math>p^o</math> is real, then it is possible to write <math>p</math> in the form given above. Each <math>\alpha_i</math> is a different root of <math>p^o</math>, counting [[Multiplicity (mathematics)#Multiplicity of a root of a polynomial|multiple roots]] as distinct.<ref name="overill1997"/> Furthermore, if all the roots of <math>p</math> (except perhaps one) lie in one half of the [[Complex number|complex plane]], then every root of <math>p^o</math> is real.<ref name="eve1964"/>

Ultimately, it may be necessary to "precondition" <math>p</math> by shifting it {{--}} by setting <math>p(x) \gets p(x + t)</math> for some <math>t</math> {{--}} to endow it with the structure that all its roots lie in one half of the complex plane. At runtime, this shift has to be "undone" by first setting <math>x \gets x - t</math>.

=== Preprocessing step ===

The following algorithm is run once for a fixedgiven polynomial <math>p</math>. Its results are used by the [[#Evaluation step|evaluation step]], and they can be reused across many calls to that step.

At this point, the values of <math>x</math> that <math>p</math> will be evaluated on are not known.

Line 53:

{{framebox|blue}}

* Let <math>r_1, \cdots, r_n \in \mathbb{C}</math> be the complex roots of <math>p</math>

* Let <math display="inline">t = \max_{i = 1}^n \text{Re}(r_i)</math><hr/>

* LetSet <math>q_0(x)p =\gets p(x + t)</math>

<hr/>

* Let <math>q_0^e</math> and <math>q_0^o</math> be the polynomials such that <math>q_0(x) = q_0^e(x^2) + x \cdot q_0^o(x^2)</math>

* Let <math>\alpha_1, \cdots \alpha_m \in \mathbb{R}p^e</math> be all the roots ofand <math>q_0p^o</math>. Allbe ofthe itspolynomials rootssuch willthat be<math>p(x) real.= p^e(x^2) + x \cdot p^o(x^2)<hr/math>

* Let <math>\alpha_1, \cdots \alpha_m \in \mathbb{R}</math> be all the roots of <math>p^o</math>. All of its roots will be real.

<hr/>

* Initialize <math>q \gets p</math>

* For <math>i \gets 1, \cdots, m</math>:

** Divide <math>q_{i-1}q</math> by <math>x^2 - \alpha_i</math> to get quotient <math>q_iq^\prime \in \mathbb{R}[x]</math> and remainder <math>\gamma_i \in \mathbb{R}</math>. The remainder will be a constant polynomial {{--}} a number.

** Set <math>q \gets q^\prime</math>
<hr/>
* ''Output:'' The derived values <math>t</math>, <math>\alpha_1, \cdots, \alpha_m</math>, and <math>\gamma_1, \cdots, \gamma_m</math>; as well as the base-case polynomial <math>q_mq</math>

{{frame-footer}}

</div>

Line 65 ⟶ 70:

=== Evaluation step ===

The following algorithm evaluates <math>p</math> at some, now known, point <math>x</math>. It consumes the output of the [[#Preprocessing step|preprocessing step]].

<div style="margin-left: 35px;">

{{framebox|blue}}

* LetSet <math>ux =\gets x - t</math>

* Let <math>s = ux^2</math>. Compute this once so it can be reused throughout the algorithm.
<hr/>

* Compute <math>y_my =\gets q_mq(ux)</math> using [[Horner's method]]

* For <math>i \gets m, \cdots, 2, 1</math>:

** Let <math>y_{i-1}y =\gets y_iy \cdot (s - \alpha_i) + \gamma_i</math>

* ''Output:'' <math>y_0y</math>

{{frame-footer}}

</div>