#include <math.h>
#include <stdio.h>
#define N 20
double f(double x){
    return 3*sin(x);
}

int main(int argc, char **argv)
{
    //Зчитуємо результат табуляції функції
    FILE * fp;
    double arr[N][2];
    fp = fopen("test.txt", "r");
    for(int j = 0;j<N;j++){
        fscanf(fp,"%lf %lf\n",&arr[j][0],&arr[j][1]);
        printf("%lf %lf\n",arr[j][0],arr[j][1]);
    }
    double a[N], b[N], c[N], d[N], h[N];
    double alph[N],beta[N],gama[N],delta[N];
    FILE  *newx, *fe, *newy;


    for(int i = 1; i < N; i++){
        h[i] = arr[i][0] - arr[i-1][0];
        a[i] = arr[i-1][1];
    }
    //Розрахунок коефіцієгтів альфа бета гама дельта
    for(int i = 2; i < N - 1 ; i++){
        alph[i] = h[i-1];
        beta[i-1] = 2*(h[i-1] + h[i]);
        gama[i] = h[i];
        delta[i-1] = 3 * (((arr[i][1] - arr[i-1][1]) / h[i]) - ((arr[i-1][1] - arr[i-2][1]) / h[i-1]));
    }

    alph[N-2] = h[N-2];
    beta[N-2] = 2*(h[N-2] + h[N-1]);
    delta[N-2] = 3 * (((arr[N-1][1] - arr[N-2][1]) / h[N-1]) - ((arr[N-2][1] - arr[N-3][1]) / h[N-2]));

    //Прямий хід методу прогонки
    double A[N], B[N];
    A[1] = -(gama[1] / beta[1]);
    B[1] = delta[1] / beta[1];
    for(int i = 2; i < N-1; i++){
        A[i] = -(gama[i] / (alph[i] * A[i-1] + beta[i]));
        B[i] = (delta[i] - alph[i]*B[i-1]) / (alph[i] * A[i-1] + beta[i]);
    }
    //Зворотній хід методу прогонки

    c[N-1] = (delta[N-1] - alph[N-1]*B[N-2]) / (alph[N-1]*A[N-2] + beta[N-1]);
    for(int i = N-1; i > 1; i--){
        c[i] = A[i-1] * c[i+1] + B[i-1];
    }
    c[1]=0;

    //Розрахунок коефіцієнтів b і d сплайнів
    for(int i = 1; i < N-1; i++){
        d[i] = (c[i+1]-c[i])/(3*h[i]);
        b[i] = (arr[i][1] - arr[i-1][1]) / h[i] - h[i] * (c[i+1] + 2*c[i])/3;}
    d[N-1] = -c[N-1]/(3*h[N-1]);
    b[N-1] = (arr[N-1][1] - arr[N-2][1]) / h[N-1] - 2/3*h[N-1]*c[N-1];

    newx=fopen("newx.txt","w");
    fe=fopen("e.txt","w");
    newy=fopen("newy.txt","w");
    double yy[20*N];
    double xx[20*N];
    double H = (arr[N-1][0] - arr[0][0])/(20*N);
    xx[0] = arr[0][0];
    int L=1;
    double z = arr[0][0];

    for(int i = 0; i < 20*N; i++){
        while(z>arr[L][0]) L++;
        yy[i] = a[L] + b[L] * (z - arr[L-1][0]) + 	c[L] * (z - arr[L-1][0]) * (z - arr[L-1][0]) + d[L] * (z - arr[L-1][0]) * (z - arr[L-1][0]) * (z - arr[L-1][0]);
        fprintf(newx, "%lf\n", z);
        fprintf(fe, "%lf\n",fabs(yy[i]-f(z)));
        fprintf(newy, "%lf\n",yy[i]);
        xx[i] = z;
        z += H;
        if(z > arr[N-1][0]) break;
    }

    fclose(newx);
    fclose(fe);
    fclose(newy);
    return 0;
}