Find the coefficients of the polynomial given its roots

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

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