/*
*   NEWTON'S INTERPOLATORY DIVIDED-DIFFERENCE FORMULA ALGORITHM 3.2
*
*   To obtain the divided-difference coefficients of the interpolatory
*   polynomial P on the (n+1) distinct numbers x(0), x(1), ..., x(n)
*   for the function f:
*
*   INPUT:   numbers x(0), x(1), ..., x(n); values f(x(0)), f(x(1)), ...,
*            f(x(n)) as the first column Q(0,0), Q(1,0), ..., Q(N,0) OF Q,
*            or may be computed if function f is supplied.
*
*   OUTPUT:  the numbers Q(0,0), Q(1,1), ..., Q(N,N) where
*            P(x) = Q(0,0) + Q(1,1)*(x - x(0)) + Q(2,2)*(x - x(0))*(x - x(1))
*            +... + Q(N,N)*(x - x(0))*(x - x(1))*...*(x - x(N - 1)).
*/

#include<stdio.h>
#include<math.h>
#include <iostream>
#include <fstream>
#define true 1
#define false 0
using namespace std;

double F(double);
void INPUT(int *, double *, double [][26], int *);
void OUTPUT(ofstream *);


main()
{
   double Q[26][26],X[26];
   int I,J,N,OK;
   ofstream myfile;

   INPUT(&OK, X, Q, &N);
   if (OK) {
      OUTPUT(&myfile);
      /* STEP 1 */
      for (I=1; I<=N; I++)
         for (J=1; J<=I; J++)
            Q[I][J] = (Q[I][J-1] - Q[I-1][J-1]) / (X[I] - X[I-J]);
      /* STEP 2 */
      if (myfile.is_open()) 
      {
          myfile << "Input data follows:\n";
          for (I=0; I<=N; I++)
               myfile << "X("<<I <<") = "<< X[I] <<"  F(X("<<I<<")) = "<< Q[I][0] <<"\n"; 
          myfile << "\nThe coefficients Q(0,0), ..., Q(N,N) are:\n";
          
          for (I=0; I<=N; I++) myfile << Q[I][I]<<"\n"; 
          
          myfile.close();
      }
      else
      {
          printf("Input data follows:\n");
          for (I=0; I<=N; I++)
              printf("X(%d) = %12.8f F(X(%d)) = %12.8f\n", I, X[I], I, Q[I][0]);
          printf("\nThe coefficients Q(0,0), ..., Q(N,N) are:\n");
          for (I=0; I<=N; I++) printf("%12.8f\n", Q[I][I]);
      }
   }
   return 0;
}

/* Change F if program is to calculate the first column of Q */
double F(double X)
{
   double f; 

   f =  1.0/X;
   return f;
}

void INPUT(int *OK, double *X, double Q[][26], int *N)
{
        int I, FLAG;
        char A;
        char NAME[30];
        FILE *INP; 

         printf("Newtons form of the interpolation polynomial\n");
         printf("Warning: total number of nodes should be less than 26\n");
         *OK = true;

         printf("Has a text file been created with the data in two columns ?\n");
         printf("Enter Y or N\n");
         scanf("\n%c", &A);
         if ((A == 'Y') || (A == 'y')) {
            printf("Input the file name in the form - ");
            printf("drive:name.ext\n");
            printf("For example:   DATA.DTA\n");
            scanf("%s", NAME);
            INP = fopen(NAME, "r");
            *OK = false;
            while (!(*OK)) {
               printf("Input n, which is number of points - 1\n");
               scanf("%d", N);
               if (*N > 0) {
                  for (I=0; I<=*N; I++) 
                     fscanf(INP, "%lf %lf", &X[I], &Q[I][0]);
                  fclose(INP);
                  *OK = true;
               }
               else printf("Number must be a positive integer\n");
            }
         }
         else {
            printf("Please create the input file in two column ");
            printf("form with the X values and\n");
            printf("F(X) values in the corresponding columns.\n");
            printf("The program will end so the input file can ");
            printf("be created.\n");
            *OK = false;
         }
}


void OUTPUT(ofstream *OUP)
{
   int FLAG;
   char NAME[30];

   printf("Select output destination\n");
   printf("1. Screen\n");
   printf("2. Text file\n");
   printf("Enter 1 or 2\n");
   scanf("%d", &FLAG);
   if (FLAG == 2) {
      printf("Input the file name in the form - drive:name.ext\n");
      printf("For example:   OUTPUT.DTA\n");
      scanf("%s", NAME);
      OUP->open(NAME);
   }
   cout<<"NEWTONS INTERPOLATION POLYNOMIAL\n\n";
}