<commit_before>#!/usr/bin/env python
#
# Author:  <NAME> 2000
# Updated Sep. 2003 (fixed bugs --- tested to be accurate)

"""
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

For a more detailed discussion see Numerical Recipies in C (p. 156) and
Abramowitz and Stegun (p. 774, 782); while our implementation differs, the
algorithm and formulas are similar.

Functions:

  gen_roots_and_weights  -- Generic roots and weights.
  j_roots                -- Jacobi
  js_roots               -- Shifted Jacobi
  la_roots               -- Generalized Laguerre
  h_roots                -- Hermite
  he_roots               -- Hermite (unit-variance)
  cg_roots               -- Ultraspherical (Gegenbauer)
  t_roots                -- Chebyshev of the first kind
  u_roots                -- Chebyshev of the second kind
  c_roots                -- Chebyshev of the first kind ([-2,2] interval)
  s_roots                -- Chebyshev of the second kind ([-2,2] interval)
  ts_roots               -- Shifted Chebyshev of the first kind.
  us_roots               -- Shifted Chebyshev of the second kind.
  p_roots                -- Legendre
  ps_roots               -- Shifted Legendre
  l_roots                -- Laguerre
"""

# Scipy imports.
import numpy as np
from numpy import all, any, exp, inf, pi, sqrt
from numpy.dual import eig

# Local imports.
import _cephes as cephes
_gam = cephes.gamma

def poch(z,m):
    """Pochhammer symbol (z)_m = (z)(z+1)....(z+m-1) = gamma(z+m)/gamma(z)"""
    return _gam(z+m) / _gam(z)

class orthopoly1d(np.poly1d):
    def __init__(self, roots, weights=None, hn=1.0, kn=1.0, wfunc=None, limits=None, monic=0):
        np.poly1d.__init__(self, roots, r=1)
        equiv_weights = [weights[k] / wfunc(roots[k]) for k in range(len(roots))]
        self.__dict__['weights'] = np.array(zip(roots,weights,equiv_weights))
        self.__dict__['weight_func'] = wfunc
        self.__dict__['limits'] = limits
        mu = sqrt(hn)
        if monic:
            mu = mu / abs(kn)
            kn = 1.0
        self.__dict__['normcoef'] = mu
        self.__dict__['coeffs'] *= kn
<commit_msg>Add bibliographical references contributed by <NAME>saac after discussion of origins of the algorithm on the mailing list.<commit_after>#!/usr/bin/env python
#
# Author:  <NAME>00
# Updated Sep. 2003 (fixed bugs --- tested to be accurate)

"""
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

For the mathematical background, see [golub.welsch-1969-mathcomp]_ and
[abramowitz.stegun-1965]_.

Functions:

  gen_roots_and_weights  -- Generic roots and weights.
  j_roots                -- Jacobi
  js_roots               -- Shifted Jacobi
  la_roots               -- Generalized Laguerre
  h_roots                -- Hermite
  he_roots               -- Hermite (unit-variance)
  cg_roots               -- Ultraspherical (Gegenbauer)
  t_roots                -- Chebyshev of the first kind
  u_roots                -- Chebyshev of the second kind
  c_roots                -- Chebyshev of the first kind ([-2,2] interval)
  s_roots                -- Chebyshev of the second kind ([-2,2] interval)
  ts_roots               -- Shifted Chebyshev of the first kind.
  us_roots               -- Shifted Chebyshev of the second kind.
  p_roots                -- Legendre
  ps_roots               -- Shifted Legendre
  l_roots                -- Laguerre


.. [golub.welsch-1969-mathcomp]
   Golub, <NAME>, and <NAME>. 1969. Calculation of Gauss
   Quadrature Rules. *Mathematics of Computation* 23, 221-230+s1--s10.

.. [abramowitz.stegun-1965]
   Abramowitz, Milton, and <NAME>. (1965) *Handbook of
   Mathematical Functions: with Formulas, Graphs, and Mathematical
   Tables*. Gaithersburg, MD: National Bureau of Standards.
   http://www.math.sfu.ca/~cbm/aands/
"""

# Scipy imports.
import numpy as np
from numpy import all, any, exp, inf, pi, sqrt
from numpy.dual import eig

# Local imports.
import _cephes as cephes
_gam = cephes.gamma

def poch(z,m):
    """Pochhammer symbol (z)_m = (z)(z+1)....(z+m-1) = gamma(z+m)/gamma(z)"""
    return _gam(z+m) / _gam(z)

class orthopoly1d(np.poly1d):
    def __init__(self, roots, weights=None, hn=1.0, kn=1.0, wfunc=None, limits=None, monic=0):
        np.poly1d.__init__(self, roots, r=1)
        equiv_weights = [weights[k] / wfunc(roots[k]) for k in range(len(roots))]
        self.__dict__['weights'] = np.array(zip(roots,weights,equiv_weights))
        self.__dict__['weight_func'] = wfunc
        self.__dict__['limits'] = limits
        mu = sqrt(hn)
        if monic:
            mu = mu / abs(kn)
            kn = 1.0
        self.__dict__['normcoef'] = mu
        self.__dict__['coeffs'] *= kn

{'fraction_non_alphanumeric': 0.10562685093780849, 'fraction_numerical': 0.021388614675880224, 'mean_word_length': 2.426719278466742, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2199903', 'n_tokens_mistral': 2234, 'n_tokens_neox': 2016, 'n_words': 643}