#include
#include
#include
#include

#define MAXN 101

int n; // maksymalny indeks węzła interpolacji,

int m; // stopien wielomianu

double da, db; // przedział

double h;

double x[MAXN]; //tablica węzłów

double f[MAXN]; //tablica warto?ci funkcji w węzłach

char *nazwa; //nazwa pliku


double sqr(double t)
{
return (t*t);
}

double funkcja(double t)
{
return (sin(t)); //funkcja podana przez prowadzšcego
}

// Wyznaczenie węzłów i warto?ci funkcji w węzłach
void generWezly()
{
int i;
h=(db-da)/n;
for (i=0; i<=n; i++)
{
x[i]=da+i*h;
f[i]=funkcja(x[i]);
printf("%f %f \n", x[i], f[i]);
}
}

double silnia(int t)
{
double wynik=1.0;
int i;
for (i=2; i<=t; i++)
wynik = wynik * i;
return wynik;
}

double symbNewtona(int k, int s)
{
return (silnia(k) / (silnia(s) * silnia (k-s)));
}


double iloczWiel(double q, int s) // iloczyn q[s]/n[s]
{
int i;
double wielQ=1.0;
double wielN=1.0;
for (i=1; i<=s; i++)
{
wielQ = wielQ * (q - i + 1);
wielN = wielN * (n - i + 1);
}
return (wielQ/wielN);
}

double wielGrama(int k, double q)
{
int s;
double wynik=0.0;

if (k==0)
return 1.0;

for (s=0; s<=k; s++)
{
if (s % 2)
wynik = wynik - symbNewtona(k, s) * symbNewtona(k+s, s) * iloczWiel(q, s);
else
wynik = wynik + symbNewtona(k, s) * symbNewtona(k+s, s) * iloczWiel(q, s);
}
return wynik;
}


double s_j(int j)
{
int q;
double wynik=0.0;
for (q=0; q<=n; q++)
wynik = wynik + sqr(wielGrama(j, q));
return wynik;
}

double c_j(int j)
{
int q;
double wynik=0.0;
for (q=0; q<=n; q++)
wynik = wynik + f[q] * wielGrama(j, q);
return wynik;
}

double funAprox(double t)
{
double wynik=0.0;
int j;
for (j=0; j<=m; j++)
wynik = wynik + (c_j(j) / s_j(j)) * wielGrama(j, (t - x[0]) / h);
return wynik;
}

double generPlik()
{
FILE *plk;
double blad=0.0;
plk = fopen( nazwa, "w" );
fprintf(plk, "xi\t f(xi)\t F(xi)\n");
for (int i=0; i<=n; i++)
{
fprintf(plk, "%f\t %f\t %f\n", x[i], f[i], funAprox(x[i]));
blad = blad + sqr(f[i] - funAprox(x[i]));
}
fprintf(plk, "*InfKon*\nV 0 0 0 %d",n+1 );
fclose(plk);
return sqrt(blad);
}

int main(int argc, char* argv[])
{
// Odczyt danych z linii polecenia

da=atof(argv[1]);
db=atof(argv[2]);
n=atoi(argv[3]);
m=atoi(argv[4]);
nazwa=argv[5];

// Wyprowadzenie na ekran danych z linii polecenia
printf("a= %f\n", da);
printf("b= %f\n", db);
printf("n= %d\n", n);
printf("m= %d\n", m);

generWezly();
printf("\nBlad sredniokwadratowy punktowy wynosi: %f", generPlik());
return 0;
}