#Nathan Fox
#Homework 10
#I give permission for this work to be posted online

#Read procedures from class
read(`C10.txt`):

Help:=proc():
    print(` BnM(M, n) , PnFast(A, B, n, n1, x, a0) `):
end:

##PROBLEM 1##

#x*Pn(x)=P_(n+1)+A(n)*P_(n)+B(n)*P_(n-1)
#Int(x*Pn(x)^2)=A(n)*Int(Pn(x)^2), so
#A(n)=Int(x*Pn(x)^2)/Int(Pn(x)^2)
#B(n)=Int(x*Pn(x)*P_(n-1)(x))/Int(P_(n-1)(x)^2)

#BnM(M,n): the first n terms of the B(n)
BnM:=proc(M, n) local i, intseq, x:
    intseq:=PnMseq(M, n, x):
    return [seq(AM(M,x*intseq[i+1]*intseq[i],x)/
        AM(M,intseq[i]^2,x), i=1..n)]:
end:

##PROBLEM 2##

#Legendre:
#A(n)=1/2
#B(n)=n^2/(4*(2*n-1)*(2*n+1))

#Laguerre:
#A(n)=2*n+2
#B(n)=n*(n+1)

#Third:
#A(n)=(2*n^2+2*(a+1)*n+(a^2+a))/((2*n+a+2)*(2*n+a))
#B(n)=(n^2*(a+n)^2)/((2*n+a+1)*(2*n+a-1)*(2*n+a)^2)

##PROBLEM 3##

#PnFast(A,B,n,n1,x,a0): Inputs expressions A, B in n, a positive
#integer n1, and a variable name x and outputs the n1-th term in
#the sequence of polynomials satisfying the recurrence 
#x*Pn(x)= Pn+1(x)+A(n)*Pn(x)+B(n)*Pn-1(x) 
#under the initial conditions 
#P0(x)=1   P1(x)=x+a0
PnFast:=proc(A, B, n, n1, x, a0)
option remember:
    if n1 = 0 then
        return 1:
    elif n1 = 1 then
        return x+a0:
    else
        return simplify((x-subs(n=n1-1, A))*
            PnFast(A, B, n, n1-1, x, a0) -
            subs(n=n1-1, B)*PnFast(A, B, n, n1-2, x, a0)):
    fi:
end:

##PROBLEM 4##

#Legendre:
#PnFast(1/2, n^2/(4*(2*n-1)*(2*n+1)), n, 1000, x, 1/2);
#Done; too ridiculous to put here

#Laguerre:
#PnFast(2*n+2, n*(n+1), n, 1000, x, 2);
#Done; too ridiculous to put here

#Third:
#PnFast((2*n^2+2*(a+1)*n+(a^2+a))/((2*n+a+2)*(2*n+a)), (n^2*(a+n)^2)/((2*n+a+1)*(2*n+a-1)*(2*n+a)^2), n, 1000, x, 2);
#Took too long