[SciPy-dev] Numerical Recipes (was tagging 0.7rc1 this weekend?)

[Fernando Perez]

Howdy,

On Tue, Dec 16, 2008 at 7:34 AM, Alan G Isaac <aisaac at american.edu> wrote:

>>> See Numerical Recipies in C, page 156 and
>>> Abramowitz and Stegun p. 774, 782
>
>
>
> Travis O. is listed as writing this orthogonal.py in 2000
> and doing bug fixes in 2003.  Note that the
> reference to Numerical Recipes is one of two
> references documenting an algebraic recursion
> relation and is not referring to an implementation.
>
> http://www.fizyka.umk.pl/nrbook/c4-5.pdf
>
> http://www-lsp.ujf-grenoble.fr/recherche/a3t2/a3t2a2/bahram/biblio/abramowitz_and_stegun/page_774.htm
>
> I do not see any problem with the code.
> I think this one should also be considered resolved.

I concur.  I also had a look at the Golub&Welsch reference pointed to
by Chuck, and I don't see any problems.  The basic algorithm may be
the same, but the actual implementation in those chapters of NR is
fairly different from what's in special/orthogonal.py.  The Scipy code
is structured as a bunch of routines that provide quite a few
different orthogonal polyinomial objects, while NR only lists code for
a few specific cases where the weights and abscissas are being
computed (gauleg, gaulag, gauher, gaujac), plus the gaucof() that
implements the Golub-Welsch algorithm for generic abscissas and
weights.

In orthogonal.py, the closest pairing is between  NR's gaucof() and
gen_roots_and_weights(), but their actual codes are fairly different
to the extent possible. That is, it's a short code that does a very
specific thing, but the scipy code does it in a very python/numpy
specific way, while the NR code is very 'C'-ish.

Other than the fact that these codes both use the same piece of
mathematics, I see no problem here whatsoever.  It might be worth
clarifying in the docstring that reads

"""
See Numerical Recipies in C, page 156 and
Abramowitz and Stegun p. 774, 782
"""

that the NR reference is to the formulas only,  not to any actual
implementation.  If nothing else, to prevent any questions in the
future...

Cheers,

f