[Arjen Markus] (3 april 2004) Knowing orthogonal polynomials like the Legendre family and others, and having played a bit with discrete Fourier transforms (see [Discrete Fourier Transform]) I thought I would try an analoguous technique polynomials.

There is one practical difference with discrete Fourier transforms: it is difficult (but not impossible) to write down closed formulae for the polynomials to be used. I started with:

   * The nodes or auxiliary points are x = k/N (where k = -m, -m+1, ..., m-1, m, N=2m+1,
and similar - but different for even N)
   * The poynomials must be orthogonal in the sense that the sum
       m
      Sum   Pi(x)*Pj(x) = 0,   if i != j
      k=-m
   * The first two are simple:
    P0(x) = 1
    P1(x) = x
   * The third becomes troublesome already because you have to distinguish between even and odd N.

So, I decided to give up on this analytic approach and instead use a practical (numerical) one:
   * I try to fit N data points
   * The nodes for the polynomials are x = k, k=1,N
   * The polynomials that I want to fit with are determined "on the fly", using a simple Gram-Schmidt orthogonalisation process and subsequent normalisation.
   * Due to the orthogonality the fitting of the data becomes almost trivial, just as with the discrete Fourier transform.

(Improvements:
   * Reduce the overhead of double computations - certain intermediate results are computed more than once
   * The toString procedure does not take care of minus signs - I overlooked that possibility and was too lazy to repair that. Probably the best solution is to modify the string afterwards.)

----