NumericalAnalysis_LAB_1

import java.io.BufferedReader;
import java.io.IOException;
import java.io.InputStreamReader;
import java.text.DecimalFormat;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collections;


// Xj: {-2 0 2 3 4}
// Yj: {3 4 1 2 -1}


public class Main {
    public static void main(String[] args) throws IOException {

        // Поток для чтения:
        BufferedReader reader = new BufferedReader(new InputStreamReader(System.in));


        // Считываем значения всех Xj-ых с клавиатуры:
        System.out.println("Введите значения Xj через пробел:");
        String[] string_xjValues = reader.readLine().split("\\s+");


        // Количество узлов интерполяции:
        int n = string_xjValues.length;


        // Занесем значения всех Xj-ых в список xjValues типа X:
        ArrayList<X> xjValues = new ArrayList<>();
        for (String string_xjValue : string_xjValues) {
            double coefficient = Double.parseDouble(string_xjValue);
            X newElemX = new X(coefficient, 0);
            xjValues.add(newElemX);
        }


        // Считываем значения всех Yj-ых с клавиатуры:
        System.out.println("Введите значения Yj через пробел:");
        String[] string_yjValues = reader.readLine().split("\\s+");


        // Преобразуем значения Yj-ых в целочисленный тип:
        double[] yjValues = new double[n];
        for (int i = 0; i < n; i++) {
            yjValues[i] = Double.parseDouble(string_yjValues[i]);
        }



        // (1). Построим полином Лагранжа:
        ArrayList<X> lagrangePolynomial = X.buildLagrangePolynomial(xjValues, yjValues, n);

        // Выведем его на экран:
        System.out.println("\nПолином Лагранжа:\n" + lagrangePolynomial + "\n");




        // (2). Построим полином Ньютона:
        ArrayList<X> newtonPolynomial = X.buildNewtonPolynomial(xjValues, yjValues, n);

        // Выведем его на экран:
        System.out.println("Полином Ньютона:\n" + newtonPolynomial);




        // (3а). Для построения полинома методом найменьших квадратов (Least Square Method, LSM),
        // сначала запросим степень полинома - m (1 <= m < n):
        int m = -1;
        do {
            System.out.print("\nВведите целое число m (Степень полинома Qm(x)), 0 <= m <= " + (n - 1) + ": ");
            m = Integer.parseInt(reader.readLine());
        } while (m < 0 || m > n - 1);

        // Построим полином методом найменьших квадратов:
        ArrayList<X> lsmPolynomial = X.buildLSMPolynomial(xjValues, yjValues, n, m);

        // Выведем его на экран:
        System.out.println("\nПолином по методу найменьших квадратов степени m = " + m + " :\n" + lsmPolynomial);

        // (3б). Теперь найдём среднеквадратическое отклонение:
        double standardDeviationValue = X.standardDeviationLSM(lsmPolynomial, xjValues, yjValues, n);
        System.out.println("\nСреднеквадратическое отклонение = " + standardDeviationValue);



        // (4). Задача экстраполирования.
        // Найдём значения полиномов в точках X(i) = (x(i) + (h(i) / 2)),
        // где h(i) = x(i+1) - x(i), при этом i є [0; n-1].
        // Также найдем значения в точках:
        // X(-1) = x(0) - h(0)/2 и X(n) = x(n) + h(n-1)/2.

        // (4.а):
        // Найдём сначала нужные нам точки:
        double[] extrapolateNodes = X.extrapolateNodes(xjValues, n);
        System.out.println("\nТеперь проверим значения в точках: " + Arrays.toString(extrapolateNodes) + "\n");

        // (4.б):
        // Значения в этих точках для полинома Лагранжа:
        double[] lagrangeExtrapolateValues = X.extrapolatePolynomialValues(lagrangePolynomial, extrapolateNodes);
        for (int i = 0; i < lagrangeExtrapolateValues.length; i++) {
            //L(Xi) = lagrangeExtrapolateValues[i]
            System.out.println("L(" + extrapolateNodes[i] + ") = " + new DecimalFormat("#.###").format(lagrangeExtrapolateValues[i]));
        }
        System.out.println();

        // Значения в этих точках для полинома Ньютона:
        double[] newtonExtrapolateValues = X.extrapolatePolynomialValues(newtonPolynomial, extrapolateNodes);
        for (int i = 0; i < newtonExtrapolateValues.length; i++) {
            //N(Xi) = newtonExtrapolateValues[i]
            System.out.println("N(" + extrapolateNodes[i] + ") = " + new DecimalFormat("#.###").format(newtonExtrapolateValues[i]));
        }
        System.out.println();

        // Значения в этих точках для полинома LSM:
        double[] lsmExtrapolateValues = X.extrapolatePolynomialValues(lsmPolynomial, extrapolateNodes);
        for (int i = 0; i < lsmExtrapolateValues.length; i++) {
            //Qm(Xi) = lsmExtrapolateValues[i]
            System.out.println("Qm(" + extrapolateNodes[i] + ") = " + new DecimalFormat("#.###").format(lsmExtrapolateValues[i]));
        }


        // Закрыть поток для чтения:
        reader.close();
    }
}

public class X {
    private double coefficient;
    private final int degree;

    //Конструктор для создания переменной x с конкретным коэф-м и показателем степени.
    public X (double coefficient, int degree) {
        this.coefficient = coefficient;
        this.degree = degree;
    }

    //Конструктор для создания обычной переменной x.
    public X () {
        this.coefficient = 1;
        this.degree = 1;
    }

    //(1). Функция, которая создаёт интерполяционный полином Лагранжа.
    public static ArrayList<X> buildLagrangePolynomial(ArrayList<X> xjValues, double[] yjValues, int n) {

        //Пока что пустой интерполяционный полином Лагранжа:
        ArrayList<X> lagrangePolynomial = new ArrayList<>();

        //Суммируем все j-е базисные полиномы, умноженные на сооотв. значения Yj-x,
        //получая, таким образом, полином Лагранжа.
        for (int j = 0; j < n; j++) {

            double denominator = 1; //знаменатель (число)
            ArrayList<X> nominator = new ArrayList<>(Collections.singletonList(new X(1, 0))); //числитель (полином)


            //Перемножение всех числителей и знаменателей от 0 до n, не включая индекс j:
            for (int m = 0; m < n; m++) {
                if (m != j) {

                    //Произведение числителей (полином n-1 степени)
                    ArrayList<X> secondBraces = new ArrayList<>(Arrays.asList(new X(), new X((-1) * xjValues.get(m).coefficient, xjValues.get(m).degree)));
                    nominator = multiplyBraces(nominator, secondBraces);

                    //Произведение знаменателей (число):
                    denominator *= (xjValues.get(j).coefficient - xjValues.get(m).coefficient);
                }
            }

            //Базисный полином Лагранжа - произведение числителя(полинома) на перевернутый знаменатель(число).
            ArrayList<X> basisPolynomial = X.multiplyNumberAndBrace((1.0 / denominator), nominator);

            //Перемножение Yj на базисный полином Лагранжа Lj(x):
            ArrayList<X> multipliedPolynomial = X.multiplyNumberAndBrace(yjValues[j], basisPolynomial);

            //Добавляем результат перемножения в выходной полином Лагранжа:
            lagrangePolynomial.addAll(multipliedPolynomial);
        }

        //Приведём абсолютно все свободные слагаемые в полиноме Лагранжа:
        cancellation(lagrangePolynomial);

        //Возвращаем полином Лагранжа:
        return lagrangePolynomial;
    }

    //(2). Функция, которая создаёт интерполяционный полином Ньютона.
    public static ArrayList<X>  buildNewtonPolynomial(ArrayList<X> xjValues, double[] yjValues, int n) {
        /*
        * Полином Ньютона находится по формуле:
        * Pn(x) = ((сумма от k=1 до k=n) произведений f(x0,...,xk)
        * на (произведения от i=0 до i=k-1) (x - xi)), где f(x0) = y0 по условию,
        * а f(x0,...,xk) - разделенная разность k-го порядка,
        * которая представляет собой число.
        */

        //Пока что пустой интерполяционный полином Ньютона:
        ArrayList<X> newtonPolynomial = new ArrayList<>();

        //Теперь рассчитаем оставшуюся сумму:
        for (int k = 1; k <= n; k++) {

            double dividedDifference = 0; // разделенная разность
            ArrayList<X> polynomial = new ArrayList<>(Collections.singletonList(new X(1, 0))); // пока что пустой полином произведения (X - Xi)

            //Поиск полинома в виде произведения (X - Xi)
            for (int i = 0; i < k - 1; i++) {
                ArrayList<X> tempPolynomial = new ArrayList<>(Arrays.asList(new X(), new X((-1) * xjValues.get(i).coefficient, xjValues.get(i).degree)));
                polynomial = multiplyBraces(polynomial, tempPolynomial);
            }

            //Поиск разделенной разности: Сумма (f(Xi) / Произведение(xi - xj))
            for (int i = 0; i < k; i++) {

                double f_xi = yjValues[i]; // Числитель = f(Xi) = Yi
                double denominator = 1; // Знаменатель = Произведение (xi - xj)

                for (int j = 0; j < k; j++) {
                    if (j != i) {
                        denominator *= (xjValues.get(i).coefficient - xjValues.get(j).coefficient);//(xi - xj)
                    }
                }

                dividedDifference += (f_xi * (1.0/denominator));
            }

            //Умножаем разделенную разность(число) на наш временный полином.
            ArrayList<X> multipliedPolynomial = multiplyNumberAndBrace(dividedDifference, polynomial);

            //Добавляем каждое слагаемое (полином (n-1)-й степени) искомой суммы в наш полином Ньютона.
            newtonPolynomial.addAll(multipliedPolynomial);

        }

        //Приведём абсолютно все свободные слагаемые в полиноме Ньютона:
        cancellation(newtonPolynomial);

        //Развернём его в обратном порядке:
        Collections.reverse(newtonPolynomial);

        //Возвращаем полином Ньютона:
        return newtonPolynomial;
    }

    //(3а). Функция, которая строит интерп. полином методом найменьших квадратов:
    public static ArrayList<X>  buildLSMPolynomial(ArrayList<X> xjValues, double[] yjValues, int n, int m) {
        //Создадим полином
        ArrayList<X> lsmPolynomial = new ArrayList<>();

        //Для его нахождения создаём матрицу
        double[][] matrix = new double[m+1][m+2];

        //Заполним матрицу, не включая столбец свободных членов, соотв. значениями Ck,
        //которые наход. по формуле:
        //Ck = Сумма от p до n (Xp ^ k), где p є [0; n], k є [0; 2m]
        double ck = 0;
        for (int i = 0; i < m + 1; i++) {
            for (int j = 0, k = i;  j < (m + 1) && k < (m + 1 + i);   j++, k++) {

                //Находим Ck-й коэффициент:
                for (int p = 0; p < n; p++) {
                    ck += Math.pow(xjValues.get(p).coefficient, k);
                }

                matrix[i][j] = ck;
                ck = 0;
            }
        }

        //Заполним столбец свободных членов коэф-ми Di,
        //которые рассчит-ся по формуле:
        //Di = Сумма от j до n (f(Xj) * (Xj ^ i)), где j є [0; n], i є [0; m]
        double di = 0;
        for (int i = 0; i < m + 1; i++) {
            for (int j = 0; j < n; j++) {
                di += (yjValues[j] * Math.pow(xjValues.get(j).coefficient, i));
            }

            matrix[i][m+1] = di;
            di = 0;
        }

        //Совершаем прямой ход по Гауссу и приводим матрицу к ступенчатому виду
        gaussForward(matrix);

        //Найдём коэф-ы ai для нашего искомого полинома обратным ходом по Гауссу:
        double[] aiCoefficients = gaussMethod(matrix);

        //Создадим полином М-й степени по полученным коэф-м ai:
        for (int i = 0; i < aiCoefficients.length; i++) {
            lsmPolynomial.add(new X(aiCoefficients[i], i));
        }

        //Развернем полином в обратном порядке для читабельности.
        Collections.reverse(lsmPolynomial);

        //Вернём полученный полином.
        return lsmPolynomial;
    }

    //(3б). Функция, считающая среднеквадратическое отклонение для полинома LSM:
    public static double standardDeviationLSM(ArrayList<X> lsmPolynomial, ArrayList<X> xjValues, double[] yjValues, int n) {

        double stdDeviation = 0;
        double tempSum;

        for (int j = 0; j < n; j++) {
            tempSum = getPolynomialValue(xjValues.get(j).coefficient, lsmPolynomial);
            stdDeviation += Math.pow((yjValues[j] - tempSum), 2);
        }

        stdDeviation = Math.sqrt(stdDeviation / n);
        return stdDeviation;
    }

    //(4а). Функция, которая находит нужные узлы экстраполяции на оси ОХ:
    public static double[] extrapolateNodes(ArrayList<X> xjValues, int n) {
        double[] extrapolateNodes = new double[n + 1];

        // X(-1) = x(0) - h(0)/2.
        extrapolateNodes[0] = xjValues.get(0).coefficient - ((xjValues.get(1).coefficient - xjValues.get(0).coefficient) / 2);

        // X(i) = (x(i) + (h(i) / 2)),
        // где h(i) = x(i+1) - x(i), при этом i є [0; n-1].
        for (int i = 0; i < n - 1; i++) {
            extrapolateNodes[i+1] = xjValues.get(i).coefficient + ((xjValues.get(i+1).coefficient - xjValues.get(i).coefficient) / 2);;
        }

        // X(n) = x(n) + h(n-1)/2.
        extrapolateNodes[n] = xjValues.get(n-1).coefficient + ((xjValues.get(n-1).coefficient - xjValues.get(n-2).coefficient) / 2);

        return extrapolateNodes;
    }

    //(4б). Функция, которая находит значения полинома во всех заданных узлах экстраполяции:
    public static double[] extrapolatePolynomialValues(ArrayList<X> polynomial, double[] extrapolateNodes) {
        double[] polynomialValues = new double[extrapolateNodes.length];

        for (int i = 0; i < extrapolateNodes.length; i++) {
            polynomialValues[i] = getPolynomialValue(extrapolateNodes[i], polynomial);
        }

        return polynomialValues;
    };

    //Получим значение полинома в заданной точке.
    private static double getPolynomialValue(double xj, ArrayList<X> polynomial) {
        double value = 0;

        for (X elem: polynomial) {
            value += elem.coefficient * Math.pow(xj, elem.degree);
        }

        return value;
    }

    //Перемножение двух полиномов(скобок):
    private static ArrayList<X> multiplyBraces(ArrayList<X> firstBraces, ArrayList<X> secondBraces) {

        //Результат перемножения 2-х скобок запишем в этот список:
        ArrayList<X> resultBraces = new ArrayList<>();

        //Перемножение двух скобок (коэффициенты перемножаем, степени складываем):
        for(X firstElem: firstBraces) {
            for (X secondElem: secondBraces) {
                double resultCoefficient = firstElem.coefficient * secondElem.coefficient;
                int resultDegree = firstElem.degree + secondElem.degree;
                X resultElemX = new X(resultCoefficient, resultDegree);
                resultBraces.add(resultElemX);
            }
        }

        //Приводим подобные слагаемые:
        cancellation(resultBraces);

        //Возвращаем полученное произведение в виде полинома(список из объектов типа X)
        return resultBraces;
    }

    //Умножение числа на полином(скобку):
    private static ArrayList<X> multiplyNumberAndBrace(double number, ArrayList<X> braces) {
        for (X elem : braces) {
            elem.coefficient *= number;
        }
        return braces;
    }

    //Приведение подобных слагаемых:
    private static void cancellation(ArrayList<X> braces) {
        for (int i = 0; i < braces.size(); i++) {
            for (int k = 0; k < braces.size(); k++) {
                if (
                        (k != i)
                        && (braces.get(k).degree == braces.get(i).degree)
                        && (braces.get(i).coefficient != 0)
                ) {
                    braces.get(i).coefficient += braces.get(k).coefficient;
                    braces.get(k).coefficient = 0;
                }
            }
        }
        //Уберём все элементы, у которых коэффициент обнулен.
        braces.removeIf(elem -> (Math.abs(elem.coefficient) <= 0.00001));
    }

    //Функция сдвига заданной строки в конец.
    private static void changeRows(double[][] matrix, int position) {
        int n = matrix.length;
        double[] temp = matrix[position];
        for (int i = position; i < n - 1; i++) {
            matrix[i] = matrix[i + 1]; //каждой строке, начиная с заданной, присвоим строку, что стоит после неё.
        }
        matrix[n-1] = temp;
    }

    //Функция прямого хода по Гауссу.
    private static void gaussForward(double[][] matrix) {
        int n = matrix.length;
        int m = matrix[0].length;

        double delta;
        for (int k = 0; k < n; k++) {
            for (int i = k + 1; i < n; i++) {
                while(matrix[k][k] == 0) {
                    changeRows(matrix, k);
                }
                delta = matrix[i][k] / matrix[k][k]; // формула (1)
                for (int j = k; j < m; j++) {
                    matrix[i][j] -= delta * matrix[k][j]; // формула (2)
                }
            }
        }
    }

    //Метод Гаусса обратный.
    private static double[] gaussMethod(double[][] matrix) {
        int n = matrix.length;
        int m = matrix[0].length;
        double[] result = new double[m-1];

        double[] b = new double[n];//Столбец свободных членов
        for (int i = 0; i < n; i++) {
            b[i] = matrix[i][m - 1];
        }

        for (int i = n - 1; i >= 0; i--) {
            double sum = 0.0;
            for (int j = i + 1; j < n; j++) {
                sum += matrix[i][j] * result[j];
            }
            result[i] = (b[i] - sum) / matrix[i][i];
        }
        return result;
    }

    @Override
    public String toString() {

        if (this.degree == 0) //Если показатель степени при Х == 0, то выводим лишь коэффициент.
            //Форматирование коэф. до максимум 5 знаков после запятой
            return "" + new DecimalFormat("#.#####").format(this.coefficient);

        /*
        else if(this.coefficient > 0) //Если коэф. больше нуля, то выводим с дополнительным знаком "+".
            //Форматирование коэф. до максимум 5 знаков после запятой
            return "+" + new DecimalFormat("#.#####").format(this.coefficient) + "x^" + this.degree;
        */

        else //В других случаях выводим всё без изменений.
            //Форматирование коэф. до максимум 5 знаков после запятой
            return new DecimalFormat("#.#####").format(this.coefficient) + "x^" + this.degree;
    }
}