#C10.txt, Feb. 24, 2014 Orthogonal polynomials in general
Help:=proc(): print(`GuessRF(L,n) , AM(M,P,x)`): 
print(`PnM(M,n,x) , PnMseq(M,n,x), WtToM(w,n,x,a,b), AnM(M,n) `):
end:

#old stuff
#GuessRF1(L,d,n): inputs a list of numbers (or expressions)
#L, a non-neg. integer d, and a (discrete) variable name
#n and tries to guess a rational function of n of degree
#d, such that L[i]-R(i)=0 for all i from 1 to nops(L)
GuessRF1:=proc(L,d,n) local R,a,b,i,j,eq,var:

if nops(L)<=2*d+6 then
 print(`Make list bigger`):
 RETURN(FAIL):
fi:

R:=(n^d+add(a[i]*n^i,i=0..d-1))/add(b[j]*n^j,j=0..d):

var:={seq(a[i],i=0..d-1), seq(b[j],j=0..d)}:

eq:={seq( numer(L[i]-subs(n=i,R)) =0,i=1..2*d+6)}:

var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
fi:

R:=subs(var , R):


if {seq( numer(L[i]-subs(n=i,R)),i=1..nops(L))}<>{0}  then
 print(`It worked up to`, 2*d+6, `terms, but that's life `):
 FAIL:
else
 R:
fi:


end:


#GuessRF(L,n): inputs a list of numbers (or expressions)
#L, and a (discrete) variable name
#n and tries to guess a rational function of n of degree
#<=(nops(L)-6)/2 such that L[i]-R(i)=0 for all i from 1 to nops(L)
GuessRF:=proc(L,n) local d, katie:

for d from 0 to (nops(L)-6)/2 do
  katie:=GuessRF1(L,d,n):

   if katie<>FAIL then
      RETURN(katie):
   fi:

od:

FAIL:

end:

#end old stuff

#AM(M,P,x): inputs a sequence of moments x^i->M[i+1]
#outputs the "integral" of polynomial P of x
AM:=proc(M,P,x) local i:

if nops(M)<degree(P,x)+1 then
 FAIL:
else
add(M[i+1]*coeff(P,x,i),i=0..degree(P,x)):
fi:

end:




#inputs a sequence of numbers "the moments" such
#that M[i+1] is what is mapped to x^i
#outputs the n-th degree MONIC polynomial in
#the sequence of OPs such that
#Umbra(P(n)*P(m))=0 if n <>m
PnM:=proc(M,n,x) local P,i,a,var,eq:
option remember:

P:=x^n+add(a[i]*x^i, i=0..n-1): 

var:={seq(a[i],i=0..n-1)}:

#AM (M,x^i*P,x)=0, i=0, ...n-1
eq:={ seq(AM(M,x^i*P,x)=0, i=0..n-1)}:

var:=solve(eq,var):

if var=NULL then
 FAIL:
else
subs(var,P):
fi:


end:


PnMseq:=proc(M,n,x) local i:
[seq(PnM(M,i,x),i=0..n)]:
end:

#WtToM(w,n,x,a,b): the moments of f(x)->int(w(x)*f(x),x=a..b);
WtToM:=proc(w,n,x,a,b) local i:
[seq(int(w*x^i,x=a..b), i=0..n)]:

end:


#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)^2)

#AnM(M,n): the first n terms of the A(n)
AnM:=proc(M,n) local i,andrew,x:
andrew:=PnMseq(M,n,x):

[seq(AM(M,x*andrew[i+1]^2,x)/AM(M,andrew[i+1]^2,x),i=1..n)]:
end: