22
$\begingroup$

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.

$\endgroup$
2
  • 6
    $\begingroup$ Perhaps polynomial chains are what you are looking for, see for example section 4.6.4 of The Art of Computer Programming. $\endgroup$ Jul 4, 2011 at 19:46
  • 2
    $\begingroup$ If the target solution is to be optimized for software running on commodity hardware, it would be wise to expand your criteria a bit, as modern CPUs are typically dominated by instruction latency as opposed to throughput, meaning that the smallest possible number of adders and multipliers will not necessarily be the fastest. Further, if the the degree is high, you will become bound by memory latency as well. K-th order Horner's form tends to work well on x86, especially with K=2. This is a classic latency vs throughput trade-off balancing problem. $\endgroup$
    – awdz9nld
    Feb 2, 2014 at 22:33

9 Answers 9

21
$\begingroup$

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.

$\endgroup$
9
  • 1
    $\begingroup$ Nathann Cohen, what do you mean by 'full information'? $\endgroup$
    – user9072
    Jul 4, 2011 at 16:20
  • 3
    $\begingroup$ My question probably was not accurate enough, so I appended a comment at the end. I am not looking for an algorithm that would be able to evaluate a polynomial, but an algorithm which -- given a polynomial -- would output "the best algorithm possible to evaluate this polynomial", using the least number of additions/muliplications possible. It would know the exact coefficients of this polynomial and can spend as much time as it likes trying to find out the best sequence of addition/multiplications to evaluate it. This way, this polynomial can be evaluated very often afterwards at no cost ! :-) $\endgroup$ Jul 4, 2011 at 16:23
  • 16
    $\begingroup$ You are basically asking about the arithmetical circuit complexity of the polynomial. For multivariate polynomials, this question leads to extremely hard problems in complexity, see en.wikipedia.org/wiki/Arithmetic_circuit_complexity for a start. Even in the univariate case, I doubt there is any general method to compute the complexity of specific polynomials. $\endgroup$ Jul 4, 2011 at 16:27
  • 3
    $\begingroup$ In many cases you might need to evaluate the polynomial many times for the same coefficient, in which cases it can pay to compile it down to a more efficient scheme than Horner's method. Clearly an algorithm to do so exist: iterate through all polynomial circuits in order of complexity until one is found to be equivalent to the problem at hand. Whether this can be done efficiently is another story. $\endgroup$
    – Arthur B
    Jul 25, 2016 at 21:23
  • 4
    $\begingroup$ You can actually improve upon Horner's scheme, reducing to n/2 multiplications. Knuth has a section on this in Volume 2, "4.6.4. Evaluation of Polynomials". $\endgroup$ Apr 23, 2021 at 23:23
11
$\begingroup$

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

$\endgroup$
5
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ The OP said that he is willing to neglect the differences between additions and multiplications, and treat them similarly. The question is answerable, and the answer is Horner's algorithm, look at the accepted answer. $\endgroup$
    – Alex M.
    Dec 19, 2021 at 8:17
  • 1
    $\begingroup$ Okay, looks like I did miss that the OP is willing to consider multiplications the same as additions. However, Horner's method is not the correct answer either. For instance it would use 256 multiplications to evaluate x^256, which is very suboptimal. $\endgroup$ Dec 20, 2021 at 14:58
4
$\begingroup$

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.

$\endgroup$
3
$\begingroup$

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

$\endgroup$
9
  • 1
    $\begingroup$ Regarding "the cheapest way": looking things up in a table comes with a cost, and the cost will generally be at least logarithmic in the table size $\endgroup$ Oct 1, 2013 at 16:19
  • $\begingroup$ Sasho, indeed it does. However, table-type access has been optimized at the compiler and machine levels, and I don't know of anything faster. Do you? Gerhard "Ask Me About System Design" Paseman, 2013.10.01 $\endgroup$ Oct 1, 2013 at 17:15
  • 1
    $\begingroup$ But that also seems to assume one knows at what points the polynomial needs to be evaluated and in general, once the number of possible points is big enough, an actual evaluation algorithm will be faster. While look up in small tables may be very fast, this is not at all the case for very large tables when random access actually becomes an unreasonable assumption. $\endgroup$ Oct 1, 2013 at 17:45
  • $\begingroup$ If we are worried about table size, we should also be worried about bit length. I suspect the claim (inferred from your comment) that for any polynomial P, there will be n_0 such that for all n >n_0 evaluating P(n) will be faster than looking up P(n) in a "Google-sized" table. (Here I mean tables accommodating 10^m values for m at most 40. m about 80 is likely physically impossible in this millenium.) Even with idealized infinite tables with extrapolated access times, I doubt the inference you suggest. Gerhard "Recall The Unlimited Preprocessing Power" Paseman, 2013.10.01 $\endgroup$ Oct 1, 2013 at 17:53
  • 1
    $\begingroup$ Yes, and if I knew how to give a better answer to the question, then I would. If you know of a better way to get fast performance given unlimited preprocessing resources, please let us know. While table lookup may be theoretically unsatisfactory to you, it may turn out to be useful to the poster, especially if there is an access pattern to use in optimization. Gerhard "Willing To Agree To Disagree" Paseman, 2013.10.02 $\endgroup$ Oct 2, 2013 at 19:00
1
$\begingroup$

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

$\endgroup$
1
  • 1
    $\begingroup$ I also cite an answer involving differencing below. Can you say how forward differencing differs from kth order differences? Gerhard "Not Sure About Differing Differences" Paseman, 2011.07.14 $\endgroup$ Jul 14, 2011 at 17:31
1
$\begingroup$

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

$\endgroup$
1
$\begingroup$

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.

$\endgroup$
5
  • 4
    $\begingroup$ The problem with your answer is that it states only vague, generic truths, and opinions. If you had read the other answers, you would have noticed that Horner's algorithm has been indicated as the fastest available algorithm (provably so), with an explicitly given computational complexity. "Sparseness", "NP-completeness"... these sound like buzzwords to me, in this context. This is a low-quality answer. $\endgroup$
    – Alex M.
    Dec 19, 2021 at 8:13
  • 1
    $\begingroup$ There is so little information in the question. For example, if the $a_i$'s are small integers the polynomial may be made "sparse", in some computationally useful sense. But the coefficients may be high precision floats for all we know. $\endgroup$ Dec 19, 2021 at 14:20
  • $\begingroup$ @Alex M. Not really. The majority of answers quote Horner's rule, but if you carefully read the question and the comments that is not the answer he was looking for. Yes, you can evaluate an arbitrary polynomial in exactly N additions and multiplications, but if you have, say a polynomial with many repeated roots like in the example provided in the question you can do better than Horner by just factoring the polynomial in advance and using exponentiation by squaring. $\endgroup$ Dec 20, 2021 at 14:05
  • $\begingroup$ @Alex M. This suggests to me that you can systematically try every way of splitting the polynomial into sums of polynomials with enough repeated roots, something that could be in principle be tried using dynamic programming for polynomials with integer coefficients and would in general be a task with exponential complexity. $\endgroup$ Dec 20, 2021 at 14:19
  • $\begingroup$ @Alex M. There are other possible factorizations too, for instance: x^3+4x^2+6x+4 = (((x+1)^2)^2 -1)/x saves operations in some sense. Of course division is usually significantly more expensive than multiplication, but there are analogous examples where the tradeoff makes sense. $\endgroup$ Dec 20, 2021 at 14:29
-2
$\begingroup$

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
$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.