http://www.velocityreviews.com/forums/(E-Mail Removed) wrote:

> I am looking for a quick C program that takes n+1 pairs of values

> (integers) (a_i, f(a_i)), i=0,...,n, generates the coefficients

> \alpha_i, i=0,...,n of the polynomial of degree n that fits these

> points and then outputs the polynomial with indeterminant in the form

> \sum_i=0^n\alpha_i X^i, where X is just some indeterminant.
In general, you need an n+1 degree polynomial to fit n points. You are

specifying the points as 0, ..., n which is already n+1 pointts. I'm

going to assume you mean to specify the points as 0, ..., n-1.

There is a well known interpolation mechanism called "Lagrange

Interpolation":

Let LF(X) = \prod_i=0^n-1\(X - a_i), and let F_i(X) = LF(X)/(X - a_i),

where F_I(a_i) is calculated by taking the limit as X->a_i; i.e.,

cancel out the factors.

Then the polynomial you are looking for is:

\sum_i=0,n-1\f(a_i)*F_i(X)/F_i(a_i)

The reason is that F_i(a_j) = 0 for any i not equal to j.

So then your problem reduces to how do you write a program to expand

the polynomials: f(a_i)/F_i(a_i) * F_i(X)

The way you would want start is to break this down into its

coefficients by writing a function that computed:

double coefficient_F(int i, int p, double * a, int n);

which is the pth coefficient of F_i as defined on the vector a of size

n. The idea is that you would do it recursively, by peeling off the

last factor one at a time:

double coefficient_F(int i, int p, double * a, int n) {

if (p == n-(i < n)) return 1.0;

if (p < n) return 0.0;

if (i = n-1) n--;

return coefficient_F(i, p-1, a, n-1)

- a[n-1] * coefficient_F(i, p, a, n-1);

}

Notice how the first termination condition basically looks for when you

have exactly p factors of (X - a_*) in this factor, the second just

deals with the case where you asked for coefficient that was negative,

or otherwise too low, and otherwise, it performs the one recursive

expansion corresponding to the last factor in the polynomial.

Then you need the actually need to be able to compute value F_i(a_i)

which is just a simple for loop (I'll leave that as an exercise for

you.)

Then you have what you need to perform a polynomial expansion of each

term of the sum, then you need to just sum the coefficients.

--

Paul Hsieh

http://www.pobox.com/~qed/
http://bstring.sf.net/