Полином Лагранжа (интерполяция Лагранжа)

Привожу свой код. Должен работать. Подробности читай в камментах :)

// Интерполяция Лагранжа.
// Компилятор: Visual C++ 6.0

#include "stdafx.h"

// таблица из книги по вычислительной математике 
float x[6] = {1.5, 1.54, 1.56, 1.60,1.63, 1.70}; 
float y[6] = {3.873,3.924,3.950,4.00,4.037,4.123}; 

/* Функция, вычисляющая коэффициенты Лагранжа
x - аргумент функции
n - степень многочлена (или число x-ов)
i - номер узла
*/

float L(float xp,int n,int i) {
// числитель и знаменатель 
float Chesl; 
float Znam;

Chesl = 1; Znam = 1;

int k;
// вычисление числителя
for (k = 0; k!= n; k++ ) {

if ( k == i ) continue; 
// убираем множитель x - x(i)
Chesl *= xp - x[k];
}
// вычисление знаменателя
for(k= 0; k!= n;k++) {

if (x[i] == x[k]) continue; 
// убираем, а то ноль в знаменателе

Znam *= x[i] - x[k];
}


return Chesl/Znam; 

}
int main(int argc, char* argv[])
{

// вычисляем степень полинома
int n = sizeof(y)/sizeof(float); 

// начальное значение 
float R = 0; 

// произвольная точка для проверки 
float px = 1.55;

// вычисляем значение интерполяционного многочлена
// в точке должно получиться 3.937 
for (int i = 0; i != n; i++) { 

R += y[i]*L(px,n,i);
}

printf("Результат : %f ",R); 

// вывод результата
return 0;
}