(*
* NEVILLE'S ITERATED INTERPOLATION ALGORITHM 3.1
*
* To evaluate the interpolating polynomial P on the
* (n+1) distinct numbers x(0), ..., x(n) at the number x
* for the function f:
*
* INPUT:   numbers x(0), ..., x(n) as XX(0), ..., XX(n);
*          number x; values of f as the first column of Q
*          or may be computed if the function f is supplied.
*
* OUTPUT:  the table Q with P(x) = Q(n+1,n+1).
*)
OK = 0;
While[OK == 0,
   FLAG=Input["This is Neville's Method.\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]<>")"];
            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",
      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 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,
   XX=Input["Input the point at which the polynomial\n 
   is to be evaluated\n"];
];
If[OK == 1,
(*STEP 1*)
   Y[0]=Evaluate[XX-X[0]];
   For[i=1,
      i<=n,
      i++,
      Y[i]=Evaluate[XX-X[i]];
      For[j=1,
         j<=i,
         j++,
         Q[i,j]=Evaluate[(Y[i]*Q[i-1,j-1]
                -Y[i-j]*Q[i,j-1])/(Y[i]-Y[i-j])];
      ];
   ];
(*STEP 2*)
   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,"Neville's Method\n"];
   Write[OUP,"Table for P evaluated at X = ",SetPrecision[XX,10],"  follows:\n"];
   Write[OUP,"Entries are XX(i), Q(i,0), ..., Q(i,i) for each i = 0, ..., n \n"]; 
   Write[OUP,"   where n = \n",n];
   For[i=0,
      i<=n,
      i++,
      Write[OUP,"X(",i,") = ",SetPrecision[X[i],10]];
      For[j=0,
         j<=i,
         j++,
         Write[OUP,"Q(",i,",",j,") = ",SetPrecision[Q[i,j],9]];
      ];
      Write[OUP,"\n"];
   ];
   If[OUP == "OutputStream[",NAME," 3]",
      Print["Output file: ",NAME," created successfully\n"];
      Close[OUP];
   ];
];
Quit[];