[scipy.special] Fast and stabler evaluation of orthogonal polynomials

[scipy.special] Fast and stabler evaluation of orthogonal polynomials


Rationale
---------

Often, one needs just values of orthogonal polynomials at a couple of points.
This calls for vectorized functions -- the routines currently in orthogonal.py
compute roots first, which is quite slow.

Sometimes high orders are also needed, so the evaluation functions need to
be sufficiently stable numerically.

Moreover, it would be useful to integrate the evaluation functions to the
current interface.


This patch
----------

- Adds vectorized eval_* routines for evaluating orthogonal polynomials.

  They are either written using stable recurrences in Cython, or related to
  hypergeometric functions.

- Hooks these evaluation routines to the orthopoly1d class, so that __call__
  uses them.

  This gains mainly numerical stability.

  orthopoly1d lose their identity on arithmetic, so after any arithmetic,
  evaluation is done via the unstable numpy.poly1d route.

- Makes *all* orthogonal routines return orthopoly1d's. I believe the
  normalisation metadata etc. are still handled properly.

- There are tests comparing the numpy.poly1d route for low polynomial
  orders.


TODO
----

- I haven't really carefully checked that eg. eval_jacobi is *really*
  numerically stabler. Quick checks seem to show that it is, but
  more tests are needed.


Notes
-----

- Many of the polynomials are related to hyp2f1. I fixed recently several bugs
  in it that hit hard the parameter ranges useful for orthogonal polynomials,
  so if you want to try this out, grab a SVN version of Scipy.

[Anne, 2009-10-26 01:20:19]

Looks good. I'd say go ahead and add it.

In the long run, I think we need something more comprehensive; as it stands a
user must remember to do a*p(x)+q(x) rather than (a*p+q)(x) or their numerical
stability will mysteriously disappear. 

My favored approach would be to allow any family of orthogonal polynomials to be
used as a basis, since (I think) there is an efficient and stable recurrence to
evaluate linear combinations of orthogonal polynomials. But this is a much
bigger piece of code, and the subject of much debate on the mailing list. This
patch is a big improvement that's ready to go right now.