newton_divided_difference.c

#include<stdio.h>
#include<stdlib.h>
#define N 1010
#define INF 1e5

typedef struct pt
{
  double x , y;
} pt;

pt a[N];
double f[N][N];
double b[N];
double A[N];

void gen(int n , double roots[N] , double k)
{
    double coeff[N];
    int i , j;
    for(i = 0 ; i<=n ; ++i)
        coeff[i] = 0;
    coeff[n] = 1;
    for(i = 1; i<=n ; ++i)
    {
        for(j = n - i - 1 ; j<n ; ++j)
        {
            coeff[j] = coeff[j] + (-1)*roots[i-1]*coeff[j+1];
        }
    }

    for(i = n ; i>=0 ; --i)
    {
        A[i] += (coeff[i]*k);
    }
}

double ff(int l , int r)
{
  double res;
  if(f[l][r] != INF)
    return f[l][r];

  if(l == r)
  {
    f[l][l] = a[l].y;
    return f[l][l];

  }

  res = (ff(l+1 , r) - ff(l,r-1)) / (a[r].x - a[l].x);
  f[l][r] = res;
  return res;
}
void calculate_b(int n)
{
  int i,k;
  for(k = 0 ; k<=n ; ++k)
  {
    b[k] = ff(0,k);
  }
}

double P(double x , int n)
{
  double res = 0 , tmp;
  int i, j;
  for(i = 0 ; i<=n ; ++i)
  {
    tmp = b[i];
    for(j = 0 ; j<i ; ++j)
    {
      tmp *= (x - a[j].x);
    }
    res+=tmp;
  }
  return res;
}

void print_polynomial(int n)
{
  int i , j;
  printf("P(x) = ");
  for(i = 0 ; i<=n ; ++i)
  {
    printf(" %.2lf " , b[i]);
    
    for(j = 0 ; j<i ; ++j)
      printf("(x - %.2lf)" , a[j].x);
    if(i != n)
      printf(" + ");
  }
  printf("\n\n");
}

void print_simplified_polynomial(int n)
{
  printf("Simplified: \n");
  int i , j;
  int count = 0;
  double roots[N];
  for(i = 0 ; i<=n ; ++i)
  {
    count = 0;
    for(j = 0 ; j<N ; ++j)
    {
      if(j < i)
      {
        roots[j] = a[j].x;
        count++;
      }
      else
      {
        roots[j] = 0;
      }
    }
    gen(count , roots , b[i]);
  }

  printf("\nP(X) = ");
  for(i = n ; i>=0 ; --i)
  {
    if(i == 0)
      printf("(%.2lf)" , A[i]);
    else if(i == 1)
      printf("(%.2lf)*(X) + " , A[i]);
    else
      printf("(%.2lf)*(X^%d) + " , A[i] , i);
  }
  printf("\n\n");

}
signed main()
{

  int n , i , j;
  double x;
  for(i = 0 ; i<N ; ++i)
    for(j = 0 ; j<N ; ++j)
      f[i][j] = INF;

  printf("\t--NEWTON DIVIDED DIFFERENCE--\n\n");
  printf("Enter number of data points: ");
  scanf("%d" , &n);
  n-=1;
  printf("\nEnter the data points:\n");
  for(i = 0 ; i<=n ; ++i)
  {
    scanf("%lf %lf" , &a[i].x , &a[i].y);
  }

  calculate_b(n);
  print_polynomial(n);
  print_simplified_polynomial(n);
  
  printf("Enter Value of x: ");
  scanf("%lf" , &x);
  printf("\n\nP(%.0lf) = %.4lf\n\n" , x , P(x , n));






  return 0;
}