Find the coefficients of the polynomial given its roots

Question

I am trying to write an algorithm which will find a(0),..., a(n-1), given the values of n, x_1, ..., x_n, a(n), such that:

a(n)*p^n + a(n-1)*p^(n-1) + ... + a(1)*p + a(0) = a(n)(p-x_1)(p-x_2)...(p-x_n)

for all real p.

After multiplying a(n)(p-x_1)(p-x_2) I've thought of using Viete's formulas to find the coefficients.

But it turns out writing the code down isn't as obvious as I expected.

I want to use only the basics in my code - that is loops, if-s addition and multiplication - no ready/ complex functions.

Here are the formulas:

First, I would like to emphasise that I only need a pseudocode, and I do not care about defining arrays for the root and coefficients. That's why I will just write a(n), xn. Oh, and I hope it won't bother you very much if I start indexing from i=1 not i=0 in order to be in synch with the mathematical notation. In order to start with i=0 I would have to renumerate the roots and introduce more brackets.

And this is what I've come up with so far:

a(n-1)=0;
for(i=1; i <= n; i++){
    a(n-1) = a(n-1) + x_i;
}

a(n-1) = -a(n)*a(n-1);
a(n-2)=0;

for(i=1; i <= n; i++){ 
    for(j=i; j <= n; j++){
        a(n-2) = a(n-2)+ x_i * x_j;
    }
}

a(n-2) = -a(n)*a(n-2);
a(n-3)=0;

for(i=1; i <= n; i++){
    for(j=i; j <= n; j++){
        for(k=j; k <= n; k++){
            a(n-3) = a(n-3)+ x_i * x_j * x_k;
        }
    }
}

a(n-3) = a(n)*a(n-3);

...

a(0)=1;
for(i=1; i<=n; i++){
    a(0) = a(0) * x_i;
}
if(n%2 == 0) a(0) = a(n) * a(0);
else a(0) = -a(n) * a(0);

As you can see, it doesn't look good.

I would like to link all those loops into one loop, because without I cannot write the full code, I cannot fill the gap between a(0) and a(n-j) for a fixed j.

Could you help me out?

This is what I have, based on Nico Schertler's answer:

for(i=1; i<=n; i++)
{a(i)=1; 
for(j=1; j <= n; j++)
{b(i)= clone( a(i) );
a(i) = a(i-1);
b(i) = x_j * b(i);
c(i) = a(i) - b(i);
}
}

Would it be the same if instead we wrote

for(i=1; i<=n; i++)
{a(i)=1; b(i)=1;
for(j=1; j <= n; j++)
{t = a(i) ;
a(i) = a(i-1);
b(i) = x_j * t;
c(i) = a(i) - b(i);
}
}

(this is how we for example swap two elements of an array, by keeping the value of a[i] in some variable t).

Which language are you writing in? When you say "As you can see, it doesn't look good.", what do you mean? Is the algorithm working as you wrote it down? — Gijs, Nov 08 '15 at 13:27
@Gijs I'm writing in C++, so there should be rectangular brackets [ ] instead of ( ). But by it doesn't look good I mean what I wrote in the line below. Precisely, I start with computing a[n-1], then a[n-2], ..., and at some point I need to find a[2], a[1], a[0] so I need to switch from a[n- something] to a[ something] — Don, Nov 08 '15 at 13:56
Another, probably less practical but interesting way to calculate your coefficients if the roots are integers is a technique called [Kronecker substitution](https://en.wikipedia.org/wiki/Kronecker_substitution). — biziclop, Nov 08 '15 at 14:29
Possible duplicate of [How to efficiently find coefficients of a polynomial from its roots?](http://stackoverflow.com/questions/33067755/how-to-efficiently-find-coefficients-of-a-polynomial-from-its-roots) — Lutz Lehmann, Nov 16 '15 at 17:41

score 7 · Accepted Answer · edited Sep 14 '16 at 19:00

7

You can create the polynomial incrementally.

Start with p = 1. I.e. a(0) = 1.

In order to add a root, you have to multiply the current polynomial by x - x_i. This is:

p * (x - x_i) = p * x - p * x_i

So you need to support three operations:

1. Multiplication by x

This is quite simple. Just shift all coefficients by one to the left. I.e.

a(i    ) := a(i - 1)
a(i - 1) := a(i - 2)
...
a(1    ) := a(0)
a(0    ) := 0

2. Multiplication by a scalar

This is equally simple. Multiply each coefficient:

a(i    ) *= s
a(i - 1) *= s
...

3. Subtraction

Just subtract the respective coefficients:

c(i    ) = a(i    ) - b(i    )
c(i - 1) = a(i - 1) - b(i - 1)
...

Altogether

Add root by root. First, clone your current polynomial. Then, do the operations as described above:

p := 1
for each root r
    p' = clone(p)
    multiply p with x
    multiply p' with r
    p := p - p'
next

edited Sep 14 '16 at 19:00

Anton Sherwood

123
6

answered Nov 08 '15 at 13:54

Nico Schertler

32,049
4
39
70

Thank you. Could you explain to me why we take p' = clone(p) and not just p' = p? – Don Nov 08 '15 at 14:07
1

@Don Because you need to multiply the original with something and the clone with something else. – biziclop Nov 08 '15 at 14:14
Of course if the multiplication doesn't happen in place but returns a copy anyway, you don't need to do a separate cloning. – biziclop Nov 08 '15 at 14:22
What would go wrong if we just wrote: p' := p; multiply p with x; multiply p' with r; p := p - p' ? – Don Nov 08 '15 at 14:23
1

@Don The result will be the constant 0 polynomial because at the end you'll do `p := p - p`. – biziclop Nov 08 '15 at 14:30
Could you take a look at the edit to my question above? – Don Nov 08 '15 at 14:58

rmojab63 · Answer 2 · 2017-04-09T10:41:52.983

A static function in c# for this purpose. The roots of x^4-11x^3+44x^2-76x+48 are {2,2,3,4} and given the argument

 roots = new Complex[4] {2, 2, 3, 4}

this function returns [48,-76,44,-11,1]

public static double[] FromRoots(Complex[] roots)
{
   int N = roots.Length;
   Complex[] coefs = new Complex[N + 1];
   coefs[0] = -roots[0];
   coefs[1] = 1.0;

   for (int k = 2; k <= N; k++)
   {
       coefs[k] = 1.0;
       for (int i = k - 2; i >= 0; i--)
       {
           coefs[i + 1] = coefs[i] - roots[k - 1] * coefs[i + 1];
       }
       coefs[0] *= -roots[k - 1];

       if (Math.IEEERemainder(k, 2) == 1)
           coefs[k] = -coefs[k];
    }

        double[] realCoefs = new double[N + 1];
        for (int i = 0; i < N + 1; i++)
            realCoefs[i] = coefs[i].Real; // Not sure about this part!

        return realCoefs;
    }

Find the coefficients of the polynomial given its roots

2 Answers2

1. Multiplication by x

2. Multiplication by a scalar

3. Subtraction

Altogether

Linked