(*
*NEWTON'S INTERPOLATORY DIVIDED-DIFFERENCE FORMULA ALGORITHM 3.2
*
* To obtain the divided-difference coeffieients 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,0), ..., 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)).
*)
OK = 0;
While[OK == 0,
   FLAG=Input["Newton's form of the interpolation polynomial\n
   \n
   Choice of input method\n
   1. Input entry by entry from keyboard\n
   2. Input data from a text file\n
   3. Generate data using a function\n"];
   If[FLAG == 1 || FLAG == 2 || FLAG == 3,
      OK = 1;
   ]; 
];
If[FLAG == 1,
   OK = 0;
   While[OK != 1,
      n = Input["Input n"];
      If[n > 0,
         OK = 1;
         For[i = 0,
            i <= n,
            i++,
            X[i] = Input["Input X("<>ToString[i]<>")"];
            X[i] = N[X[i]];
            Q[i,0] = Input["Input F(X("<>ToString[i]<>"))"];
         ],
         Input["Number must be a positive integer.\n
         Enter 1 [enter] to continue\n"]
      ];
   ];
];
If[FLAG == 2,
   A = InputString["Has a text file been created with the data\n
            in two columns?\n\n
            Enter 'Yes' or 'No'\n"];
   If[A == "y" || A == "Y" || A == "Yes" || A == "yes",
      NAME = InputString["Input the file name\n
               In the form - name.ext\n
               For example:   output.dta\n"]; 
      INP = OpenRead[NAME];
      OK = 0;
      While[OK == 0,
         n = Input["Input n\n"];
         If[n > 0,    
            For[i = 0,
               i <= n,
               i++,
               X[i] = Read[INP,Number];
               Q[i,0] = Read[INP,Number];
            ];
            Close[INP];
            OK = 1,
            Input["Number must be a positive integer\n
            \n
            Press 1 [enter] to continue\n"];
         ];
      ],
      Input["Please create the input file in two column form\n
      with the X values and the F(X) values in the\n
      corresponding columns.\n
   The program will end so the input file\n
   can be created.\n
      Press 1 [enter] to continue\n"];
      OK = 0;
   ];
];  
If[FLAG == 3,
   TEMP = Input["Input the function F(X) in terms of x.\n
   \n
   For example: Cos[x]\n"];
   F[x_] := Evaluate[TEMP];
   OK = 0;
   While[OK == 0,
   n = Input["Input n"];
      If[n > 0,
         OK = 1;
         For[i = 0,
            i <= n,
            i++,
            X[i] = Input["Input X("<>ToString[i]<>")"];
            X[i]=N[X[i]];  
            Q[i,0] = N[F[X[i]]];
         ];
         OK = 1,
         Input["Number must be a positive integer.\n
         Enter 1 [enter] to continue\n"];
      ];
   ];
];
If[OK == 1,
   FLAG = Input["Select output destination\n
                 1. Screen\n
                 2. Text file\n
                 Enter 1 or 2\n"];
   If[FLAG == 2,
      NAME = InputString["Input the file name\n
                    For example:   output.dta\n"];
      OUP = OpenWrite[NAME,FormatType->OutputForm],
      OUP = "stdout";
      Write[OUP,"NEWTON'S INTERPOLATION POLYNOMIAL\n"];
   ];
(*STEP 1*)
   For[i = 1,
      i <= n,
      i++,
      For[j = 1,
         j <= i,
         j++,
         Q[i,j] = N[(Q[i,j-1]-Q[i-1,j-1])
                       /(X[i]-X[i-j])];
      ];
   ];
(*STEP 2*)
   Write[OUP,"Data input follows:\n"];
   For[i = 0,
      i <= n,
      i++,
      Write[OUP,"X(",i,") = ",SetPrecision[X[i],10],
    "    F(X(",i,")) = \n",SetPrecision[Q[i,0],9]];
   ];
   Write[OUP,"The coefficients Q(0,0), ..., Q(n,n) are:\n"];
   For[i = 0,
      i <= n,
      i++,
      Write[OUP,"\n",SetPrecision[Q[i,i],9]];
   ];
   If[OUP == "OutputStream[",NAME," 3]",
      Print["Output file: ",NAME," created successfully\n"];
      Close[OUP];
   ];
];
Quit[];