#include <iostream>
using namespace std;
double x[115];
void Gauss(double** A, int n)
{
    for (int i = 0; i < n; i++)
        x[i] = 0;

    for (int j = 0; j < n; j++) {
        for (int i = 0; i < n; i++) {
            if (i>j) {
                double c = A[i][j] / A[j][j];
                for (int k = 0; k < n + 1; k++) {
                    A[i][k] = A[i][k] - c*A[j][k];
                }
            }
        }
    }
    x[n - 1] = A[n - 1][n] / A[n - 1][n - 1];
    for (int i = n - 1; i >= 1; i--) {
        double sum = 0;
        for (int j = i + 1; j <= n; j++) {
            sum = sum + A[i - 1][j - 1] * x[j - 1];
        }
        x[i - 1] = (A[i - 1][n] - sum) / A[i - 1][i - 1];
    }
}
double Wielomian(double a, double b, double c, double d,double S, int i, int staly)
{
    S = a + b*(x[i] - x[staly]) + c*(x[i] - x[staly])*(x[i] - x[staly]) + d*(x[i] - x[staly])*(x[i] - x[staly])*(x[i] - x[staly]);
    return S;
}
int main()
{
    double freq[20] = { 2.160913 ,  2.184642 , 2.208656 , 2.232956,  2.257543 , 2.282417 , 2.307579 , 2.333029, 2.358767 , 2.384794,  2.411110  ,2.437714 , 2.464608 , 2.491789 , 2.519259 , 2.547017 , 2.575062 , 2.603393 , 2.632010 , 2.660913 };
    double S21[20] = { 0.0385736053 ,  0.0418602337 , 0.0406238069,  0.0222182757 , 0.0681165327,  0.5274988276 , 0.9992090072 , 0.9877714749 , 0.9997958915,  0.9872738013, 0.9872738013,  0.9997958915 , 0.9877714749 , 0.9992090072,  0.5274988276, 0.0681165327,  0.0222182757 , 0.0406238069 , 0.0418602337, 0.0385736053 };
    double h[20];
    //obliczanie h
    for (int i = 1; i < 20; i++) h[i] = freq[i] - freq[i - 1];
    h[0] = 0.0;

    double mi[20];
    double lambda[20];
    double delta[20];
    //obliczanie mi , lambda, delta
    for (int i = 1; i < 20; i++)
    {
        mi[i - 1] = h[i - 1] / (h[i] + h[i - 1]);
        lambda[i - 1] = lambda[i - 1] / (lambda[i - 1] + lambda[i]);
        delta[i - 1] = 6 / (h[i - 1] + h[i])*(((S21[i] - S21[i]) / h[i]) - ((S21[i - 1] - S21[i - 2]) / h[i - 1]));
    }
    //warunki brzegowe
    mi[19] = 0;
    delta[0] = 0;
    delta[19] = 0;
    lambda[0] = 0;
    //wypelnianie macierzy
    double macierz[20][21];
    for (int i = 0; i < 20; i++) for (int j = 0; j < 21; j++) macierz[i][j] = 0.0;
    for (int i = 0; i < 20; i++) macierz[i][i] = 2.0;
    for (int i = 0; i < 19; i++) macierz[i][i + 1] = lambda[i];
    for (int i = 0; i < 20; i++) macierz[i + 1][i] = mi[i];
    for (int i = 0; i < 20; i++) macierz[i][20] = delta[i];

    //przepisywanie macierzy
    double** macierz2 = new double*[20];
    for (int i = 0; i < 20; i++)
        macierz2[i] = new double[21];
    for (int i = 0; i < 20; i++)
        for (int j = 0; j < 21; j++)
            macierz2[i][j] = macierz[i][j];
    //odpalamy algorytm Gaussas
    Gauss(macierz2,20);

    double a[20], b[20], c[20] , d[20], S[115];
    //obliczanie wspolczynnikow a,b,c,d
    for (int i = 1; i < 20; i++)
    {
        a[i - 1] = S21[i - 1];
        b[i - 1] = (((S21[i] - S21[i - 1]) / h[i]) - ((2 * x[i - 1] + x[i])*h[i] / 6));
        c[i - 1] = x[i - 1] / 2;
        d[i - 1] = (x[i] - x[i - 1]) / (6 * h[i]);
    }
    //zamiana indeksow na nowe, przygotowane do rozszerzenia
    double xstare[20];
    for (int i = 0; i < 20; i++) xstare[i] = x[i];
    for (int i = 0; i < 115; i++) x[i] = 0;
    for (int i = 0; i < 20; i++) x[i * 6] = xstare[i];
                                //przepisane wyniki gaussa na odpowiedznie indeksy , co 6
                                //tablica x jest gotowa do wypełniania na miejscach od 0 do 19 stare iksy,
                                //0 , 6 , 12 , 18 , 24 , 30 , 36 , 42 , 48...





        for (int j = 0; j < 20; j++)
            delete[] macierz2[j];
        delete[] macierz2;

    return 0;
}