using Polynomials, PyPlot, QuadGK, LinearAlgebra

# compute the dot product ⟨p,q⟩ = ∫p(x)q(x) on [-1,1]
polydot(p::Polynomial, q::Polynomial) = integrate(p*q, -1,1)

# force IJulia to display as LaTeX rather than HTML
Base.showable(::MIME"text/html", ::Polynomial) = false

Polynomial([3,4,5,6])

Polynomial([3,4,5,6])^2

p0 = a0 = Polynomial([1//1])

a1 = Polynomial([0, 1//1])

p1 = a1 - p0 * polydot(p0, a1) // polydot(p0, p0)
p1 = p1 / p1(1)

polydot(p0, a1)

a2 = Polynomial([0, 0, 1//1])

p2 = a2 - p0 * polydot(p0, a2) // polydot(p0, p0) -
          p1 * polydot(p1, a2) // polydot(p1, p1)
p2 = p2 / p2(1)

function legendre_gramschmidt(n)
    legendre = [Polynomial([1//1])]
    for i = 1:n
        p = Polynomial([k == i ? 1//1 : 0//1 for k=0:i])
        for q in legendre
            p = p - q * (polydot(q, p) // polydot(q,q))
        end
        push!(legendre, p / p(1))
    end
    return legendre
end

L = legendre_gramschmidt(5)

display.(L);

leg = []
x = range(-1, 1, length=300)
for i in eachindex(L)
    plot(x, L[i].(x), "-")
    push!(leg, L"P_{%$(i-1)}(x)")
end
plot(x, x * 0, "k--")
legend(leg)
xlabel(L"x")
title("Legendre polynomials")

p = Polynomial([1,3,4,7,2,5])

α = [polydot(q,p)/polydot(q,q) for q in L]

sum(α .* L) # α[1]*L[1] + α[2]*L[2] + ... + α[6]*L[6]

sum(α .* L) - p

polydot(p::Polynomial, f::Function) = quadgk(x -> p(x)*f(x), -1,1, atol=1e-13, rtol=1e-11)[1]

coeffs = [polydot(p,exp)/polydot(p,p) for p in L]

p = sum(coeffs .* L)

plot(x, exp.(x), "r-")
plot(x, p.(x), "b--")
legend([L"e^x", L"fit $p(x)$"])
xlabel(L"x")
title(L"fitting $e^x$ to a degree-5 polynomial")

plot(x, exp.(x), "k--")
for n = 1:4
    plot(x, sum([polydot(p,exp)/polydot(p,p) for p in L[1:n]] .* L[1:n]).(x), "-")
end
legend([L"e^x", ["degree $i" for i=0:3]...])
xlabel(L"x")
title(L"fitting $e^x$ to polynomials")

t(x) = 1 - abs(x)

L16 = legendre_gramschmidt(16) # compute a few more terms

plot(x, t.(x), "k--")
N = [1:2:16;]
for n in N
    plot(x, sum([polydot(p,t)/polydot(p,p) for p in L16[1:n]] .* L16[1:n]).(x), "-")
end
legend([L"t(x)", ["degree $i" for i in N .- 1]...])
xlabel(L"x")
title("fitting a triangle-shape function to polynomials")