/*  Lagrangian Differentiation  */

#include <stdio.h>
#include <math.h>

#define PI  3.1415926535

double lagdif(double*, double*, int, double);

// µ¥ÀÌÅÍ°ª 1
double x1[16]={  0.3378,   1.0134,  1.6890,  3.043,   4.3915,
                 5.7427,   7.0939,  8.4451,  9.7964, 11.1476,
                17.9037,  24.6598, 31.4159, 38.1720, 44.9282,
                51.6843, };

double y1[16]={ 0.,    0.,   2.05, 11.95,
               19.75, 23.,  24.2,  24.9,
               25.1,  25.2, 25.4,  25.5,
               25.75, 26.2, 26.4, 23.};   

// µ¥ÀÌÅÍ°ª 2
double x2[29]={  0.3378,   1.0134,   1.6890,   3.043,     4.3915,
                 5.7427,   7.0939,   8.4451,   9.7964,   11.1476,
                17.9037,  24.6598,  31.4159,  38.1720,   44.9282,
                51.6843,  58.4404,  65.1965,  71.9526,   78.7087,
                85.4648,  92.2209,  98.9771, 105.7332,  112.4893,
               119.2454, 126.0015, 132.7576, 139.5137 };

double y2[29]={   0.,   0.,    0.1,   3.2,   5.8,
                  6.9,  7.75,  7.85,  8.1,   8.35,
                  9.2,  9.9,  10.5,  11.,   11.55,
                 12.,  12.4,  12.8,  13.05, 13.5,
                 13.8, 14.15, 14.45, 14.75, 15.05,
                 15.3, 15.65, 15.85, 16.1 };

void main()
{
   double result;

   for(int i=0; i<16; i++) {
      result = lagdif(x1, y1, i, x1[i]);
      printf("%10.4e   ", result);
      if(!( (i+1)%5) ) printf("\n");
   }

   for(int j=0; j<29; j++) {
      result = lagdif(x2, y2, j, x2[i]);
      printf("%10.4e   ", result);
      if(!( (j+1)%5) ) printf("\n");
   }
}

double lagdif(double x[], double y[], int n, double xx)
{
   int i, j, k;
   double xi, sm0, sm1, sm2, w;

   if(n < 2 || xx < x[0] || xx > x[n]) return(999);
   for(i = 0; i < n; i++)
      if(x[i+1] <= x[i]) return(999);
   w = 0.0;
   for(i = 0; i <= n; i++){
      sm0 = 0.;
      sm2 = 1.;
      xi = x[i];
      for(j = 0; j <= n; j++) {
         if(i != j) {
            sm2 *= (xi - x[j]);
            sm1 = 1.;
            for(k = 0; k <= n; k++)
               if((k != i) && (k != j))
                  sm1 *= (xx - x[k]);
                  sm0 += sm1;
         }
      }
      w += y[i] * sm0 / sm2;
   }
   return w;
        
}