Symbolic Polynomial Interpolator in C

Discussion in 'C Programming' started by daniel.wolff@csfb.com, Aug 10, 2005.

  1. Guest

    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
     
    , Aug 10, 2005
    #1
    1. Advertising

  2. In article <>,
    <> 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.


    >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
    --
    Entropy is the logarithm of probability -- Boltzmann
     
    Walter Roberson, Aug 10, 2005
    #2
    1. Advertising

  3. 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.
    >
    > 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
     
    William Hughes, Aug 10, 2005
    #3
  4. Guest

    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.
     
    , Aug 10, 2005
    #4
  5. In article <>,
    <> wrote:
    :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
    --
    Feep if you love VT-52's.
     
    Walter Roberson, Aug 10, 2005
    #5
  6. Guest

    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.
     
    , Aug 10, 2005
    #6
  7. wrote:
    > 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.
     
    William Hughes, Aug 10, 2005
    #7
  8. Guest

    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.
     
    , Aug 10, 2005
    #8
  9. Markus Moll Guest

    [OT] Re: Symbolic Polynomial Interpolator in C

    Hi

    wrote:

    > 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
     
    Markus Moll, Aug 10, 2005
    #9
  10. Guest

    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.
     
    , Aug 10, 2005
    #10
  11. Guest

    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/
     
    , Aug 10, 2005
    #11
  12. Eric Sosman Guest

    wrote:
    > 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.


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

    --
    Eric Sosman
    lid
     
    Eric Sosman, Aug 11, 2005
    #12
  13. 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.
    >
    > 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

    ^^^^^^^^^^^^^^^^^^
    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.
     
    Julian V. Noble, Aug 12, 2005
    #13
    1. Advertising

Want to reply to this thread or ask your own question?

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. Manfred Balik

    polynomial division remainder

    Manfred Balik, May 12, 2004, in forum: VHDL
    Replies:
    5
    Views:
    3,322
    Mike Treseler
    May 18, 2004
  2. S.Muralidharan

    polynomial

    S.Muralidharan, Nov 4, 2004, in forum: VHDL
    Replies:
    1
    Views:
    632
    Tuukka Toivonen
    Nov 4, 2004
  3. aDawg

    Linear Interpolator

    aDawg, Jun 23, 2006, in forum: VHDL
    Replies:
    0
    Views:
    598
    aDawg
    Jun 23, 2006
  4. Grant Edwards

    Looking for triangulator/interpolator

    Grant Edwards, May 27, 2006, in forum: Python
    Replies:
    13
    Views:
    606
    Travis E. Oliphant
    May 28, 2006
  5. Chen L.
    Replies:
    2
    Views:
    886
    Dik T. Winter
    Jul 6, 2004
Loading...

Share This Page