Repository: scipy/scipy
File: scipy/special/_orthogonal.py (Python)

Code before:
    >>> import numpy as np
    >>> from scipy.special import sh_legendre, legendre
    >>> from scipy.integrate import trapezoid
    >>> x = np.arange(0.0, 1.0, 0.01)
    >>> n = 3
    >>> np.allclose(sh_legendre(n)(x), legendre(n)(2*x - 1))
    True

    The polynomials :math:`P_n^*` satisfy a recurrence
    relation obtained by the change of variables
    :math:`t = 2x - 1` in the standard Legendre recurrence:
    
    .. math::
    
        (n+1) P_{n+1}^*(x) = (2n+1)(2x-1)\,P_n^*(x) - n\,P_{n-1}^*(x).

    This can be easily checked on :math:`[0, 1]`
    for :math:`n = 3`:

    >>> x = np.linspace(0.0, 1.0, 101)
    >>> lhs = (3 + 1) * sh_legendre(4)(x)
    >>> rhs = (
    ...     (2*3 + 1) * (2*x - 1) * sh_legendre(3)(x)
    ...     - 3 * sh_legendre(2)(x)
    ... )
    >>> np.allclose(lhs, rhs)
    True

    Orthogonality over :math:`[0,1]` with weight 1 can be
    checked numerically; for example, :math:`P_2^*`
    is orthogonal to :math:`P_3^*`:

    >>> x = np.linspace(0.0, 1.0, 400)
    >>> y = sh_legendre(2)(x) * sh_legendre(3)(x)
    >>> np.isclose(trapezoid(y, x), 0.0, atol=1e-12)
    True

    See Also
    --------
    scipy.special.legendre
    scipy.special.roots_sh_legendre
    scipy.integrate.trapezoid
    """
    if n < 0:
        raise ValueError("n must be nonnegative.")

    def wfunc(x):
        return 0.0 * x + 1.0
    if n == 0:
        return orthopoly1d([], [], 1.0, 1.0, wfunc, (0, 1), monic,
                           lambda x: _ufuncs.eval_sh_legendre(n, x))

Review comment (suggestion):
```suggestion
    >>> n = 3
    >>> x = np.linspace(0.0, 1.0, 101)
    >>> lhs = (n + 1) * sh_legendre(n + 1)(x)
    >>> rhs = (
    ...     (2*n + 1) * (2*x - 1) * sh_legendre(n)(x)
    ...     - n * sh_legendre(n - 1)(x)
    ... )
```
final suggestion

Code after:
    >>> import numpy as np
    >>> from scipy.special import sh_legendre, legendre
    >>> from scipy.integrate import trapezoid
    >>> x = np.arange(0.0, 1.0, 0.01)
    >>> n = 3
    >>> np.allclose(sh_legendre(n)(x), legendre(n)(2*x - 1))
    True

    The polynomials :math:`P_n^*` satisfy a recurrence
    relation obtained by the change of variables
    :math:`t = 2x - 1` in the standard Legendre recurrence:
    
    .. math::
    
        (n+1) P_{n+1}^*(x) = (2n+1)(2x-1)\,P_n^*(x) - n\,P_{n-1}^*(x).

    This can be easily checked on :math:`[0, 1]`
    for :math:`n = 3`:

    >>> n = 3
    >>> x = np.linspace(0.0, 1.0, 101)
    >>> lhs = (n + 1) * sh_legendre(n + 1)(x)
    >>> rhs = (
    ...     (2*n + 1) * (2*x - 1) * sh_legendre(n)(x)
    ...     - n * sh_legendre(n - 1)(x)
    ... )
    >>> np.allclose(lhs, rhs)
    True

    Orthogonality over :math:`[0,1]` with weight 1 can be
    checked numerically; for example, :math:`P_2^*`
    is orthogonal to :math:`P_3^*`:

    >>> x = np.linspace(0.0, 1.0, 400)
    >>> y = sh_legendre(2)(x) * sh_legendre(3)(x)
    >>> np.isclose(trapezoid(y, x), 0.0, atol=1e-12)
    True

    See Also
    --------
    scipy.special.legendre
    scipy.special.roots_sh_legendre
    """
    if n < 0:
        raise ValueError("n must be nonnegative.")

    def wfunc(x):
        return 0.0 * x + 1.0
    if n == 0:
        return orthopoly1d([], [], 1.0, 1.0, wfunc, (0, 1), monic,
                           lambda x: _ufuncs.eval_sh_legendre(n, x))