Интерполяция Лагранжа

 //преобразование точек c постоянной разностью по времени к функции методом интерполяции Лагранжа
        public double Points2LagFunc(PointF[] points,double X)
        {
            //функция лагранжа
            double L_y = 0;
 
            double h=points[1].X-points[0].X;
 
            double Yi = (X - points[0].X) / h;
 
            for (int i = 0; i < points.Length; i++)
            {
                double Li_x = 1;
                    for (int j = 1; j < points.Length; j++)
                    {
                           if (j != i)
                        Li_x *= (X - points[j].X) / (points[i].X - points[j].X);
                           //промежуточное значение функции
                    }
                    L_y += (Li_x * Yi);
            }
            return L_y;
        }

 //вспомогательный метод интерполяции лагранжа
        private double L(float xp,int i,PointF[] mas) 
        {
         
            int n = mas.Length;
            // числитель и знаменатель 
            double Chesl; 
            double Znam;
            Chesl = 1; Znam = 1;
            int k;
 
            // вычисление числителя
            for (k = 0; k!= n; k++ ) 
            {
            if ( k == i ) continue;
 
            double temp_chesl = xp - mas[k].X;
            Chesl =Chesl*temp_chesl;
            }
 
            // вычисление знаменателя
            for(k= 0; k!= n;k++) 
            {
            if (mas[i].X == mas[k].X) 
                continue; 
            // убираем, а то ноль в знаменателе
            float temp_znam = mas[i].X - mas[k].X;
            Znam =Znam* temp_znam;
            }
       
            double rez=(Chesl / Znam);
        return rez;
 
        }
        
        //основной метод интерполяции лагранжа
        public double InterpLagr(PointF[] mas,float X)
        {
            // начальное значение 
            double R = 0;
 
            // вычисляем значение интерполяционного многочлена
            for (int i = 0; i != mas.Length; i++)
            {
                R = R + (mas[i].Y * L(X, i, mas) );
            }
 
            return R;
        }

//вспомогательный метод интерполяции лагранжа
        private double L(float xp,int i,PointF[] mas) 
        {
            int tenmulcz = 0;
            int tenmulcch = 0;
            int n = mas.Length;
            // числитель и знаменатель 
            double Chesl; 
            double Znam;
 
            Chesl = 1; Znam = 1;
            int k;
            // вычисление числителя
            for (k = 0; k!=n; k++ ) 
            {
 
            if ( k == i ) continue;
 
            double temp_chesl = xp - mas[k].X;
            
                if (Chesl != 0)
            do
            {
                Chesl *= 10;
                tenmulcch++;
            }
            while (Math.Abs(Chesl) < 1 );
            //убираем множитель x - x(i)
            Chesl =Chesl* temp_chesl;
            }
 
            // вычисление знаменателя
            for(k= 0; k!= n;k++) 
            {
 
            if (mas[i].X == mas[k].X) 
                continue; 
            // убираем, а то ноль в знаменателе
            float temp_znam = mas[i].X - mas[k].X;
                if (Znam!=0)
            do
            {
                Znam *= 10;
                tenmulcz++;
            }
            while (Math.Abs(Znam) < 1);
            //убираем множитель x - x(i)
                Znam *= temp_znam;
            }
 
            int ten_stepen = tenmulcz - tenmulcch;
            double rez = (Chesl / Znam);
 
            if (ten_stepen < -50)
                return 0;
            else
            {
                rez*=Math.Pow(10, (double)ten_stepen);
                return rez;
            }
        }