// ConsoleApplication10.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include<iostream>
#include<fstream>
#include<string>
#include<cmath>

#define n 21 // rozmiar tabeli - stopien wielomianu + 1, wynika to z przesuniecia w tabeli
#define e 50 // liczba wezlow
#define m 20 // stopien wielomianu
#define a -2.0 // dolna granica przedzialu
#define b 2.0 // gorna granica przedzialu

double x[e];
double y[e];// tablica argumentow
double yy[e];// tablica argumentow

double A[n + 1][n + 1]; // macierz wspolczynnikow A
double B[n+1]; // macierz wyrazow wolnych
double **U; // macierz U
double **L; // macierz L
double *X; // wektor wspolczynnikow A
double *Y;
using namespace std;
double blad;

void GenerujTablice()// generuje tablice z argumentami i wartosciami funkcji aproksymowanej w tych punktach
{       double w = a; // dolna granica przedzialu
    double odstep = (double)(abs(a)+abs(b))/(e-1); // odstep miedzy kolejnymi punktami ->dlugoscprzedzialu przez liczbeelementow - 1
    for(int i = 0; i <e; i++)
    {   x[i] = w;
        y[i] = sin(x[i])*sin(3*x[i]);
        w += odstep;
    }
}

int macierzeUL() // obliczanie macierzy posrednich L i U
{
    U = new double* [n+1];
    for(int i = 0; i<= n; i++)
    {
        U[i] = new double[n+1]; //Alokacja pamieci
        memset( U[i], 0, (n+1) * sizeof(double));
    }  L = new double* [n+1];
    for(int i = 0; i<= n; i++)
    {   L[i] = new double[n+1]; //Alokacja pamieci
    memset( L[i], 0, (n+1) * sizeof(double) );
    }  for(int i = 1; i <= n; i++)
{
    for (int j = 1; j <= n; j++) {
        U[i][j] = A[i][j]; //Przepisujemy zawartosc macierzy A do macierzy U
        L[i][j] = 0;  // zerowanie macierzy L
    }  }  for(int i=1;i<=n;i++)
    {   L[i][i]=1;
    if( U[i][i] == 0 ) //Sprawdzanie czy na głównej przekątnej nie ma 0
        return 2;
    }
    for(int k=1;k<n;k++)
    {   for(int i=k+1;i<=n;i++)
        {    if (U[k][k]==0)
            return 2; // sprawdzenie czy dzielimy przez 0
            L[i][k]=U[i][k]/U[k][k]; // obliczanie macierzy L
            for(int j=k+1;j<=n;j++)  // obliczanie macierzy U
            {     U[i][j]=U[i][j] - ((U[i][k]*U[k][j])/U[k][k]);
            }    U[i][k]=0;
        }
    }
    return 0;
}

void rozwiaz() // rozwiazanie macierzy trojkatynychgornych
{
        X= new double [n+1];
        Y= new double [n+1];
        Y[1]=B[1];
        for( int i= 2; i<=n; i++)
        {
            double suma=0.0;
            for(int k=1; k<=i-1; k++)
                suma+=L[i][k]* Y[k];
        Y[i]=B[i] - suma;
        }
        X[n]= Y[n] / U[n][n];
        for (int i=n-1; i>=1; i--)
        {   double suma=0.0;
        for(int k= i+1; k<=n; k++)
            suma+=U[i][k]*X[k];
        X[i]=(Y[i]- suma) / U[i][i];  }
}

void GenerujMacierz() // funkcja odpowiadajaca za wygenerowanie macierzy współczynników
 {
            for (int j = 0; j <= m; j++)
            for (int k = 0; k <= m; k++)
            {    A[k+1][j+1] = 0;
            for (int i = 0; i <e; i++)
                A[k+1][j+1] += pow(x[i],(k+j));
            }
}

double LiczR(int j) // funkcja wykorzystywana podczas generowanie kolumny wyrazów wolnych B ( w skryptowych wzorach wektor R)
{
        double w = 0.0;
        for(int i = 0; i <e; i++)
            w += y[i]*pow(x[i],j);
        return w;
}

    void GenerujB() // funkcja odpowiadająca za wygenerowanie kolumny wyrazów wolnych ( w skryptowych wzorach wektor R
    {

        for (int i = 0; i <= m; i++)
            B[i + 1] = LiczR(i);
}

void WartosciFunkcjiAproksymujacej() // funkcja obliczajacawartosc funkcji aproksymowanej w punktach
{  for (int i = 0; i <e; i++)
{   yy[i] = 0;
    for (int j = 1; j <= m+1; j++)
    yy[i] += X[j] * pow(x[i], j-1);
}
}

void ObliczBlad() // funkcja obliczajcablad
{
    double tmp = 0;
    for (int i = 0; i <e; ++i)
        tmp += pow(yy[i] - y[i], 2);
    blad = sqrt(tmp / e);
}

void SaveToFile()
{
    string filename;
    cout << "Podaj nazwe pliku gdzie zapisac raport: ";
    cin >> filename;
    ofstream output (filename.c_str(),ios_base::out);
    output << "Liczba wezlow: " << e;
    output << "\nStopien wielomianu: " << m;
    output << "\nPrzedzial: [" << a << ", " << b << "]";
    output.setf(ios_base::scientific);
    output.precision(12);
    output << "\nBlad aproksymacji: " << blad;
    output.setf(ios_base::fixed, ios_base::floatfield);
    output << "\nTabela wartości funkcji stablicowanej i aproksymujacej (dokladnosc do 12 miejsc po przecinku):\n";
    output << "Wezel Wartoscdokladna Wartosc funkcji aproksymujacej\n";
    for (long i = 0; i <e; ++i)
    {
        output << " ";
        output << x[i];
        output << " ";
        output << y[i];
        output << " ";
        output << yy[i] << endl;
    }
    for (long i = 0; i < e; ++i)
    {
        output << yy[i] << endl;
    }
    output.unsetf(ios_base::fixed);
    output.setf(ios_base::scientific);
    output << "\nWspolczynniki A\n";
    for (int i = 0; i <= m; ++i)
    {
        output.width(3);
        output << i << ": ";
        output.width(21);
        output << X[i] << endl;
    }
    output.close();
}

int main() {

    GenerujTablice();
    GenerujMacierz();
    GenerujB();
    macierzeUL();
    rozwiaz();
    WartosciFunkcjiAproksymujacej();
    ObliczBlad();
    SaveToFile();
    getchar();
}