Algorithm in C to calculate coefficients of polynomial using Lagrange interpolation

Question

I've been stuck on this for a while now. I'm writing an algorithm in C to pull out the coefficients of a polynomial using Lagrange's interpolation method.

My code partially works, for instance if we do the first example here http://en.wikipedia.org/wiki/Lagrange_polynomial#Example_1 then the code can print out the first 2 coefficients (0 and 4.834848)

Similarly with example 3 on that article, it will print the 2 coefficients 6 and -11.

I need to be able to accurately get all the coefficients from any set of points. Please advise on the alterations required to the code.

Thanks in advance!

Updated with latest code, 7:57PM, GMT on August 5th. 9 coefficients now working, getting ugly looking. Will investigate iterative process for n degrees tomorrow!

#include<ncursesw/ncurses.h>
#include<math.h>
#include <stdio.h>
#include <string.h>
#include <errno.h>
#include <stdlib.h>
#define MAX 200
float coeff[MAX], coefftwo[MAX], coeffthree[MAX], coefffour[MAX];
int count;
void main()
{
int n,i,j ;
char ch;
float x[MAX],y[MAX],fp2, coeff1, coeff2;
printf("\n\nn = ");
scanf("%i", &count);

    for(i=0; i < count; i++)
    {
        printf("\n\n The value of x%i= ", i);
        scanf("%f",&x[i]);
        printf("\n The value of f(x%i)= ", i);
        scanf("%f",&y[i]);
    }
    for(i=0;i<count;i++)
    {
        coeff1 = 1.0;
        coeff2 = 0.0;
        coeff3 = 0.0;
        coeff4 = 0.0;
        coeff5 = 0.0;
        coeff6 = 0.0;
        coeff7 = 0.0;
        coeff8 = 0.0;
        coeff9 = 0.0;
        for(j=0; j<count; j++)
        {
            if(i!=j) {
                coeff1 = coeff1 * (array[i]-array[j]);
                coeff2 -= array[j];
                for (int k=j; k < count; k++) {
                if ((j!=k) && (k!=i)) {
                coeff3 += array[j] * array[k];
                    for(int l=k; l < count; l++) {
                    if ((l!=j) && (l!=k) && (l!=i))     {
    coeff4 -= array[j] * array[k] * array[l];

    for (int m = l; m < count; m++) {
        if ((m!=l) && (m!=k) && (m!=j) && (m!=i)) {                 coeff5 += array[j] * array[k] * array[l] * array[m];
            for (int n = m; n < count; n++) {
            if ((n!=m) && (n!=l) && (n!=k) && (n!=j) && (n!=i)) {
            coeff6 -= array[j] * array[k] * array[l] * array[m] * array[n];
            for (int o = n; o < count; o++) {
            if ((o!=n) && (o!=m) && (o!=l) && (o!=k) && (o!=j) && (o!=i)) {
            coeff7 += array[j] * array[k] * array[l] * array[m] * array[n] * array[o];
            for (int p = o; p < count; p++) {
            if ((p!=o) && (p!=n) && (p!=m) && (p!=l) && (p!=k) && (p!=j) && (p!=i)) {
            coeff8 -= array[j] * array[k] * array[l] *array[m] *array[n] * array[o] * array[p];
            for (int q = p; q < count; q++) {
            if ((q!=p) && (q!=o) && (q!=n) && (q!=m) && (q!=l) && (q!=k) && (q!=j) && (q!=i)) {
            coeff9 += array[j] * array[k] * array[l] * array[m] * array[n] * array[o] * array[p] * array[q];
            }
            }
            }
            }
            }
            }
            }
            }
                    }               
        }
                    }
                }
            }
        }    
    }
}       
        coeff[i] = y[i] / coeff1;
        coefftwo[i] = y[i] * coeff2 / coeff1;   
        coeffthree[i] = y[i] * coeff3 / coeff1;
        coefffour[i] = y[i] * coeff4 / coeff1;
        coefffive[i] = y[i] * coeff5 / coeff1;
        coeffsix[i] = y[i] * coeff6 / coeff1;
        coeffseven[i] = y[i] * coeff7 / coeff1;
        coeffeight[i] = y[i] * coeff8 / coeff1;
        coeffnine[i] = y[i] * coeff9 / coeff1;
    }
float coefficientone = 0.0;
float coefficienttwo = 0.0;
float coefficientthree = 0.0;
float coefficientfour = 0.0;
float coefficientfive = 0.0;
float coefficientsix = 0.0;
float coefficientseven = 0.0;
float coefficienteight = 0.0;
float coefficientnine = 0.0;
for (int i = 0; i< count; i++){
        coefficientone = coefficientone + coeff[i];
        coefficienttwo = coefficienttwo + coefftwo[i];
        coefficientthree = coefficientthree + coeffthree[i];
        coefficientfour = coefficientfour + coefffour[i];
        coefficientfive = coefficientfive + coefffive[i];
        coefficientsix = coefficientsix + coeffsix[i];
        coefficientseven = coefficientseven + coeffseven[i];
        coefficienteight = coefficienteight + coeffeight[i];
        coefficientnine = coefficientnine + coeffnine[i];
}
printf("coefficient 1 = %f\n", coefficientone);
printf("coefficient 2 = %f\n", coefficienttwo);
printf("coefficient 3 = %f\n", coefficientthree);
printf("coefficient 4 = %f\n", coefficientfour);
printf("coefficient 5 = %f\n", coefficientfive);
printf("coefficient 6 = %f\n", coefficientsix);
printf("coefficient 7 = %f\n", coefficientseven);
printf("coefficient 8 = %f\n", coefficienteight);
printf("coefficient 9 = %f\n", coefficientnine);

}

Answer 1

Your algebra is simply wrong, and that fact is hidden by poorly chosen variable names.

When you calculate the contribution of the i^th basis polynomial (never mind y for now) what variable represents the coefficient of the x² term? It's coeff3. And you don't calculate it correctly.

Take a simpler case. Suppose you want to work out (x+a)(x+b)(x+c)(x+d). The first term is x⁴, easy. The next is (a+b+c+d)x³, not too bad. The next is (ab + ac + ad + bc + bd + cd)x², and now it's clear that a single loop won't do the job. It's worth taking the time to make sure you can write code that handles the simple problem correctly, before you try the more complex one. You need something like this:

for(unsigned int j=0 ; j<count ; ++j)
{
  ...
  coeff2 -= x[j];
  for(unsigned int k=j ; k<count ; ++k)
  {
    if(j!=k && k!=i)
      coeff3 += x[j] * x[k];
    ...
  }
}

That should be enough to get you started.