Programmer 7500/MAX_PROGRAMMERS.iso/PASCAL/TPTOOL4.ZIP/POLY.INC

Text File | 1987-03-28 | 5.5 KB | 242 lines

#log Polynomial interpolation

(*
 *  This program fits least-squares polynomials to bivariate
 *  data, using an orthogonal polynomial method.
 *
 *  Equation form:
 *       Y := C0 + C1 * X + C2 * X^2 + C3 * X^3 ...  + Cn * X^n
 *
 *  The most common use of this technique is for computerizing
 *  graphical or tabular data.  the equation requires much less
 *  space and no interpolating.
 *
 *  see louis g. kelly -- handbook of numerical methods and
 *                        applications, p. 68
 *
 *
 * usage:
 *    include poly.inc
 *    declare X,Y as data_array
 *    declare C as const_array
 *    put datapoints into X,Y
 *    call solve_polynomial(degree,X,Y,count,C)
 *    returns array of C constants
 *    degree = 0 finds C[0]
 *    degree = N finds C[0]..C[n]
 *
 *)


const
   max_degree = 30;   {highest constant subscript}
   max_data = 50;    {maximum number of datapoints}

type
   const_array = array[0..max_degree] of real;  {array of constants}
   data_array  = array[0..max_data] of real;    {array of datapoints}


(* use current interpolation data to find Y for a given X *)
function poly_lookup(x:      real;
                     var C:  const_array;
                     degree: integer): real;
var
   i:  integer;
   y0: real;
   x0: real;

begin
   x0 := 1;
   y0 := 0;
   for i := 0 to degree do
   begin
      y0 := y0 + x0 * C[i];

      if abs(x0) < 1e17 then   {prevent overflow}
         x0 := x0 * x;
   end;

   poly_lookup := y0;
end;



(* solve the polynomial equation for a given table of data
   input in X,Y,count;  passes back array of constants in C;
   calculates constants up to C[min_degree]. *)

procedure solve_polynomial(min_degree:  integer;       {input}
                           var X,Y:     data_array;    {input}
                           count:       integer;       {output}
                           var C:       const_array);  {output}
var
   degree:  integer;
   points,u,v,x0: real;
   A,D: const_array;
   P: array[0..max_degree] of const_array;


   (* subroutine for the computation of P(degree,J), J := 0,...,degree,      *)
   (* from U[degree],V[degree-1],&P[I] := (X - U[I])*P[I-1] - V[I-1]*P[I-2]  *)
   procedure calc1;
   var
      i7,i8: integer;
      j:     integer;

   begin
      i8 := degree - 1;
      i7 := 0;

      if (degree <> 1) then
         i7 := i8 - 1;

      P[degree][0] := - u * P[i8][0] - v * P[i7][0];

      if (i8 <> 0) then
         for j := 1 to i8 do
            P[degree][j] := P[i8][j-1] - u * P[i8][j] - v * P[i7][j];

      P[degree][degree] := P[i8][i8];
   end;


   (* subroutine for the evaluation of P0 := P[I0](X0) -- P[I] is *)
   (* the i-th polynomial.                                       *)
   function calc2(i0: integer; x0: real): real;
   var
      j:  integer;
      p0: real;

   begin
      p0 := 1;

      if (i0 <> 0) then
      begin
         p0 := P[i0][i0];
         for j := i0-1 downto 0 do
            p0 := p0 * x0 + P[i0][j];
      end;

      if p0 > 1e18  then p0 := 1e18;
      if p0 < -1e18 then p0 := -1e18;
      calc2 := p0;
   end;


   (* evaluation of p(i,j), the j-th coefficient of the i-th   *)
   (* orthogonal polynomial. *)
   procedure find_poly_term;
   var
      j,k: integer;
      p0:  real;
      p1:  real;

   begin
      (* computation of U[I] *)
      u := 0;

      for j := 1 to count do
      begin
         x0 := X[j];
         p0 := calc2(degree-1,x0);

         if abs(p0) < 1e17 then       {prevent overflow}
            u := u + x0 * p0 * p0;
      end;

      u := u / D[degree-1];
      v := 0;

      if (degree <> 1) then
      begin

         (* computation of V[I-1] *)
         for j := 1 to count do
         begin
            x0 := X[j];
            p1 := calc2(degree-1,x0);

            (* computation and storage of P(I,J), J := 0,...,I *)
            p0 := calc2(degree-2,x0);

            if abs(p0) < 1e17 then       {prevent overflow}
               v := v + x0 * p1 * p0;
         end;

         v := v / D[degree-2];
      end;

      calc1;

      (* computation of A(I), the coefficients of the orthogonal *)
      (* basis which determine the fit.  *)
      A[degree] := 0;
      D[degree] := 0;

      for j := 1 to count do
      begin
         x0 := X[j];
         p0 := calc2(degree,x0);
         A[degree] := A[degree] + Y[j] * p0;
         D[degree] := D[degree] + p0 * p0;
      end;

      A[degree] := A[degree] / D[degree];

      (* evaluation of the coefficients C(I) of the fit polynomial. *)
      for k := 0 to degree do
         C[k] := C[k] + A[degree] * P[degree][k];
   end;


   procedure solve_linear_fit;
   var
      i,j: integer;
      mean: real;

   begin
      for i := 0 to max_degree do
      begin
         A[i] := 0;
         C[i] := 0;
         D[i] := 0;
         for j := 0 to max_degree do
            P[i][j] := 0;
      end;

      mean := 0;
      for i := 1 to count do
         mean := mean + Y[i];

      mean := mean / points;

      P[0][0] := 1;
      A[0] := mean;
      C[0] := mean;
      D[0] := points;
   end;



(* body of solve_polynomial *)

begin

   if (min_degree >= count) then
      min_degree := count-1;

   if (min_degree > max_degree-1) then
      min_degree := max_degree-1;

   points := int(count);

   solve_linear_fit;

   for degree := 1 to min_degree do
      find_poly_term;

end;