Efficient calculation of polynomial coefficients from its roots

Question

I have the roots of a monic polynomial, i.e.

p(x) = (x-x_1)*...*(x-x_n)

and I need the coefficients a_n, ..., a_0 from

p(x) = x^n + a_{n-1}x^{n-1} + ... + a_0.

Does anybody know a computationally efficient way of doing this? If anybody knows a C/C++ implementation, this would actually be the best. (I already had a look at GSL but it did not provide a function.)

Of course, I know how to to it mathematically. I know, that the coefficient a_i is the sum of all products of the subsets with n-i elements. But if I would do it the dumb way, this means to iterate across all subsets, I would need

sum^{n-1}_{k=1} ( k choose n) * (k-1)

multiplications and

sum^n_{k=0} ( k choose n) - n

additions. Hence both terms grow with O(n!), which is too much computation to transform a list of n root into a list of n coefficients. I believe there must be some intelligent way to reuse most of the intermediate results, but I do not find one.

Friederike

You can recursively build the polynomial by convolution. If this is a very large polynomial, at some point FFT will beat the O(n^2) method. – Barbarese

Answer 1

You can do this in O(n^2) very easily if you incrementally build your polynomial. Let's define:

p_k(x) = (x-x_1)*...*(x-x_k)

That is p_k(x) is the multiplication of the first k (x-x_i) of p(x). We have:

p_1(x) = x-x_1

In other words the array of coefficients (a) would be (indices start from 0 and from left):

-x_1 1

Now assume we have the array of coefficients for p_k(x):

a_0 a_1 a_2 ... a_k

(side note: a_k is 1). Now we want to calculate p_k+1(x), which is (note that k+1 is the index, and there is no summation by 1):

p_k+1(x) = p_k(x)*(x-x_k+1)
=> p_k+1(x) = x*p_k(x) - x_k+1*p_k(x)

Translating this to the array of coefficients, it means that the new coefficients are the previous ones shifted to the right (x*p_k(x)) minus the k+1th root multiplied by the same coefficients (x_k+1*p_k(x)):

           0   a_0 a_1 a_2 ... a_k-1 a_k
- x_k+1 * (a_0 a_1 a_2 a_3 ... a_k)
-----------------------------------------
-x_k+1*a_0 (a_0-x_k+1*a_1) (a_1-x_k+1*a_2) (a_2-x_k+1*a_3) ... (a_k-x_k+1*a_k-1) a_k

(side note: and that is how a_k stays 1) There is your algorithm. Start from p_1(x) (or even p_0(x) = 1) and incrementally build the array of coefficients by the above formula for each root of the polynomial.

Paske

Umpf! If earth does not decide to swallow me up I volunteer to crawl under a rock and die. What ever, thank you. :-) – Friederike

@user2690527, this is really just a simple for loop inside another. Just 3 or 4 lines of code. Don't give up! – Paske

I know. I believe you misinterpreted my 1st comment to your solution.My comment is meant to to say, shame on me. Such an easy solution and I did not come up with it by myself. Two master degrees (computer science and mathematics), but I spent my whole afternoon on that problem. – Friederike

hahahaha, I understand now. Being a non-native speaker, sometimes some expressions turn out to have a different meaning than what I interpreted. Anyway, glad I could help. – Paske

Just my 2 ¢: shouldn't there be parenthesises around (x_k+1)? Multiplication has a higher precedence and multiplying by 1 doesn't change the value. SCNR. – Cirrostratus

No, that's x indexed by k+1. SO doesn't have support for tex math, so formatting math is hard. Technically, I could have written x_(k+1), but there is already enough parentheses there that more would just look more confusing. – Paske

In the answer, I had also specified this: note that k+1 is the index, and there is no summation by 1 – Paske