(* :Name: NumericalMath`PolynomialFit` *)

(* :Title: Curve Fitting over Polynomial Vector Spaces in a Numerically
	Stable Way *)

(* :Author: Jerry B. Keiper *)

(* :Summary:
This package does polynomial least squares fitting of data and avoids
the numerical instability that often results if the basis functions are
chosen to be simple powers of a variable.  Shifted Chebyshev
polynomials of the first kind are used as a basis for the vector space
of polynomials in a single variable.  The shift is chosen based on the
distribution of the abscissas of the data.  The result is returned as a
pure function, partly for convenience and partly to discourage the
user from trying to rewrite it as a linear combination of simple
powers.  Such "simplification" would destroy much of the numerical
stability gained by this method.
*)

(* :Context: NumericalMath`PolynomialFit` *)

(* :Package Version: 2.0 *)

(* :Copyright: Copyright 1990-1999, Wolfram Research, Inc.
*)

(* :History:
	Originally written by Jerry B. Keiper, March 1989.
	Revised by Jerry B. Keiper, December 1990.
*)

(* :Keywords: Fit, least squared error, polynomial approximation *)

(* :Source:
	G. Dahlquist & A. Bjorck, Numerical Methods, Prentice-Hall,
		Englewood Cliffs, NJ, 1974.
*)

(* :Warnings: None. *)

(* :Mathematica Version: 2.0 *)

(* :Limitations:
	Although much more stable than Fit, PolynomialFit is still
	affected by floating-point round-off.  It allows the user to evaluate
	the result at a symbolic argument (say x) and expand to get a result
	expressed in terms of basis {1, x, x^2, ..., x^n}.  Such abuse
	is strongly discouraged due to the numerical stability that
	results.
*)


BeginPackage["NumericalMath`PolynomialFit`"]

FittingPolynomial::usage =
"FittingPolynomial[<>, n] represents a polynomial of degree n.  It is returned
by PolynomialFit."

PolynomialFit::usage =
"PolynomialFit[data, n] gives the least squares polynomial fit to data
much as Fit[data,{1,x,x^2,...,x^n}, x] would do, except that the calculation
is done in a way that is numerically stable and the result is given as a
FittingPolynomial that is numerically stable.  Changing the form of the result
is NOT advised for reasons of numerical stability."

Begin["NumericalMath`PolynomialFit`Private`"]

PolynomialFit[data_List, n_Integer] :=
    Module[{len = Length[data], answer},
	FittingPolynomial[answer, n] /;
		(answer = ChebyFit[data, len, Min[n, len-1]];
				answer =!= $Failed)
	];


FittingPolynomial[f_Function, n_][x_] := f[x];

Format[FittingPolynomial[f_Function, n_]] :=
	SequenceForm["FittingPolynomial[<>, ", Round[n], "]"];

Cheby[n_, x_?NumberQ] := ChebyshevT[n, x];

chebycoeff[Times[coeff_, cheby_Cheby]] := {coeff, cheby[[1]]};
chebycoeff[0.] = {0,0};

PolynomialFit::baddata = "`1` is not a list of ordered pairs of numbers."

ChebyFit[data_, len_, n_] :=
    Module[{i, x, t, a, b, c, chebylist, temp, prec},
	If[ n < 0, Return[$Failed]];
	If[VectorQ[data],
		If[ n == 0, Return[Function @@ {Apply[Plus,data]/len}]];
		t = (2x-len-1)/(len-1),
	    (* else *)
		a = Transpose[data];
		If[n == 0, Return[Function @@ {Apply[Plus, a[[2]]]/len}]];
		a = a[[1]];
		If[!VectorQ[a],
			Message[PolynomialFit::baddata, data];
			Return[$Failed]];
		b = Apply[Plus,a]/len;
		If[!NumberQ[b],
			Message[PolynomialFit::baddata, data];
			Return[$Failed]];
		a = 2 Sqrt[(a . a)/len - b b];
		t = (x-b)/a
		];
	prec = Precision[data];
	If[prec === Infinity, prec = $MachinePrecision];
	t = N[t, prec];
	chebylist = Fit[data, Table[Cheby[i, t], {i, 0, n}], x];
	If[Head[chebylist] =!= Plus, Return[$Failed]];

	(* form a simple list of coefficients including zeros *)
	chebylist = Map[chebycoeff, Apply[List, chebylist]];
	temp = Table[0, {n+1}];
	Do[temp[[chebylist[[i,2]]+1]] = chebylist[[i,1]],
						{i,Length[chebylist]}];

	Clear[a,b,c];	(* These variables will be used to smuggle information
			   inside of Function[ ] which is HoldAll. *)
	t = t /. x -> #;

	(* Construct the algorithm for evaluating the polynomial, but
		don't evaluate anything. *)
	Function[Module[{ta = 0, tb = 0, tc, tt = a, clist = b},
			(* use Clenshaw's algorithm *)
			Do[tc = 2 tt tb - ta + clist[[ti]];
				ta = tb;
				tb = tc,
				{ti, Round[c], Round[2], Round[-1]}];
			tt tb - ta + First[clist]]] /.
					{a -> t, b -> temp, c -> n+1}
	];

End[]  (* NumericalMath`PolynomialFit`Private` *)

EndPackage[] (* NumericalMath`PolynomialFit` *)

