#hw9.txt, 2/23/14, Anthony Zaleski

with(PolynomialTools):

Help:=proc()
print(`ApplyUmbra(P,x,U,m),PnG(n,x,U,m),GuessCoeff(n,i,U,m)`):
print(`GuessRF(L,n),Pn(n,x)`):
end:

##########Problem 1: Evaluating linear functionals given action on basis##########

#ApplyUmbra(P,x,U,m) plugs the polynomial P(x)
#into the linear functional sending x^m -> U(m).

ApplyUmbra:=proc(P,x,U,m) local L,i:

L:=CoefficientList(P,x):
add(L[i]*subs(m=i-1,U),i=1..nops(L)):

end:

##########Problem 2: Generalizing Pn for arbitrary functionals##########

#PnG(n,x,U,m) finds the unique monic poly P(x) 
#of degree n s.t. U(P*x^i)=0 for i=0..n-1.

PnG:=proc(n,x,U,m) local P,i,a,eq,var:

P:=add(a[i]*x^i,i=0..n-1)+x^n:
eq:={seq(ApplyUmbra(P*x^i,x,U,m),i=0..n-1)}:
var:={seq(a[i],i=0..n-1)}:

sort(subs(solve(eq,var),P)):

end:

#It's easy to check Pn(n,x)=PnG(n,x,1/(m+1),m) for n=1..10.

##########Problem 3: GuessCoeff(n,i,U,m)##########

#GuessCoeff(n,i,U,m) inputs a discrete variable n, 
#a non-negative integer i, an expression U 
#(describing an umbra) in m, and a variable name m, 
#and guesses a rational function in n for the 
#coefficient of x^(n-i) (for n>=i) in PnG(n,x,U,m). 

GuessCoeff:=proc(n,i,U,m) local N,L,nn,x,Exp:

N:=30: 
L:=[0$N]:

for nn from i to i+N-1 do
	L[nn-i+1]:=coeff(PnG(nn,x,U,m),x,nn-i):
od:

Exp:=GuessRF(L,x): # note x=1 corresponds to n=i
subs(x=n-i+1,Exp):

end:



#Here is GuessCoeff in action!  Beware, for VERY long expressions follow!



(*
for i from 0 to 7 do GuessCoeff(n,i,1/(m+1),m); print(`_______________________`); od;
1
_______________________
-(1/2)*n
_______________________
((n-1)^3+(n-1)^2)/(8*n-4)
_______________________
((n-2)^4+3*(n-2)^3+2*(n-2)^2)/(-48*n+24)
_______________________
((n-3)^6+7*(n-3)^5+17*(n-3)^4+17*(n-3)^3+6*(n-3)^2)/(384*(n-3)^2+1536*n-3168)
_______________________
((n-4)^7+11*(n-4)^6+45*(n-4)^5+85*(n-4)^4+74*(n-4)^3+24*(n-4)^2)/(-3840*(n-4)^2-23040*n+58560)
_______________________
((n-5)^9+18*(n-5)^8+132*(n-5)^7+510*(n-5)^6+1119*(n-5)^5+1392*(n-5)^4+908*(n-5)^3+240*(n-5)^2)/(46080*(n-5)^3+483840*(n-5)^2+1647360*n-6422400)
_______________________
((n-6)^10+24*(n-6)^9+240*(n-6)^8+1302*(n-6)^7+4179*(n-6)^6+8106*(n-6)^5+9260*(n-6)^4+5688*(n-6)^3+1440*(n-6)^2)/(-645120*(n-6)^3-8709120*(n-6)^2-38545920*n+175392000)
_______



for i from 0 to 7 do GuessCoeff(n,i,1/(m+3),m); print(`_______________________`); od;
1
_______________________
(n^2+2*n)/(-2*n-2)
_______________________
((n-1)^3+4*(n-1)^2+3*n-3)/(8*n+4)
_______________________
((n-2)^4+7*(n-2)^3+14*(n-2)^2+8*n-16)/(-48*n-24)
_______________________
((n-3)^6+13*(n-3)^5+63*(n-3)^4+143*(n-3)^3+152*(n-3)^2+60*n-180)/(384*(n-3)^2+2304*n-3552)
_______________________
((n-4)^7+18*(n-4)^6+127*(n-4)^5+450*(n-4)^4+844*(n-4)^3+792*(n-4)^2+288*n-1152)/(-3840*(n-4)^2-30720*n+62400)
_______________________
((n-5)^9+27*(n-5)^8+306*(n-5)^7+1902*(n-5)^6+7089*(n-5)^5+16203*(n-5)^4+22124*(n-5)^3+16428*(n-5)^2+5040*n-25200)/(46080*(n-5)^3+622080*(n-5)^2+2753280*n-9774720)
_______________________
((n-6)^10+34*(n-6)^9+494*(n-6)^8+4024*(n-6)^7+20237*(n-6)^6+65086*(n-6)^5+133636*(n-6)^4+168456*(n-6)^3+117792*(n-6)^2+34560*n-207360)/(-645120*(n-6)^3-10644480*(n-6)^2-57899520*n+243613440)
_______________________



*)


#Note: Using U=m! took too long.







###########Relevant Old Functions##########

#Pn(n,x): the unique momic polynomial of x of degree n
#such that int(Pn(n,x)*x^i,x=0..1)=0 for i=0..n-1
Pn:=proc(n,x) local eq,var,i,P,a:
option remember:
P:=x^n+add(a[i]*x^i,i=0..n-1):

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

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

sort(subs(solve(eq,var),P)):
end:


#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 of old stuff###########