#C14.txt: March 11, 2019, The P-recursive  (holonomic) "ansatz"

Help:=proc(): print(` SeqToP1(L,n,k,d), SeqToP1([seq(i!,i=1..20)],n,1,1); , SeqToP(L,n,MaxC), `):
print(`EstAsy(L,n), PtoSeq(P,N) , Asy(rec) `): end:

##from C12.txt
#C12.txt: March 4, 2019, (virtual class), Experimental Mathematics class, Math640 (Rutgers University)(Dr. Z.)





#Recall that a sequence {a(n), n=0..infinity} is C-finite (or order k) if there exist  NUMBERS c1, ..., ck such that


#a(n)=c1*a(n-1)+ ...+ ck*a(n-k)

#One you know c1, c2, ..., ck, and the k-1 initial values a(0), ..., a(k-1), the sequence is uniquely determined
#Our data structure for C-finite sequence is 
#C=[INI,Rec] where
#INI=[a(0),..., a(k-1)] and Rec=[c1,.., ck]. See
#http://sites.math.rutgers.edu/~zeilberg/EM19/C10.txt
#The imporant procedures were:

#CtoSeq(C,N) inputs a C-finite sequence (in our data-structure) and outputs the first N+1 terms of the sequence
#SeqToC(L): that inputs a list L and tries to find a linear-recurrence with CONSTANT coefficients satisfied by it
#(with enough `safey data'). Of course, if the list is not long enough, or the sequence is not C-finite it would return FAIL.

#The  NATURAL generalization of C-finite sequences is P-recursive sequences (aka discrete holonomic sequences)
#that satisfy a linear recurrence with POLYNOMIAL coefficients. If its order is k, then for some polynomials
#P_0(n),..., P_k(n), we have
#P_0(n)*a(n)+P_1(n)*a(n-1)+ ... + P_k(n)*a(n-k)=0
#Of course we also need initial conditions, so our data structure for them would be [INI,Rec]
#but now Rec would be written as the triple, [n,[P_0,P_1,..., P_k]]. the first entry ,n, just denotes
#the symbol used to describe the polynomial coefficients of the recurrence. 
#For example the sequence a(n)=n! satisfies (by definition of factorial)
#a(n)-n*a(n-1)=0, a(1)=1, so its representation is going to be
#[[1],[n,[1,-n]]. Also C-finite sequences can be represented in this data structure (now the polynomials P0,P1,..., Pk
#have degree 0, so the representation of the Fibonaii sequence is
#[[1,1],[n,[1,-1,-1]] ]  or [[1,1],[n, [-1,1,1]] (why?)
#The representation of the sequence a(n)=2^n is
#[[2],[n,[1,-2]] ], etc.
#Note that our convention is that INI starts at n=1, rather than n=0, since it is more convenient, because in Maple
#L[i] is the i-th member of the list

#Our first goal is to write a procedure SeqToP(L,n), that inputs a sequence (usually of numbers, but of course, they may
#be symbolic), and a discrete variable name, n, and guesses a P-recursive description. Of course, we must start
#with an auxiliary procedure SeqToP1(L,n,k,d) that looks for such a description with order k and degree of the
#k+1 polynomials P0(n), ..., Pk(n)  equal to d. This requires (k+1)*(d+1) degrees of freedom, and to be safe
#against `over-fitting' let's require 4 extra data points,
#and since we can only start collecting data for n>=k, we need the `training data', the input list L to be at least
#of length (k+1)*(d+1)+k+4. If it is less, we would get an ERROR message.

#SeqToP1(L,n,k,d): Inputs a list of length >= (k+1)*(d+1)+k+4, outputs a P-recursive representation (if it exists) that
#fits the list. For example, try:
#SeqToP1([seq(i!,i=1..20)],n,1,1);
SeqToP1:=proc(L,n,k,d) local eq,var,P,a,i,j,n1,Rec,c:
#Some checking for the input

if not (type(L,list) and type(n,symbol) and type(k,integer) and type(d,integer) and k>=1 and d>=0) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

if nops(L)< (k+1)*(d+1)+k+4 then
 print(`List too short`):
 RETURN(FAIL):
fi:

#end checking the input

#Setting up the k+1 generic polynomials P[i], i=0..,k, each of degree d using a[i,j] for the coefficient of n^j in P[i]
for i from 0 to k do
 P[i]:=add(a[i,j]*n^j,j=0..d):
od:

#Since the set of unknowns (undetermined coefficients) is {a[i,j], i=0..k,j=0..d)} and we are going use linear algebra (the solve command)
#we need to define it
var:={seq(seq(a[i,j],j=0..d),i=0..k)}:

#We now need to set up the set of equations, let's call it eq, by plugging-in into the tentative recurrence the values of the sequence
#we use n1 to access the terms of the sequence
eq:={seq(add(subs(n=n1,P[i])*L[n1-i],i=0..k),n1=k+1..nops(L))}:


#we now solve the system, renaming var as the solution
var:=solve(eq,var):

#We now plug-in into the generting Rec
Rec:=[seq(subs(var,P[i]),i=0..k)]:


#Since the 0 recurrence is always a recurrence, we have to check that it is not 0

if Rec=[0$(k+1)] then
 RETURN(FAIL):
fi:

#since our system is homogeneous, we will get an annoying aribitrary factor, so let's divide by it
c:=lcoeff(Rec[1],n):
Rec:=[seq(normal(Rec[i]/c),i=1..nops(Rec))]:

#We now return the sequence [INI,[n,Rec]] in our data structure
[[op(1..k,L)],[n,Rec]]:

end:

#Now we are ready to keep trying SeqToP1(L,n,k,d) for the `smallest' k and d until it is no longer true that
#nops(L)>=(k+1)*(d+1)+k+4 and then we have to admit failure. Let's define the "complexity" of the P-recurrence
#sequence as k+d, but let's prefer low order to lower degree.


#SeqToP(L,n,MaxC): Inputs a list outputs a P-recursive representation (if it exists), and maximum complexity  that
#fits the list. The output is the lowest Degree+Order that fits it, and we prefer lower order. For example, try:
#SeqToP([seq(i!,i=1..20)],n,2);
SeqToP:=proc(L,n,MaxC) local k,d,C,tem,L1,Rec,i,n1:


for C from 1 to MaxC do
 for k from 1 to C do
  for d from 0 to C-k while nops(L)>=(k+1)*(d+1)+k+4 do
   #For the sake of efficiency, we truncate the list to shorter one
   L1:=[op(1..(k+1)*(d+1)+k+4,L)]:
  
   tem:=SeqToP1(L1,n,k,d):
   if tem<>FAIL then
     #To make sure, we check that the guessed recurrence holds for the full sequence
       Rec:=tem[2][2]:
     if normal({seq(add(subs(n=n1,Rec[i+1])*L[n1-i],i=0..k),n1=k+1..nops(L))})<>{0} then
      print(tem, `did not work out`):
       RETURN(FAIL):
    else
     RETURN(tem):
    fi:
  fi:
  od:
 od:
od:

#If nothing worked, we return FAIL
FAIL:

end:

##end from C12.txt
  



#PtoSeq(P,N): inputs a P-recursive sequence in the format [INI,[n,[P0,.., Pk]] where INI is the list of initial values for n=1, ...,k and
#P0,..., Pk are polynomials in n, representing the coefficients of the linear recurrence with positive coefficients
#P0(n)*x(n)+ ...+ Pk(n)*x(n-k)=0 with those initial conditions, and outputs the first N terms of the sequence
#(stolen from Victoria Chayes' hw12.txt)
PtoSeq:=proc(P,N) local L,n, nt, i, l,k,P0:
L:=P[1]:
n:=P[2][1]:
l:=nops(P[2][2])-1:
k:=nops(L)+1:
P0:=-1/P[2][2][1]:

while k<=N do
  nt:=subs(n=k,P0)*add(subs(n=k, P[2][2][i+1])*L[-i],i=1..l):
  L:=[op(L),nt]:
  k:=k+1:
od:

L:
end:

#EstAsy(L,n): inputs a long list of pos. numbers L and a variable n, estimates
#C,mu,theta, such that L[n] is asympt. to C*mu^n*n^theta
#log(L[n])= logC+n*logmu+theta*log(n). Outputs [mu,theta] 
EstAsy:=proc(L,n) local logC,logmu,theta,eq,var,n1,JS,mu,C:
n1:=nops(L):
var:={logmu,logC,theta}:
JS:=logC+n*logmu+theta*log(n):
eq:={seq(evalf(log(L[i]))=subs(n=i,JS),i=n1-2..n1)}:
var:=solve(eq,var):
mu:=exp(subs(var,logmu)):
theta:=subs(var,theta):
C:=exp(subs(var,logC)):

[C,mu,theta]:
#C*mu^n*n^theta:

end:

 #P0(n)*x(n)+P1(n)*x(n-1)+ ...+Pk(n)*x(n-k)=0 
 #Asy(rec): inputs a P-recursive sequence rec=[INI, [n, [P0,...,Pk]]]
 #and outputs mu and theta such L[n]=C*mu^n*n^theta for some C
 Asy:=proc(rec) local n,rec0,d,INI,i,rec00,mu,p,mu0,LS,x,theta:
 INI:=rec[1]:
 n:=rec[2][1]:
 rec0:=rec[2][2]:
  d:=max(seq(degree(rec0[i],n),i=1..nops(rec0))):
 rec00:=[seq(coeff(rec0[i],n,d),i=1..nops(rec0))]:
 #add(rec00[i+1]*mu^(n-i),i=0..nops(rec00)-1)=0
 p:=numer(add(rec00[i+1]/mu^i,i=0..nops(rec00)-1)):
 mu0:=max(solve(p,mu)):
  if not evalf(mu0)> 0 then
       RETURN(FAIL):
   fi:
 #y(n)=x(n)/mu0^n:
 #rec0:=[rec0[1],rec0[1]/mu, ..., rec0[k]/mu^k]
 rec0:=[seq(rec0[i+1]/mu0^i,i=0..nops(rec0)-1)]:
 #x(n)=n^theta:
 LS:=add(rec0[i+1]*(1-i/n)^theta/n^d,i=0..nops(rec0)-1):
 LS:=normal(subs(n=1/x,LS)):
 LS:=taylor(LS,x=0,3);
  if coeff(LS,x,0)<>0 then
     RETURN(FAIL):
   fi:
 LS:=coeff(LS,x,1):
 [mu0,solve(LS,theta)]:
 end:

  
 #BetterAsy(rec,n,K): inputs a P-recursive sequence rec=[INI, [n, [P0,...,Pk]]]
 #and outputs the asympotic expansion of thesequence to order K
  #L[n]=C*mu^n*n^theta*(1+a1/n+a2/n^2+...+aK/n^K)
# BetterAsy:=proc(rec,n,K) local n,rec0,d,INI,i,rec00,mu,p,mu0,LS,x:
 #to be hw14.txt
# end: