Fast evaluation of polynomials

Question

Hello everybody !

I was reading a book on geometry which taught me that one could compute the volume of a simplex through the determinant of a matrix, and I thought (I'm becoming a worse computer scientist each day) that if the result is exact this may not be the computationally fastest way possible to do it.

Hence, the following problem : if you are given a polynomial in one (or many) variables $\alpha_1 x^1 + \dots + \alpha_n x^n$, what is the cheapest way (in terms of operations) to evaluate it ?

Indeed, if you know that your polynomial is $(x-1)^{1024}$, you can do much, much better than computing all the different powers of $x$ and multiply them by their corresponding factor.

However, this is not a problem of factorization, as knowing that the polynomial is equal to $(x-1)^{1024} + (x-2)^{1023}$ is also much better than the naive evaluation.

Of course, multiplication and addition all have different costs on a computer, but I would be quite glad to understand how to minimize the "total number of operations" (additions + multiplications) for a start ! I had no idea how to look for the corresponding litterature, and so I am asking for your help on this one :-)

Thank you !

Nathann

P.S. : I am actually looking for a way, given a polynomial, to obtain a sequence of addition/multiplication that would be optimal to evaluate it. This sequence would of course only work for THIS polynomial and no other. It may involve working for hours to find out the optimal sequence corresponding to this polynomial, so that it may be evaluated many times cheaply later on.

Answer 1

If the polynomial is given as $\alpha_0x^0+\dots+\alpha_nx^n$ and you do not know a priori anything about the $\alpha_i$’s, then you can’t do better than Horner’s scheme (which takes $n$ additions and multiplications). If you know that the polynomial is sparse and you are given a list of nonzero coefficients, you can evaluate the individual terms using repeated squaring (this takes about $k$ additions and $O(k\log n)$ multiplications, where $k$ is the number of nonzero terms). Other information about the polynomial may also help in principle, such as some sort of symmetries in the coefficient list.

Answer 2

The simplest version of this question is: what is the quickest way to evaluate $x^n?$ For $n = 2^k,$ $k$ repeated squarings is obviously best, but for more complicated $n$ I believe that finding the optimum is very hard -- see Knuth, vol 2 for (much) more on these so-called "multiplication trees".

Answer 3

EDIT: Looks like I overlooked that the OP stipulated he wants to minimize the total number of additions and multiplications. (Although he said he wanted to do that "to start," so arguably the below is still relevant.)

However, to address the question as stated, what you are essentially looking for is the https://en.wikipedia.org/wiki/Arithmetic_circuit_complexity

---

Consider the polynomial $f(x) = nx$, where $n$ is an integer. Here are two algorithms which will evaluate this polynomial:

Algorithm 1. Multiply $n$ by $x$.

Algorithm 2. Calculate $x + x + \ldots + x$.

Which is more efficient? Given fixed $n$, this depends on your processor architecture. And this is just about the simplest case imaginable -- we only have one variable, the polynomial is linear, and we're not even thinking about pipelined calculations yet. Also, as mentioned before, you are going to have to formalize the problem in some way which eliminates the "algorithm" consisting of a table giving the value at each machine-sized number. As stated, I don't think the question is answerable.

Answer 4

If we are on an architecture which has multiple cores, CPU pipelines and multi-media extensions (MME) then Horner's method really doesn't have to be best.

If the polynomial is large, you can split into as many bins as you have processor cores. If you have two cores, for example $$f(x)=\sum_k\alpha_kx^k=x\left(\sum\alpha_{2k+1} (x^2)^k\right) + \sum\alpha_{2k}(x^2)^k$$ Each sum evaluated in a separate thread. Then of course you can apply any optimization by treating as a separate polynomial if you want.

---

Another idea is to precalculate $x,x^2,x^4,\cdots$ and store them in a table, then you can evaluate say $x^{43}$ with help of $43 = 32+8+2+1$, so you do $x \cdot x^2 \cdot x^8 \cdot x^{32}$ by grabbing those four numbers from the table. That's 4 multiplications instead of 43 - in general will be logarithmic. Worst case would be 2-logarithm rounded up number of multiplications. This can be useful as parallellized MME-instructions are common on modern CPUs.

Answer 5

The cheapest way of finding the value of a polynomial, given unlimited preprocessing resources, is to look up the precalculated value in the table. However, if you know you are going to need several more values evaluated at successive intervals, you might try a method similar to that desired by Charles Babbage: differences. Namely, store the value and the the n kth order differences (similar to evaluations at derivatives) for point x, and then use n additions to derive the differences and value for the polynomial at the point x+1. If you need to loop through to evaluate the polynomial at successive integers, this gets those values with O(n) additions per evaluation point.

(Of course needing random or real access to the polynomial will require something different, but you might find storing values at derivatives useful for evaluating the polynomial at near by points, especially if multiplication is expensive..)

Gerhard "Email Me About System Design" Paseman, 2011.07.04

Answer 6

If you want to evaluate the polynomial at a lot of equidistant points, you can do "forward differencing"; here are 3 slides explaining the method: http://zach.in.tu-clausthal.de/teaching/info2_11/folien/evaluating%20a%20polynomial%20at%20equidistant%20points.pdf (they are in German, but I believe you'll still get it).

Answer 7

For evaluation of a general polynomial in one variable, the provably fastest method is Horner's scheme as Emil has pointed out. It is worth mentioning that this scheme has a more popular face in the form of little Bézout's theorem paired with Synthetic division as is often taught in Precalculus courses in the U.S.A. This is implicitely present in the Wikipedia article for Horner's method, but the relationship is not well explained. A convenient consequence of this fact is that on computer algebra systems, storing polynomials in Horner normal form adds no performance benefit when it comes to evaluation.
To see that these two algorithms are the same one only needs to verify that the base step and the recursive procedure for both match. I do this for evaluation at $x_0$ of $$f(x)=\sum_{i=0}^n a_ix^i = \left((\dots(a_nx+a_{n-1})x+\dots)x+a_1\right)x+a_0,$$ where the center expression is the general form of a polynomial in one variable and the right-hand-side (RHS) is its Horner normal form. In synthetic division as in evaluation of the RHS following the order of operations, one begins by multipling $x_0\cdot a_n$ which will be denoted $y_n$. For the recursive step in synthetic division one sets $y_{i-1}=x_0\cdot y_i+a_{i}$ for $n>i\geq0$. For $i>0$, this is precisely the calculation of the $i^{th}$ set of parentheses counting out-to-in and thus can be thought of as the $(n-i)^{th}$ iteration of evaluating the RHS of the euqation.
In this example $f(x_0)=y_0$.

Answer 8

If I understood your question correctly, you are willing to do an arbitrarily large amount of precomputing on your polynomial in order to make the evaluation process at run time as fast as possible.

In other words, you want to find a some way of making your polynomial evaluation "sparse" in some sense.

My first intuition is that this is not possible for an arbitrary polynomial, i.e. the set of polynomials for which a "sparsification" is possible is a set of measure zero. Which is not to say it is not possible for certain specific polynomials.

In general, "sparseness" usually indicates some sort of underlying mathematical structure, which suggests you should attempt to understand said structure first.

Otherwise I believe the problem to be NP-complete for the general case. I will edit this answer with a proof if I can think of one.

Answer 9

it's probably easier to see just see it in action instead of tryign to explain the symbols and notation : here, I've automated a process that incorporates the Homer nested formatting with the binary-squaring algorithm.

The generated algebraic statement could be directly plugged into bc or Wolfram Alpha to confirm the result. Modify those ^ to ** if you wanna plug them into python or perl.

---

• all the "1" bits reside on the left half of the formula, and
• all the squarings reside on the right

I preferred using LSB bit-string instead of MSB in this structure because it's a lot easier to read -

---

  17^1
| expn LSB bitstring := 1
| 17

  17^2
| expn LSB bitstring := 01
| (17)^2

  17^3
| expn LSB bitstring := 11
| 17*(17)^2

  17^5
| expn LSB bitstring := 101
| 17*((17)^2)^2

  17^7
| expn LSB bitstring := 111
| 17*(17*(17)^2)^2

  17^11
| expn LSB bitstring := 1101
| 17*(17*((17)^2)^2)^2

  17^13
| expn LSB bitstring := 1011
| 17*((17*(17)^2)^2)^2

  17^17
| expn LSB bitstring := 10001
| 17*((((17)^2)^2)^2)^2

  17^19
| expn LSB bitstring := 11001
| 17*(17*(((17)^2)^2)^2)^2

  17^23
| expn LSB bitstring := 11101
| 17*(17*(17*((17)^2)^2)^2)^2

  17^29
| expn LSB bitstring := 10111
| 17*((17*(17*(17)^2)^2)^2)^2

  17^31
| expn LSB bitstring := 11111
| 17*(17*(17*(17*(17)^2)^2)^2)^2

  17^37
| expn LSB bitstring := 101001
| 17*((17*(((17)^2)^2)^2)^2)^2

  17^41
| expn LSB bitstring := 100101
| 17*(((17*((17)^2)^2)^2)^2)^2

  17^43
| expn LSB bitstring := 110101
| 17*(17*((17*((17)^2)^2)^2)^2)^2

  17^47
| expn LSB bitstring := 111101
| 17*(17*(17*(17*((17)^2)^2)^2)^2)^2

  17^53
| expn LSB bitstring := 101011
| 17*((17*((17*(17)^2)^2)^2)^2)^2

  17^59
| expn LSB bitstring := 110111
| 17*(17*((17*(17*(17)^2)^2)^2)^2)^2

  17^61
| expn LSB bitstring := 101111
| 17*((17*(17*(17*(17)^2)^2)^2)^2)^2

  17^67
| expn LSB bitstring := 1100001
| 17*(17*(((((17)^2)^2)^2)^2)^2)^2

  17^71
| expn LSB bitstring := 1110001
| 17*(17*(17*((((17)^2)^2)^2)^2)^2)^2

  17^73
| expn LSB bitstring := 1001001
| 17*(((17*(((17)^2)^2)^2)^2)^2)^2

  17^79
| expn LSB bitstring := 1111001
| 17*(17*(17*(17*(((17)^2)^2)^2)^2)^2)^2

  17^83
| expn LSB bitstring := 1100101
| 17*(17*(((17*((17)^2)^2)^2)^2)^2)^2

  17^89
| expn LSB bitstring := 1001101
| 17*(((17*(17*((17)^2)^2)^2)^2)^2)^2

  17^97
| expn LSB bitstring := 1000011
| 17*(((((17*(17)^2)^2)^2)^2)^2)^2

  17^101
| expn LSB bitstring := 1010011
| 17*((17*(((17*(17)^2)^2)^2)^2)^2)^2

  17^103
| expn LSB bitstring := 1110011
| 17*(17*(17*(((17*(17)^2)^2)^2)^2)^2)^2

  17^107
| expn LSB bitstring := 1101011
| 17*(17*((17*((17*(17)^2)^2)^2)^2)^2)^2

  17^109
| expn LSB bitstring := 1011011
| 17*((17*(17*((17*(17)^2)^2)^2)^2)^2)^2

  17^113
| expn LSB bitstring := 1000111
| 17*((((17*(17*(17)^2)^2)^2)^2)^2)^2

  17^127
| expn LSB bitstring := 1111111
| 17*(17*(17*(17*(17*(17*(17)^2)^2)^2)^2)^2)^2

---

    # gawk profile, created Wed Nov  1 18:44:07 2023
    # Rule(s)

    32  _ = NF { # 32

    32      printf("\f %11s^%-5s \f\r\t| expn "\
                   "LSB bitstring := %*s\f\f\r\t| ", 
                   __ = $!!_, $!_ = $_, _ = !_, $_ = ____($_)) >> ("/dev/stdout")

    32  print($!(NF=(gsub(/[^01]+/,"")^_ * gsub(_,
              "(") + gsub(!_,(__)"*(")   * sub(/[^0-9]*$/, "")) "\n") >> ("/dev/stdout")
    }
    32  NF { # 32
    32      print
    }

    # Functions, listed alphabetically

   176  function ____(__, ___, _) {

   176      if (_ == "") {
    32          ___ = length(__ = (_) __)
    32          
                if ((__ *= _++ < (__ = int(__))) < ++_)
     1               return __

    31           _^= (___ += __+__ < (_+_)^_^_ \
                         ? _*--___ \
                         : _*___-_^(__ < (_*_^_^_)^_))

    31          if (__ < _)
                    ___ -= __ < (__ += __)

    31          if (_+_ <= __)
                    while ((_ += _+!++___)+_ <= __) { }

    31          return (_<__ ? ____(__ -= _, ___, _) \
                             : sprintf("%.*d", ___, !_)) !!_
            }

   144      return \
            (!--___ ? !! +__ : (__ += __)<_ \
                    ? ____(__,    ___, _) !_ \
                    : ____(__-=_, ___, _) !!_ )
    }

---

 1  664612181399115898708880371971190911115999473:   17^127
 2  845524284376730618271605950983326957831601937:   17^113
 3  064844117339164721813894369152081312807052497:   17^109
 4  484653439852384653016657073941702703504522673:   17^107

 5  861142522879753024803977642234575037448108913:   17^103
 6  861111219802352086591017223675552162759336017:   17^101
 7  0473579831553978069612272255088225422973416977:  17^97
 8  8773613853135920674790271577278638959930577297:  17^89
 9  5753941864386337047277655624320393544131928113:  17^83

10  5405894971825820887406128729967604798715821553:  17^79
11  3616579424626190476175164425470051684475963537:  17^73
12  7348154254064450486076730672752491528319986033:  17^71

13  6645123360044229169316532090500266429898209073:  17^67
14  1303399553000691993639838764066592228085041617:  17^61
15  7305548095339102740462421587418915544041816753:  17^59
16  4628544709931149045668299921272999257668224337:  17^53

17  9730995471356071936289402038303854098687372273:  17^47
18  0641272625991922463475488672334362661535406513:  17^43
19  5607755268602048174614102036928492604365174417:  17^41

20  3362095853201812742282475234995233875224247377:  17^37
21  139288917338851014461418017489467720433:         17^31
22  481968572106750915091411825223071697:            17^29
23  19967568900859523802559065713:                   17^23

24  239072435685151324847153:                        17^19
25  827240261886336764177:                           17^17
26  9904578032905937:                                17^13
27  34271896307633:                                  17^11

28  410338673:                                       17^7
29  1419857:                                         17^5
30  4913:                                            17^3
31  289:                                             17^2

32  17:                                              17

---

I kept all the extry layers of nesting brackets for clarity purposes only, and there are obvious ways to further streamline the statements :

say

   17^89   :=   LSB : 1 001 1 01
                          \
                MSB : 10 1 100 1
                      /  |  \
               [4] __/  [2]  \__ [8]

   LSB | 1  0..0..1   1  0..1 
        -----------------------
     => 17*(  (  (17*(17*( (17) ^ 2) ^ 2) ^ 2) ^ 2) ^ 2 ) ^ 2
     => 17*(      17*(17*  (17) ^  [4]  ) ^ 2) ^ [8]

     => 17*(17*(17*(17)^4)^2)^8