#Sowmya Srinivasan hw10.txt, 24 Feb 2014 #BnM(M,n) finds the first n terms of the B(n) sequence BnM:=proc(M,n) local i,s,x: s:=PnMseq(M,n,x): [seq(AM(M,x*s[i+1]*s[i],x)/AM(M,s[i]^2,x),i=1..n)] : end: # for the Legendre polynomials, A(n) = 2(n^2+2n+1)/(4n^2+8n+3) # and B(n) = 1/4 n(n+1)/(4n^2 +4n+1) # for Leguerre polynomials, A(n) = 2n+2 and B(n) = n^2+n # and for M=[seq(1/(i+a),i=1..1000)] we have #A(n) = (a^2+2an+2n^2+a+2n)/(a^2+4an+4n^2+2a+4n) and # B(n) = (n^2(a^2+2an+n^2))/(a^4+8a^3n+24a^2 n^2 +32a n^3 +16n^4-a^2-4an-4n^2) PnFast:=proc(A,B,n,n1,x,a0) local P,i: option remember : P[0]:=1 : P[1]:=x+a0 : for i from 2 to n1 do P[i]:= expand((x-subs(n=i-1,A))*P[i-1]-subs(n=i-1,B)*P[i-2] ) : od: P[n1]: end: #My computer seemed too slow to test for the 1000th (or even 50th) term #of either of the three moment sequences at first. Then I changed P[i] so #that it was expanding the output, and then the procedure was able to find P_1000