#hw9.txt Ross Berkowitz 2/24/14
#Posting Permission Granted

Help:=proc():
print(`ApplyUmbra(P,x,U,m), PnG(n,x,U,m)`):
end:

#################################  1  ################################
#ApplyUmbra(P,x,U,m) applies the linear map given by the action x^m->U(m)
#to a given polynomial P in the variable x

ApplyUmbra:=proc(P,x,U,m) local i:
add(coeff(P,x,i)*subs(m=i,U),i=0..degree(P,x)):
end:


################################  2  ################################

#PnG(n,x,U,m) inputs a degree n, variable x, umbra U and discrete variable m
#and outputs the unique polynomial P such that U(P*x^i)=0 for all i from 0 to n

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

################################  3   ################################


#GuessCoeff(n,u,U,m) inputs a discrete variable n, a coefficient i and the
#umbra U in the discrete variable m and guess a rational function (in n) describing
#the coefficient of x^(n-i) in the "orthogonal" polynomial PnG(n,x,U,m)

GuessCoeff:=proc(n,i,U,m) local bound,j,hope,L:
option remember:
bound:=1:
hope:=FAIL:
while bound<6 and hope=FAIL do
	L:=[seq(coeff(PnG(j,x,U,m),x,j-i),j=i..2^bound)]:
	hope:=GuessRF(L,n):
	bound:=bound+1:
od:
hope:
end:


#seq(GuessCoeff(n,i,1/(m+1),m),i=0..7)=
#1,
# -(1/2)*n, (n^3+n^2)/(4+8*n),
# (n^4+2*n^2+3*n^3)/(-72-48*n),
# (n^6+6*n^2+17*n^3+17*n^4+7*n^5)/(1440+1536*n+384*n^2), #(n^7+24*n^2+74*n^3+85*n^4+45*n^5+11*n^6)/(-33600-23040*n-3840*n^2),
#(n^9+240*n^2+908*n^3+1392*n^4+1119*n^5+510*n^6+132*n^7+18*n^8)/(1814400+1647360*n+483840*#n^2+46080*n^3),
 #(n^10+1440*n^2+5688*n^3+9260*n^4+8106*n^5+4179*n^6+1302*n^7+240*n^8+24*n^9)/(-55883520-3#8545920*n-8709120*n^2-645120*n^3)

#seq(GuessCoeff(n,i,1/(m+3),m),i=0..7)=
#1,
# (n^2+2*n)/(-2-2*n),
# (n^3+3*n+4*n^2)/(12+8*n),
# (n^4+8*n+14*n^2+7*n^3)/(-120-48*n), #(n^6+60*n+152*n^2+143*n^3+63*n^4+13*n^5)/(3360+2304*n+384*n^2),
# (n^7+288*n+792*n^2+844*n^3+450*n^4+127*n^5+18*n^6)/(-60480-30720*n-3840*n^2), #(n^9+5040*n+16428*n^2+22124*n^3+16203*n^4+7089*n^5+1902*n^6+306*n^7+27*n^8)/(3991680+2753#280*n+622080*n^2+46080*n^3), #(n^10+34560*n+117792*n^2+168456*n^3+133636*n^4+65086*n^5+20237*n^6+4024*n^7+494*n^8+34*n^#9)/(-103783680-57899520*n-10644480*n^2-645120*n^3)
#
#GuessCoeff(n,i,eval(m!),m) did not terminate successfully.








################################ CLASS STUFF  ################################

#C6.txt: Feb. 10, 2014, the joy of guessing
Help6:=proc():  print(`Hn(n) , HnS(N), Hnr(n,r), HnrS(r,N)`):
print( ` Ratio1(L) , GuessRF1(L,d,n)`): 
print(`GuessRF(L,n) `):
end:

with(linalg):

#Hnr(n,r): the determinant of the generalized Hilbert n by n matrix
#whose (1,1) entry is not 1 by 1/(r+1)
Hnr:=proc(n,r) local i,j:
det([seq([seq(1/(i+j-1+r),j=1..n)],i=1..n)]):
end:

#Hn(n): the determinant of the famous Hilbert n by n matrix
Hn:=proc(n) local i,j:
det([seq([seq(1/(i+j-1),j=1..n)],i=1..n)]):
end:

#The first N terms of Hn, starting at n=1
HnS:=proc(N) local n:

[seq(Hn(n), n=1..N)]:

end:

#The first N terms of Hn, starting at n=1
HnrS:=proc(r,N) local n:

[seq(Hnr(n,r), n=1..N)]:

end:
#Ratio1(L): inputs a list L and outputs the list of L[i+1]/L[i]
Ratio1:=proc(L) local i:

[seq(L[i+1]/L[i],i=1..nops(L)-1)];
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: