from sage.all import *

def legendre_p(n, x, a = -1, b = 1):
    """
    Returns the ``n``th Legendre polynomial of the first kind over the
    interval [a, b] with respect to ``x``.

    When [a,b] is not [-1,1], we scale the standard Legendre
    polynomial (which is defined over [-1,1]) via an affine map. The
    resulting polynomials are still orthogonal and possess the
    property that `P(a) = P(b) = 1`.

    INPUT:

      * ``n`` -- The index of the polynomial.

      * ``x`` -- Either the variable to use as the independent
        variable in the polynomial, or a point at which to evaluate
        the polynomial.

      * ``a`` -- The "left" endpoint of the interval. Must be a real number.

      * ``b`` -- The "right" endpoint of the interval. Must be a real number.

    OUTPUT:

    If ``x`` is a variable, a polynomial (symbolic expression) will be
    returned. Otherwise, the value of the ``n``th polynomial at ``x``
    will be returned.

    EXAMPLES:

    Create the standard Legendre polynomials in `x`::

        sage: legendre_p(0,x)
        1
        sage: legendre_p(1,x)
        x

    Reuse the variable from a polynomial ring::
        sage: P.<t> = QQ[]
        sage: legendre_p(2,t).simplify_rational()
        3/2*t^2 - 1/2

    If ``x`` is a real number, the result should be as well::

        sage: legendre_p(3, 1.1)
        1.67750000000000

    Similarly for complex numbers::

        sage: legendre_p(3, 1 + I)
        7/2*I - 13/2

    Even matrices work::

        sage: legendre_p(3, MatrixSpace(ZZ, 2)([1, 2, -4, 7]))
        [-179  242]
        [-484  547]

    TESTS:

    We should agree with Maxima for all `n`::

        sage: eq = lambda k: bool(legendre_p(k,x) == legendre_P(k,x))
        sage: all([eq(k) for k in range(0,20) ]) # long time
        True

    We can evaluate the result of the zeroth polynomial::

        sage: f = legendre_p(0,x)
        sage: f(x=10)
        1

    We should have |P(a)| = |P(b)| = 1 for all a,b::

        sage: a = RR.random_element()
        sage: b = RR.random_element()
        sage: k = ZZ.random_element(20)
        sage: P = legendre_p(k, x, a, b)
        sage: abs(P(x=a)) # abs tol 1e-12
        1
        sage: abs(P(x=b)) # abs tol 1e-12
        1

    Two different polynomials should be orthogonal with respect to the
    inner product over `[a,b]`. Note that this test can fail if QQ is
    replaced with RR, since integrate() can return NaN::

        sage: a = QQ.random_element()
        sage: b = QQ.random_element()
        sage: j = ZZ.random_element(20)
        sage: k = j + 1
        sage: Pj = legendre_p(j, x, a, b)
        sage: Pk = legendre_p(k, x, a, b)
        sage: integrate(Pj*Pk, x, a, b) # abs tol 1e-12
        0

    The first few polynomials shifted to [0,1] are known to be::

        sage: p0 = 1
        sage: p1 = 2*x - 1
        sage: p2 = 6*x^2 - 6*x + 1
        sage: p3 = 20*x^3 - 30*x^2 + 12*x - 1
        sage: bool(legendre_p(0, x, 0, 1) == p0)
        True
        sage: bool(legendre_p(1, x, 0, 1) == p1)
        True
        sage: bool(legendre_p(2, x, 0, 1) == p2)
        True
        sage: bool(legendre_p(3, x, 0, 1) == p3)
        True

    The zeroth polynomial is an element of the ring that we're working
    in::

        sage: legendre_p(0, MatrixSpace(ZZ, 2)([1, 2, -4, 7]))
        [1 0]
        [0 1]

    """
    if not n in ZZ:
        raise TypeError('n must be a natural number')

    if n < 0:
        raise ValueError('n must be nonnegative')

    if not (a in RR and b in RR):
        raise TypeError('both `a` and `b` must be real numbers')

    if n == 0:
        # Easy base case, save time. Attempt to return a value in the
        # same field/ring as `x`.
        return x.parent()(1)

    # Even though we know a,b are real we use the symbolic ring. This
    # lets us return pretty expressions where possible.
    a = SR(a)
    b = SR(b)
    n = ZZ(n) # Ensure that 1/(2**n) is not integer division.
    dn = 1/(2**n)

    def phi(t):
        # This is an affine map from [a,b] into [-1,1] and so preserves
        # orthogonality.
        return (2 / (b-a))*t + 1 - (2*b)/(b-a)

    def c(m):
        return binomial(n,m)*binomial(n, n-m)

    def g(m):
        # As given in A&S, but with `x` replaced by `phi(x)`.
        return ( ((phi(x) - 1)**(n-m)) * (phi(x) + 1)**m )

    # From Abramowitz & Stegun, (22.3.2) with alpha = beta = 0.
    P = dn * sum([ c(m)*g(m) for m in range(0,n+1) ])

    return P