Symbolic Polynomial Interpolator in C

[daniel.wolff]

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.

Any hints?

Thanks

 

[Walter Roberson]

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.

Any hints?

You don't want a symbolic manipulation for that.

http://mathworld.wolfram.com/Condensation.html
http://mathworld.wolfram.com/Gauss-JordanElimination.html
http://mathworld.wolfram.com/GaussianElimination.html

 

[William Hughes]

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.

Any hints?

This is not really a C question.

Are you sure you want to do this?

Are you sure you want the interpolating polynomial
not a fitted polynomial?

Do you know that the calculated values of the coefficients
are very likely to be inexact and may be meaningless?

If so you can find code that does this (e.g. the routine
polcos from Numerical Recipes in C, (you need to write the
output bit but this is easy)).

On the whole I doubt that you really want to do this,
but I could be wrong.

I suggest you ask for help in sci.num-analysis, but be
sure to ask about what you want to do, not how you want
to do it.

-William Hughes

 

[daniel.wolff]

Thanks for the links. I have in mind Newton's divided difference
approach, or Lagrange's formula and I am not sure how calculating the
determinant is useful here but I will think about it. Anyway, the major
difficulty I am having is formating the output. Newton's method, which
appears to be the computably nicer of the two, quickly gives the
coefficients \lambda_i, where the poly which fits the (a_i,f(a_i)) is
given by \lamda_0+\lambda_1 (X-a_0)+\lambda_2
(X-a_0)(X-a_1)+...+\lambda_n (X-a_0)...(X-a_{n-1}). But to get it in
the form \sum_i=0^n\alpha_i X^i I need to multiply through and simplify
the above. (I would could hard code this in for this particular
example. But, a side question is of course has anyone seen C code that
could handle this symbolically? It is likely overkill for this example,
but I am curious in general.) Okay, but the embarassing question is
once we have the \alpha_i, how to we format it pretty in the desired
form? I would love to see some actual code.
Thanks.

 

[Walter Roberson]

:Thanks for the links. I have in mind Newton's divided difference
:approach, or Lagrange's formula and I am not sure how calculating the
:determinant is useful here but I will think about it.

Looks like I may have misunderstood the question.

If you are going to use Newton's divided difference, you likely do not
need symbolic calculations. See e.g.,

http://kr.cs.ait.ac.th/~radok/math/mat7/pseudo.htm#DIVIDED

 

[daniel.wolff]

Thanks William, I am trying to implement Kronecker's algorithm for
factoring integer polys in multi indeterminants. I would then like to
see if an extension is possible for some standard finite extension
rings. (I know it is slow as hell but I just want to see it built.) The
first step is to do this in a single indeterminants but I would like
keep the code as general as possible for later steps.

 

[William Hughes]

Thanks William, I am trying to implement Kronecker's algorithm for
factoring integer polys in multi indeterminants. I would then like to
see if an extension is possible for some standard finite extension
rings. (I know it is slow as hell but I just want to see it built.) The
first step is to do this in a single indeterminants but I would like
keep the code as general as possible for later steps.

Wow this is now so off topic for comp.lang.c that it makes
discusions about how to clear the screen seem on topic [1].

I'm afraid you've lost me. Try sci.math

I will point out that your interpolating polynomials
will not in general have integer coefficients.
Does this matter? (I don't know.)

As far as output, if you have a floating point vector a of
n coefficients in increasing order you want something like

printf("%f",a[0]);

for(i=1;i<n;i++)
{
printf(" + %f X^%d",a,i);
}

printf("\n");

which will produce something like

3.2734522 + 5.8012738 X^1 + 7.9245723 X^2 + 4.8274643 X^3

Not very pretty, but about the best you can do without using
extensions. On the other hand I still am not sure why
you would want to output this in the first place.


-William Hughes



[1] Using your basic "I've seen hills that make that hill
look like a valley" logic.

 

[daniel.wolff]

Thanks William, I will give the other groups a try. Although my goal
may be a bit outside of this group, what I am looking for, for now,
might be in here. Anyway, thanks for the tip. I am just looking for
easily extendable interpolation code and some pretty output stuff
(thanks for the above). And, about the coefficients: it is easy to show
that interpolating polys should have coefficients within the original
ring from which the pairs are selected. (So, interpolating integer
points gives an integer polynomial.) The only possible trouble in
practice that might arise is overflow. Cheers.

 

[Markus Moll]

Hi

And, about the coefficients: it is easy to show
that interpolating polys should have coefficients within the original
ring from which the pairs are selected. (So, interpolating integer
points gives an integer polynomial.) The only possible trouble in
practice that might arise is overflow. Cheers.

Maybe I have misunderstood you, but I'd like to see that proof.
Trying to interpolate the points (0,0) and (5,1) by a line yields the
polynomial p(x) = 1/5 x, I think...

Cheers
Markus

 

[daniel.wolff]

Sorry. You're right, if you start with a field, the coefficients will
remain in that field. But, if you start with an integral domain D, then
the best you can show, in general, is the coefficients are in quotient
field Q_D. Thanks.

 

[websnarf]

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.

 

[Eric Sosman]

In general, you need an n+1 degree polynomial to fit n points.

ITYM a degree n polynomial to fit n+1 points.

(Consider an easy case: To fit two points you need a
straight line, that is, a 1-degree polynomial. Using a
cubic would be redundantly superfluous overkill ...)

 

[Julian V. Noble]

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.

Any hints?

Thanks

Use Gram polynomials, which for integers in a given range (that is,
equally spaced x-values) are actually Chebyshev polynomials. See
Ralston, Intro to Numerical Analysis.

The chief virtue of Gram polynomials is you don't have to invert
matrices. The disadvantage is that the result is always in the
form

f\left( x \right) \approx
\sum\limits_{k = 1}^m {a_k C_k \left( x \right)}

--
Julian V. Noble
Professor Emeritus of Physics
(e-mail address removed)
^^^^^^^^^^^^^^^^^^
http://galileo.phys.virginia.edu/~jvn/

"For there was never yet philosopher that could endure the
toothache patiently."

-- Wm. Shakespeare, Much Ado about Nothing. Act v. Sc. 1.