using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using System.IO;

namespace Interpolation
{
    class Program
    {
        #region Нужные переменные

        enum Variant
        {
            A = 9, B = 29, C = 25
        }

        // номер варианта
        static Variant variant;

        const int N = 6;
        // -----------------------------------------
        //матрица разделённых разностей размерность 6 X 7
        static double[,] divided_differences = new double[N,N + 1];
        static double x = 1;
        static double f;

        // для работы с файлами
        public static StreamWriter sw = new StreamWriter("C:\\Users\\Diana\\Desktop\\Learning\\Численные методы\\Interpolation_Output.txt");
        #endregion

        #region Вычисление значений функций и производных

        // функция
        public static double Function(double cur_x)
        {
            switch (variant)
            {
                case (Variant.A):
                    {
                        return (Math.Pow(Math.E, cur_x) - 6*cur_x + 3);
                    }
                    break;

                case (Variant.B):
                    {
                        return (Math.Pow(3, cur_x) + 2*cur_x*cur_x - 3*cur_x + 1);
                    }
                    break;

                case (Variant.C):
                    {
                        return (Math.Pow(3, cur_x) + 2 - cur_x);
                    }
                default:
                    return Int32.MaxValue;
                    break;
            }
        }

        // шестая производная
        public static double Sixth_Derivative(double cur_x)
        {
            switch (variant)
            {
                case (Variant.A):
                    {
                        return (Math.Pow(Math.E, cur_x));
                    }
                    break;

                case (Variant.B):
                    {
                        return (Math.Pow(Math.Log(3), 6) * Math.Pow(3, cur_x));
                    }
                    break;
                case (Variant.C):
                    {
                        return (Math.Pow(Math.Log(3), 6) * Math.Pow(3, cur_x));
                    }
                default:
                    return Int32.MaxValue;
                    break;
            }
        }
        #endregion

        #region Вывод матрицы

        public static void Output_Matrix()
        {
            int j = 0;
            for (int i = 0; i < N; i++)
            {
                while (j < N + 1)
                {
                    sw.Write(String.Format("{0,11:0.0000000}", divided_differences[i, j]) + "\t");
                    j++;
                }
                sw.WriteLine();
                j = 0;
            }
        }
        #endregion

        #region Интерполяционная формула Ньютона

        // матрица разделённых разностей
        public static void Divided_Differences_Matrix()
        {
            // заполняем первый и второй столбцы
            // первый столбик - значения иксов от 1 до 2 с шагом 0.2
            // второй столбик - значения функции от каждого значения из первого столбца соответственно
            for (int i = 0; i < N; i++, x += 0.2)
            {
                divided_differences[i, 0] = x;
                divided_differences[i, 1] = Function(x);
            }

            // заполняем по столбцам (в матрице 7 столбцов и 6 строк)
            for (int column = 2; column < divided_differences.GetLength(1); column++)
                for (int str = 0; str < divided_differences.GetLength(0) - column + 1; str++)
                    divided_differences[str, column] = (divided_differences[str + 1, column - 1] - divided_differences[str, column - 1])
                                                     / (divided_differences[column + str - 1, 0] - divided_differences[str, 0]);
        }

        // многочлен Ньютона
        public static double NewtonPolinomial(double cur_x)
        {
            double newton_polinomial = 0;
            // добавляем f(1) к многочлену
            newton_polinomial += divided_differences[0, 1];

            // для каждого значения из первой строки таблицы разделённых разностей
            // добавляем соответствующий многочлен
            for (int column = 2; column < divided_differences.GetLength(1); column++)
                newton_polinomial += divided_differences[0, column] * PolinomialForEachColumn(column, cur_x);

            return newton_polinomial;
        }

        // подсчёт полинома для каждого столбца матрицы разделённых разностей
        public static double PolinomialForEachColumn(int cur_column, double cur_x)
        {
            int count_of_polinomials = 0;
            int str = 0;
            double new_polinomial = 1;

            while (count_of_polinomials != cur_column - 1)
            {
                new_polinomial *= (cur_x - divided_differences[str, 0]);
                str++;
                count_of_polinomials++;
            }

            return new_polinomial;
        }

        // для оценки погрешности
        public static double GetWPolinomial(double cur_x)
        {
            double w_polinomial = 1;

            for (int i = 0; i < N - 1; i++)
                w_polinomial *= cur_x - divided_differences[i, 0];

            return Math.Abs(w_polinomial);
        }

        // нахождение максимального значения шестой производной на отрезке от 1 до 2
        public static double MaxSixthDerivative()
        {
            double cur_x = 1;
            double max_sixth_derivative = Int32.MinValue;
            double sixth_derivative = 0;
            for (int i = 0; i < N - 1; i++, cur_x += 0.2)
            {
                sixth_derivative = Math.Abs(Sixth_Derivative(cur_x));
                if (sixth_derivative > max_sixth_derivative)
                    max_sixth_derivative = sixth_derivative;
            }

            return sixth_derivative;
        }

        #endregion

        static void Main(string[] args)
        {
            // менять вариант вручную! A - 9, B - 29
            variant = Variant.C;

            // --------- Интерполяционная формула Ньютона ------------
            sw.WriteLine("ИНТЕРПОЛЯЦИОННАЯ ФОРМУЛА НЬЮТОНА");

            // вывод матрицы разделённых разностей
            Divided_Differences_Matrix();
            sw.WriteLine("Матрица разделённых разностей имеет вид:");
            Output_Matrix();
            sw.WriteLine();
            x = 1.1;
            // вывод таблицы  x - f(x) - Pn(x) - Delta - Оценка
            sw.WriteLine("x - f(x) - Pn(x) - Delta - Оценка");
            for (int i = 0; i < N - 1; i++, x += 0.2)
            {
                double function = Function(x);
                double polinomial_newton = NewtonPolinomial(x);
                double delta = Math.Abs(Function(1 + (i + 0.5) * 0.2) - NewtonPolinomial(1 + (i + 0.5) * 0.2));
                double error_of_estimate = (MaxSixthDerivative() * GetWPolinomial(x)) / 720;
                sw.Write(x + "\t" + String.Format("{0,11:0.00000}", function) + "\t"
                                  + String.Format("{0,11:0.00000}", polinomial_newton) + "\t"
                                  + String.Format("{0,15:#.########E-00}", delta) + "\t"
                                  + String.Format("{0,15:#.########E-00}", error_of_estimate));
                sw.WriteLine();
            }
            // -------------------------------------------------------



            sw.Close();
        }
    }
}