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

[Jarrod Millman]

Given the discussion about Numerical Recipes derived-code, I did a
little grepping and found 7 instances where code is referencing
Numerical Recipes.  Before branching for 0.7, I would like to see this
resolved.  I am not sure what we should do, but the license for using
Numerical Recipes derived code clearly does not permit us to include
it in scipy:
  http://www.nr.com/licenses/redistribute.html

Below I include all the instances I could easily find.  Thoughts?

--------------------------------------------------

In scipy/optimize/zeros.py:

def ridder(f, a, b, args=(),
           xtol=_xtol, rtol=_rtol, maxiter=_iter,
           full_output=False, disp=False):
    """Find root of f in [a,b]

    Ridder routine to find a zero of the function
    f between the arguments a and b. f(a) and f(b) can not
    have the same signs. Faster than bisection, but not
    generaly as fast as the brent rountines. A description
    may be found in a recent edition of Numerical Recipes.
    The routine here is modified a bit to be more careful
    of tolerance.

def brentq(f, a, b, args=(),
           xtol=_xtol, rtol=_rtol, maxiter=_iter,
           full_output=False, disp=False):
    """Find root of f in [a,b]

    The classic Brent routine to find a zero of the function
    f between the arguments a and b. f(a) and f(b) can not
    have the same signs. Generally the best of the routines here.
    It is a safe version of the secant method that uses inverse
    quadratic extrapolation. The version here is a slight
    modification that uses a different formula in the extrapolation
    step. A description may be found in Numerical Recipes, but the
    code here is probably easier to understand.

In scipy/signal/medianfilter.c:

/*
 *  This Quickselect routine is based on the algorithm described in
 *  "Numerical recipies in C", Second Edition,
 *  Cambridge University Press, 1992, Section 8.5, ISBN 0-521-43108-5
 */

In  scipy/special/orthogonal.py:

"""
A collection of functions to find the weights and abscissas for
Gaussian Quadrature.

These calculations are done by finding the eigenvalues of a
tridiagonal matrix whose entries are dependent on the coefficients
in the recursion formula for the orthogonal polynomials with the
corresponding weighting function over the interval.

Many recursion relations for orthogonal polynomials are given:

a1n f_n+1 (x) = (a2n + a3n x ) f_n (x) - a4n f_n-1 (x)

The recursion relation of interest is

P_n+1 (x) = (x - A_n) P_n (x) - B_n P_n-1 (x)

where P has a different normalization than f.

The coefficients can be found as:

A_n = -a2n / a3n

B_n = ( a4n / a3n sqrt(h_n-1 / h_n))**2

     where
             h_n = int_a^b w(x) f_n(x)^2
assume:
P_0(x) = 1
P_-1(x) == 0

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

In scipy/stats/stats.py:

def ttest_ind(a, b, axis=0):
    """Calculates the t-obtained T-test on TWO INDEPENDENT samples of scores
    a, and b.  From Numerical Recipies, p.483. Axis can equal None (ravel
    array first), or an integer (the axis over which to operate on a and b).

def ttest_rel(a,b,axis=None):
    """Calculates the t-obtained T-test on TWO RELATED samples of scores, a
    and b.  From Numerical Recipies, p.483. Axis can equal None (ravel array
    first), or an integer (the axis over which to operate on a and b).

def ks_2samp(data1, data2):
    """ Computes the Kolmogorov-Smirnof statistic on 2 samples.  Modified
    from Numerical Recipies in C, page 493.  Returns KS D-value, prob.  Not
    ufunc- like.