from sage.all import *
from sage.symbolic.expression import is_Expression

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

    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.

    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.

    TESTS:

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

        sage: eq = lambda k: bool(standard_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 = standard_legendre_p(0,x)
        sage: f(x=10)
        1

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

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

    # Ensures that 1/(2**n) is not integer division.
    n = ZZ(n)
    dn = 1/(2**n)

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

    def g(m):
        return ( ((x - 1)**(n-m)) * (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) ])

    # If `x` is a symbolic expression, we want to return a symbolic
    # expression (even if that expression is e.g. `1`).
    if is_Expression(x):
        P = SR(P)

    return P


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

    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.

    TESTS:

    We agree with ``standard_legendre_p`` when `a=b=1`::

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

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

        sage: a = QQ.random_element()
        sage: b = QQ.random_element()
        sage: k = ZZ.random_element(20)
        sage: P = legendre_p(k, x, a, b)
        sage: bool(abs(P(x=a)) == 1)
        True
        sage: bool(abs(P(x=b)) == 1)
        True

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

    a = RR(a)
    b = RR(b)
    t = SR.symbol('t')
    P = standard_legendre_p(n, t)

    # This is an affine map from [a,b] into [-1,1] and so preserves
    # orthogonality. If we define this with 'def' instead of a lambda,
    # Python segfaults as we evaluate P.
    phi = lambda y: (2 / (b-a))*y + 1 - (2*b)/(b-a)

    P_tilde = P(t = phi(x))

    return P_tilde