---

#include <iostream>
#include <vector>
#include <cmath>
#include <iomanip>
#include <algorithm>

using namespace std;

const string RESET = "\033[0m";
const string GREEN = "\033[32m";
const string RED = "\033[31m";
const string YELLOW = "\033[33m";

double targetFunction(double x) {
    return 2.1 * sin(0.37 * x);
}

double getM3() {
    return 2.1 * pow(0.37, 3);
}

vector<double> generateYNodes(const vector<double>& X) {
    vector<double> Y;
    for (double x_val : X) {
        Y.push_back(targetFunction(x_val));
    }
    return Y;
}

double getTargetPoint() {
    double x_target;
    cout << "Введите значение точки x(0<=x<=3): ";
    cin >> x_target;
    return x_target;
}

vector<int> findNearestIndices(const vector<double>& X, double x_target) {
    int n = X.size();
    int i0, i1, i2;

    if (x_target <= X[1]) {
        i0 = 0; i1 = 1; i2 = 2;
    }

    else if (x_target >= X[n - 2]) {
        i2 = n - 1; i1 = i2 - 1; i0 = i2 - 2;
    }

    else {

        int interval_idx = 1;
        for (int i = 1; i < n - 2; ++i) {
            if (x_target >= X[i] && x_target <= X[i + 1]) {
                interval_idx = i;
                break;
            }
        }

        i0 = interval_idx - 1;
        i1 = interval_idx;
        i2 = interval_idx + 1;
        
    }

    return { i0, i1, i2 };
}

double calculateLagrange(double x_target, const vector<double>& X, const vector<double>& Y, const vector<int>& idx) {
    int i0 = idx[0]; int i1 = idx[1]; int i2 = idx[2];

    double l0 = ((x_target - X[i1]) * (x_target - X[i2])) / ((X[i0] - X[i1]) * (X[i0] - X[i2]));
    double l1 = ((x_target - X[i0]) * (x_target - X[i2])) / ((X[i1] - X[i0]) * (X[i1] - X[i2]));
    double l2 = ((x_target - X[i0]) * (x_target - X[i1])) / ((X[i2] - X[i0]) * (X[i2] - X[i1]));

    return Y[i0] * l0 + Y[i1] * l1 + Y[i2] * l2;
}

void printResults(double x_target, double y_approx, const vector<int>& idx, const vector<double>& X) {
    double y_exact = targetFunction(x_target);

    double actual_error = abs(y_exact - y_approx);

    double omega = abs((x_target - X[idx[0]]) * (x_target - X[idx[1]]) * (x_target - X[idx[2]]));
    double theoretical_error = (getM3() / 6.0) * omega;

    string error_color = (actual_error <= theoretical_error) ? GREEN : RED;

    cout << "\nИспользованы узлы с индексами: " << idx[0] << ", " << idx[1] << ", " << idx[2] << "\n";
    cout << "Узлы X: [" << X[idx[0]] << "; " << X[idx[1]] << "; " << X[idx[2]] << "]\n";
    cout << "------------------------------------------------------\n";
    cout << "Приближенное значение P2(x):       " << y_approx << "\n";
    cout << "Точное значение функции f(x):      " << y_exact << "\n";
    cout << "Фактическая абс. погрешность:      " << error_color << actual_error << RESET << "\n";
    cout << "Теоретическая оценка погрешности:  " << YELLOW << theoretical_error << RESET << "\n";
}

int main() {
    setlocale(LC_ALL, "Russian");
    cout << fixed << setprecision(6);

    cout << "=== Интерполяция методом Лагранжа (Квадратичная) ===\n\n";

    vector<double> X = { 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0 };

    vector<double> Y = generateYNodes(X);
    double x_target = getTargetPoint();
    vector<int> nearest_indices = findNearestIndices(X, x_target);
    double y_approx = calculateLagrange(x_target, X, Y, nearest_indices);

    printResults(x_target, y_approx, nearest_indices, X);

    return 0;
}