#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <cstring>
#include <fstream>
#include <iomanip>
#include <ios>
#include <iostream>
#include <numeric>
#include <string>
#include <utility>
#include <vector>

using namespace std::literals;

using points_t = std::pair<std::vector<double>, std::vector<double>>;
using matrix_t = std::vector<std::vector<double>>;

std::vector<double> newton_differences(points_t points) {
        std::vector<double> factors{};
        auto &[x, fx] = points;
        int n = points.first.size();
        for (int i = 1; i <= n; ++i) {
                factors.push_back(fx[0]);
                for (int j = 0; j < (n - i); ++j) {
                        fx[j] = (fx[j + 1] - fx[j]) / (x[j + i] - x[j]);
                }

        }
        return factors;
}

double newton_polynomial(points_t points, std::vector<double> factors, double x) {
        auto xs = points.first;
        int n = points.second.size() - 1;
        double val = 0;
        for (int k = 0; k <= n; ++k) {
                double acc = 1;
                for(int i = 0; i < k; ++i) {
                        acc *= (x - xs[i]);
                }
                val += acc * factors[k];
        }
        return val;
}

void calculate_natural_cubic_spline_matrix(points_t points, matrix_t &mat) {
        auto &[xs, fx] = points;
        int n = xs.size();

        // Construção da matriz
        for (int i = 1; i < n - 1; ++i) {
                mat[i][i - 1] = (xs[i] - xs[i - 1])/6;
                mat[i][i]     = (xs[i + 1] - xs[i - 1])/3;
                mat[i][i + 1] = (xs[i + 1] - xs[i])/6;
                mat[i][n]     = (fx[i + 1] - fx[i])/(xs[i + 1] - xs[i]) - (fx[i] - fx[i - 1])/(xs[i] - xs[i - 1]);
        }
        mat[0][0] = 1;
        mat[0][n - 1] = 0;
        mat[n - 1][n - 1] = 1;
        mat[n - 1][n] = 0;

        std::ofstream matriz("matriz.txt", std::ios_base::app);

        matriz << "mat:\n";
        for (auto i : mat){
	        for (auto j = i.begin(); j != i.end() - 1; ++ j)
		        matriz << *j << ",";
	        matriz << "\n";
        }

        matriz << "b:\n";
        for (auto i : mat){
	        matriz << i[n] << '\n';
        }

        // Passar para a forma triangular
        for (int k = 0; k < n; ++k) {
                for (int i = k + 1; i < n; ++i) {
                        if (mat[k][k] != 0) {
                                double mul = mat[i][k]/mat[k][k];
                                for (int j = k; j < n; ++j) {
                                        mat[i][j] -= mul * mat[k][j];
                                }
                                mat[i][n] -= mul * mat[k][n];
                        }
                }
        }

        // Resolução da matriz
        for (int i = n - 1; i > 0; --i) {
                if (mat[i][i] != 0) {
                        double mul = mat[i - 1][i]/mat[i][i];
                        for (int j = 1; j < n + 1; ++j) {
                                mat[i - 1][j] -= mul * mat[i][j];
                        }
                        mat[i - 1][i] = 0;
                        mat[i][n] /= mat[i][i];
                        mat[i][i] = 1;
                }
        }

        matriz << "x:\n";
        for (auto i : mat){
	        matriz << i[n] << '\n';
        }
}

double natural_cubic_spline(points_t points, matrix_t &mat, double x) {

        auto &[xs, fx] = points;
        int n = xs.size();

        int i = 0;
        for (int i_ = 0; i_ < n; ++i_) {
                if (xs[i_] > x){
                        i = i_;
                        break;
                }
        }

        double hi = xs[i] - xs[i - 1];
        return mat[i - 1][n] * std::pow((xs[i] - x), 3)/(6 * hi) +
               mat[i][n] * std::pow((x - xs[i - 1]), 3)/(6 * hi) +
               (fx[i - 1] - mat[i - 1][n] * (hi * hi)/6)*(xs[i] - x)/hi +
               (fx[i] - mat[i][n] * (hi * hi)/6)*(x - xs[i - 1])/hi;

}

void exercise_a(points_t points) {
        auto &[xs, fx] = points;
        auto factors = newton_differences(points);
        std::ofstream poly{"a_polynomial.txt"};
        for(double x = 0; x < 4; x += 0.001)
                poly << x << " " << newton_polynomial(points, factors, x) << '\n';

        std::ofstream spline{"a_spline.txt"};
        unsigned long n = xs.size();
        matrix_t mat(n, std::vector<double>(n + 1));
        calculate_natural_cubic_spline_matrix(points, mat);
        for(double x = 0; x < 4; x += 0.001)
	        spline << x << " " << natural_cubic_spline(points, mat, x) << '\n';
}

void exercise_b() {
        points_t points;
        auto &[xs, fx] = points;
        auto f = [](double x) { return 4 * std::pow(x, 2) + std::sin(9 * x); };

        for (double x = -1; x <= 1; x += (1 - -1)/8.0) {
                xs.push_back(x);
                fx.push_back(f(x));
        }

        std::ofstream poly{"b_polynomial.txt"};
        auto factors = newton_differences(points);
        for(double x = -1; x < 1; x += 0.001)
                poly << x << " " << newton_polynomial(points, factors, x) << '\n';

        std::ofstream spline{"b_spline.txt"};
        unsigned long n = xs.size();
        matrix_t mat(n, std::vector<double>(n + 1));
        calculate_natural_cubic_spline_matrix(points, mat);
        for(double x = -1; x < 1; x += 0.001)
                spline << x << " " << natural_cubic_spline(points, mat, x) << '\n';

        std::cout << "Erros poly: "
                  << std::abs(std::round((f(0.3) - newton_polynomial(points, factors, 0.3)) * 100.0) / 100.0) << " "
                  << std::abs(std::round((f(0.83) - newton_polynomial(points, factors, 0.83)) * 100.0) / 100.0) << '\n';

        std::cout << "Erros spline: "
                  << std::abs(std::round((f(0.3) - natural_cubic_spline(points, mat, 0.3))*100.0)/100.0) << " "
                  << std::abs(std::round((f(0.83) - natural_cubic_spline(points, mat,0.83))*100.0)/100.0) << '\n';

        std::cout << "p(0.3) = " << newton_polynomial(points, factors, 0.3) << ' '
                  << "p(0.83) = " << newton_polynomial(points, factors, 0.83) << '\n';

        std::cout << "p(0.3) = " << natural_cubic_spline(points, mat, 0.3) << ' '
                  << "p(0.83) = " << natural_cubic_spline(points, mat, 0.83) << '\n';
}

int main(int argc, char **argv) {
        std::istream *s;
        std::ifstream file;
        if (argc == 2) {
                if (argv[1] == "-"s) {
                        s = &std::cin;
                } else {
                        file.open(argv[1]);
                        s = &file;
                }
        } else{
                std::cerr << "An argument required\n";
                return 1;
        }

        points_t points;
        std::pair<double, double> p;
        while(true) {
                *s >> p.first >> p.second;
                if(s->eof())
                        break;

                points.first.push_back(p.first);
                points.second.push_back(p.second);
        }

        exercise_a(points);
        exercise_b();

}