using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace Lab1_16
{
    class Spline
    {
        double[,] nodes; // [номер точки, (0 -- точки X, 1 -- точки Y)] 
        double[] b; // Вільні члени в системі рівнянь відносно інтерполяційних коефіцієнтів
        double[,] a; // Числа при коефіцієнтах сплайну
        double[] s; // Коефіцієнти сплайну
        double[,] dd; // [номер точки, (0 -- точки X, 1 -- значення Y'')]
        int n; // Кількість вузлів
        public Spline(ref double[,] nodes, ref double[,] dd) // Конструктор
        {
            n = nodes.GetLength(0);
            a = new double[n, n];
            b = new double[n];
            s = new double[n];
            this.nodes = (double[,])nodes.Clone(); // Копіюємо масив вузлів
                this.dd = (double[,])dd.Clone(); // Копіюємо масив з додатковими умовами
            
            SplineCoeffsCalc();
        }
        private void SplineCoeffsCalc() // Розрахунок коефіцієнтів сплайну
        {
            int info;
            alglib.densesolverreport rep;
            /* Розрахунок коефіцієнтів за основники рівняннями */
            
            for (int k = 0, i = 1; k < n ; k++, i++) // Розрахунок коефіцієнтів за формулами (11.69), (11.71), (11.72)[Ю.А. Дубинский, Н.В. Копченова "Вычислительные методы для инженеров", с. 336]
            {
                if (k < (n - 2)) // (11.69)
                {
                    a[k, i - 1] = 1 / h(i);
                    a[k, i] = 2 * (1 / h(i) + 1 / (h(i + 1)));
                    a[k, i + 1] = 1 / h(i + 1);
                    b[k] = 3 * ((nodes[i, 1] - nodes[i - 1, 1]) / (h(i) * h(i)) + (nodes[i + 1, 1] - nodes[i, 1]) / (h(i + 1) * h(i + 1)));
                }
                else if (k == (n - 2)) // (11.71)
                {
                    a[i - 1, 0] = -4 / h(1);
                    a[i - 1, 1] = -2 / h(1);
                    b[i - 1] = dd[0, 1] - 6 * (nodes[1, 1] - nodes[0, 1]) / (h(1) * h(1));
                }
                else if (k == (n - 1)) // (11.72)
                {
                    a[i - 1, (int)dd[1, 0] - 1] = 2 / h((int)dd[1, 0]);
                    a[i - 1, (int)dd[1, 0]] = 4 / h((int)dd[1, 0]);
                    b[i - 1] = dd[1, 1] + 6 * (nodes[(int)dd[1, 0], 1] - nodes[(int)dd[1, 0] - 1, 1]) / (h((int)dd[1, 0]) * h((int)dd[1, 0]));
                }
            }
            alglib.rmatrixsolve(a, 4, b, out info, out rep, out s); // Вирішення отриманої системи рівнянь
        }
        public void PrintMatrix() // Виводить матрицю системи рівнянь на екран
        {
            for (int i = 0; i < n; i++)
            {
                for (int k = 0; k < n; k++)
                {
                    Console.Write("{0:0.000}\t", a[i, k]);
                }
                Console.WriteLine("| {0:0.000}", b[i]);
            }
        }
        public void PrintCoeffs() // Виводить на екран коефіцієнти інтерполяційного сплайну
        {
            for (int i = 0; i < s.Length; i++)
                Console.WriteLine("s_{0} : {1:0.000}",i , s[i]);
        }
        private double h(int i) // Повертає значення h, після формули (11.64)
        {
            return nodes[i, 0] - nodes[i - 1, 0];
        }
        public double intFun(double x) // Повертає значення функції, що отримане за допомогою інтерполяційного сплайну
        {
            int i = 1;
            double result = 0;
            if (x >= nodes[0, 0] & x <= nodes[n - 1, 0]) // Якщо точка між першим та останнім вузлом
            {
                for (int k = 0; k < (n - 1); k++) // Пошук сусідніх з точкою вузлів
                {
                    if (x >= nodes[k, 0] & x <= nodes[k + 1, 0]) // Сусідні вузли знайдено
                    {
                        i = k + 1;
                        break;
                    }
                }
            }
            else if (x <= nodes[0, 0]) // Якщо точка перед першим вузлом
                i = 1;
            else // Якщо точка після останнього вузла
                i = n - 1;

            /* Обчислення за формулою 11.64 */

            result += ((x - nodes[i, 0]) * (x - nodes[i, 0]) * (2 * (x - nodes[i - 1, 0]) + h(i))) * nodes[i - 1, 1] / (h(i) * h(i) * h(i));
            result += ((x - nodes[i - 1, 0]) * (x - nodes[i - 1, 0]) * (2 * (nodes[i, 0] - x) + h(i)) * nodes[i, 1]) / (h(i) * h(i) * h(i));
            result += (x - nodes[i, 0]) * (x - nodes[i, 0]) * (x - nodes[i - 1, 0]) * s[i - 1] / (h(i) * h(i));
            result += (x - nodes[i - 1, 0]) * (x - nodes[i - 1, 0]) * (x - nodes[i, 0]) * s[i] / (h(i) * h(i));
            return result;
        }
    }
}
