Edison.X. Blog

/*
     Largrange Interpolation Algorithm
     (1) give n, give x[0..n], give y[0..n]
     (2) Ln,k(x) = (x-x0)*(x-x1)*...*(x-xn) / (xk-x0)*(xk-x1)-...*(xk-xn)
     (3) P(x) = Ln1(x)f(x1)+Ln2(x)f(x2)+...+Lnn(x)f(xn)
     assume function = -4x^4 + 2x^3-3x^2+4x-1

                                                      -by EdisonX
*/
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#define N 4

double func(double x)
{
        double sum=0.0;
        double x2=x*x, x3=x2*x,x4=x3*x;
        sum = -4*x4+2*x3-3*x2+4*x-1;
        return sum;
}

double LargrangeInterpolation(double x,double *xn, double *yn, unsigned n)
{
        double numerator, denomenator; // 分子,分母
        double P; // result
        unsigned i, j;

        for(i=0, P=0.0; i<n; i++){
                if(yn[i]==0.0) continue;
                numerator=1.0, denomenator=1.0;
                for(j=0; j<n&& numerator!=0.0; j++){
                        if(i!=j) {
                                numerator*=(x-xn[j]);
                                denomenator*=(xn[i]-xn[j]);
                        }
                }
                P += numerator*yn[i]/denomenator ;
        }
        return P;
}

int main()
{
        unsigned i;
        double x=1.5;
        double xn[N+1]={1.234, 3.782, 4.563, 5.987, -1.157};
        double yn[N+1];

        for(i=0; i<N+1; i++) yn[i]=func(xn[i]); // reference point
        printf("correct: %lf\n", func(x));
        printf("Lagrange: %lf\n",LargrangeInterpolation(x, xn, yn, N+1));
        return 0;
}

output :

correct: -15.250000
Lagrange: -15.250000