#include <iostream>
#include <vector>
#include <math.h>

// Auswertung von p_{i,k}(t) mit Neville Schema bei t=x
template<typename REAL>
REAL evaluateNeville(REAL x, std::vector<REAL> &xs, std::vector<REAL> &ys, int i, int k) {
    if(k == 0) {
        return ys[i];
    } else {
        // nevile rekursionsformel
        return ((x - xs[i])*evaluateNeville(x, xs, ys,i+1,k-1) - (x - xs[i+k])*evaluateNeville(x, xs, ys, i, k-1))/(xs[i+k]-xs[i]);
    }
}

// Auswertung eines Interpolationspolynoms p(t) bei t=x
template<typename REAL>
REAL evaluate(REAL x, std::vector<REAL> &xs, std::vector<REAL> &ys) {
    // nutze neville verfahren mit p(x) = p_{0,n}(x)
    evaluateNeville(x, xs, ys, 0, xs.size()-1);
}

// Auswertung eines Interpolationspolynoms einer Funktion bei vorgegebenen Stuetzstellen
template <typename REAL>
REAL evaluateFunction(REAL x, std::vector<REAL> &xs, REAL(*f)(REAL)) {
    std::vector<double> ys(xs.size());
    for (int i = 0; i < xs.size(); i++) {
        // berechne stuetzstellen
        ys[i] = f(xs[i]);
    }
    // verwende polynominterpolation
    return evaluate(x, xs, ys);
}

// Auswertung eines Interpolationspolynoms einer Funktion bei äquidistanten Stützstellen
template <typename REAL>
REAL evaluateFunctionEqualDist(REAL x, REAL a, REAL b, REAL h, REAL(*f)(REAL)) {
    std::vector<double> xs(std::floor((b-a)/h));
    std::vector<double> ys(xs.size());
    for (int i = 0; i < xs.size(); i++) {
        // berechne stuetzstellen
        xs[i] = a + i*h;
        // und funktionswerte
        ys[i] = f(xs[i]);
    }
    // werte polynom aus
    return evaluate(x, xs, ys);
}

// beispiel funktionen
template <typename REAL>
REAL f1(REAL x) {
    return 1/(1+std::pow(x,2));
}

template <typename REAL>
REAL f2(REAL x) {
    return std::sqrt(std::abs(x));
}

int main() {
    // polynom grad
    int n = 5;
    // abstand zwischen stützstellen
    double h = 2.0 / n;

    // schrittweite
    double dx = 2.0 / 1000;
    double x = 0;
    for(int i = 0; i < 1000; i++) {
        x = -1.0 + i*dx;
        // werte polynom mit äquidistanten stuetzstellen zwischen -1 und 1 mit abstand h bei t=x aus
        std::cout << x << " " << evaluateFunctionEqualDist(x, -1.0, 1.0, h, f2) << std::endl;
    }
}