taylor cos

//    #############################################################
//    Este programa calcula o erro percentual da Série de Taylor
//    da função cos(bx), expandida no ponto x0
#define DIM1 50
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <iostream>
    main()
    {
      double a[DIM1];
      double x0,x,b,fxa,fxt,erro_p,Deltax,xa,xb,h,Pi;
      float fat(int n);
      FILE *ent;
      FILE *sai;
      int i,j,k,n,icont,nx,nt,aux,grau;
      ent = fopen("taylor150510v2.1","r");
      sai = fopen("taylor150510v2.s.2","w");
      fscanf(ent,"%lf",&x0);           /* Ponto de expansão da série */
      fscanf(ent,"%i",&nt);            /* grau + 1  */
      fscanf(ent,"%lf",&b);            /* Coeficiente multiplicador b */
      fscanf(ent,"%i",&nx);            /* Numero de abscissas  */
      fscanf(ent,"%lf   %lf",&xa,&xb); /* Dominio xa <= x <= xb */
      grau= nt-1;
      Deltax=xb-xa;
      h=Deltax/(nx-1);
      for ( n = 0 ; n <= grau ; n++)
      {
        i=1;
        if(n%2==0)
        {                                // Par
          j=0;
          k=1;
          if(n%4==0) i=2;                // múltiplo de 4
        }
        else
        {                               //  Ímpar
          j=1;
          k=0;
          if((n-1)%4==0) i=2;             // múltiplo de 4 + 1
        }
        a[n]=pow(-1.,i)*pow(b,n)*pow(sin(b*x0),j)*pow(cos(b*x0),k);
      }
      x=xa;
      /*      Cabeçalho     */
      fprintf(sai,"      x            f(x)             F(x)          Ep(%)\n");
      for (icont=1;icont<=nx;icont++)
      {
        fxt=0.;
        for (n=0;n<=grau;n++)
        {
          fxt=fxt+a[n]*pow(x-x0,n)/fat(n);  /* Solução numérica F(x) */
        }
        fxa=cos(b*x);                       /* Solução analítica  f(x) */

        if(fxa != 0.)
        {
          erro_p=fabs((fxt-fxa)/fxa)*100.;  /*  Erro percentual  */
          fprintf(sai,"%12.8f %16.8f %16.8f %16.8e  \n",x,fxa,fxt,erro_p);
        }
        else
        {
          fprintf(sai,"%12.8f %16.8f %16.8f     -  \n",x,fxa,fxt);
        }
        x = x +h;
      }
      fclose(ent);                         /*  Fechando os arquivos   */
      fclose(sai);
      system ("pause");
    }
// ***********************************************************************
   double fat(int n)
   {
     int i;
     double xi,fat;
     fat=n;
     for (i=n; i>=2; i--)
     {
       xi = i ;
       fat=fat*(xi-1.);
     }
     if(n==0) fat=1.;
     return(fat);
   }