#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <conio.h>

int n,m;
double x[40], y[40], a[40], b[40], r[40], mat[40][40], sol[40], sol2[40][5];

void sistema()
{
    int i,j;
    for(i=0; i<n; i++)
    {
        for(j=0; j<n; j++)
        {
            mat[i][j] = pow(x[i],j);
        }
        r[i] = y[i];
    }
}

void resolversistema()
{
    int i,k,j;
    double s;
    for(i=0;i<n;i++)
    {
    sol[i]=x[i];
    }
    sol[n-1]=r[n-1]/mat[n-1][n-1];
    for(k=n-2; k>=0; k--)
    {
        s=0;
        for(j=k+1; j<n; j++)
            s+=mat[k][j]*sol[j];
        sol[k]=(r[k]-s)/mat[k][k];
    }
}

void escalonamento()
{
    int i,j,k;
    double mult;
    for(k=0; k<(n-1); k++)
    {
        for(i=k+1; i<n; i++)
        {
            mult=mat[i][k]/mat[k][k];//multiplicador
            for(j=k+1; j<n; j++)
            {
                mat[i][j] = mat[i][j]-mult*mat[k][j];
                //imprimeMatriz();
            }
            r[i]=r[i]-mult*r[k];
            mat[i][k]=0;
        }
    }
}

double fx(double x)
{
        int i;
        double f=0;
        for(i=0;i<n;i++)
        {
                f = f+ (sol[i]*pow(x,i));
        }
        return f;
}

double integralsimpson(double a, double b, int partes)
{
        double h, s, f, rr, fa, fb, x;
        int i;
        h=(b-a)/partes;
        x=a;
        s=0;
        for(i=0;i<(partes-1)/2;i++)
        {
                x=x+h;
                f=fx(x);
                s=s +(4*f);
                x=x+h;
                f=fx(x);
                s=s+ (2*f);
        }
        fa=fx(a);
        fb=fx(b);
        s = s + fa + fb;
        rr=(h/3)*s;
        printf("\n%lf", rr);
        return rr;
}

double integraltrapezio(double a, double b, double partes)
{
        double h, s, f, rr, fa, fb, x;
        int i;
        h=(b-a)/partes;
        x=a;
        s=0;
        for(i=0;i<partes-1;i++)
        {
                x=x+h;
                f=fx(x);
                s=s+2*f;
        }
        fa=fx(a);
        fb=fx(b);
        s=s+ fa +fb;
        rr=(h/2)*s;
        printf("\n%lf\n", rr);
        return rr;
}
main()
{
        int i, j,V;
        double aux;
        FILE *dados,*saida;
        dados=fopen("entrada.txt", "r");
        saida=fopen("resultado.txt", "w");
        if(dados==NULL)
        {
                printf("\n Erro ao abrir arquivo");
        }
        else
        {
                fscanf(dados, "%d", &n);
                for(i=0;i<n;i++)
                {
                        fscanf(dados, "%lf", &x[i]);
                        fscanf(dados, "%lf", &y[i]);
                }
                fscanf(dados, "%d", &m);
                printf("valor de m: %d", m);
                for(i=0;i<m;i++)
                {
                        fscanf(dados, "%lf", &a[i]);
                        fscanf(dados, "%lf", &b[i]);
                }
            for(i=0;i<m;i++)
              printf("%lf ", b[i]);
                sistema();
                escalonamento();
                resolversistema();
                fprintf(saida,"POLINOMIO INTERPOLADOR\np(x)= ");
                for(i=n-1;i>=0;i--)
                {
                   fprintf(saida,"%lf", sol[i]);

                   if(i)
                     fprintf(saida, " X");
                   if(i>1)
                      fprintf(saida, "^%d", i);
                   if(i<=n-1 && i!=0)
                      fprintf(saida, " + ");
            }

            fprintf(saida,"\n\nAPROXIMACOES\n");

            for(i=0;i<m;i++)
            {
                sol2[i][0]=integraltrapezio(a[i],b[i], 10);
            }

            for(i=0;i<m;i++)
            {
                sol2[i][1]=integraltrapezio(a[i],b[i], 100);
            }

            for(i=0;i<m;i++)
            {
                sol2[i][2]=integraltrapezio(a[i],b[i], 1000);
            }
            for(i=0;i<m;i++)
            {
                sol2[i][3]=integralsimpson(a[i],b[i], 10);
            }

            for(i=0;i<m;i++)
            {
               sol2[i][4]=integralsimpson(a[i],b[i], 100);
            }

            for(i=0;i<m;i++)
            {
                 sol2[i][5]=integralsimpson(a[i],b[i], 1000);
            }

            printf("%lf", sol2[0][0]);

            for(i=0;i<m;i++)
            {
                for(j=0;j<6;j++)
                {
                        printf("%lf ",&sol2[i][j]);
                }
                printf("\n");
            }
            printf("%lf", sol2[0][0]);
            for(i=0;i<m;i++)
            {
                V=10;
                for(j=0;j<6;j++)
                {
                        if(j==0 || j==1 || j==2)
                        {

                         fprintf(saida,"{T(%d,%d)}=%lf  ",i+1,V,sol2[i][j]);
                         V=V*10;
                        }

                        else
                        {
                            V=pow(10,(j-2));
                            fprintf(saida,"{S(%d,%d)}=%lf  ",i+1,V,sol2[i][j]);

                        }

                }
                fprintf(saida,"\n");
            }
            fprintf(saida,"%lf", sol2[0][0]);
            printf("\nTudo calculado. Pode fechar, amigo.");
    }
}