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

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



#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: