Help := proc(); 
print(`agen(n,k,i,j), alingen(n,k,i,j,coeff)`);
print(`guessLinN(List,N,a), GuessLin(List,a), GuessLinShifts(L,S,a)`);
print(`CheckQuadLin(k,i,j,a), GuessPolyXYn(X,Y,n,x,a), GuessPolyXY(X,Y,x,a)`);
print(`Programs adapted from Professor Zeilberger's EMILIE.txt package:`);
print(`KevenIoneProof(K,n), GlobalPaul(K,n), Paul(K,n)`);
print(`GuessPol(L,n,s0), GuessPWP(L,n), GuessPol1(L,d,n,s0), GrabPolHead(L,n)`);
print(`EmilieRecG1(K,z), EmilieRec(K,N), findCrec(f,N), findrec1(f,DEGREE,ORDER,n,N)`);
end:

#--------------------------------------------------------------
#agen(n,k,i,j): inputs integers n,k,i,j and outputs the nth term
#of the sequence generated by the quadratic recurrence relation
#a(n)*a(n-k)=a(n-i)*a(n-k+i)+a(n-j)+a(n-j+i).


agen := proc(n,k,i,j) option remember;
if n <= k then 
	return 1;
else
	return simplify((agen(n-i,k,i,j)*agen(n-k+i,k,i,j)+agen(n-j,k,i,j)+agen(n-k+j,k,i,j))/agen(n-k,k,i,j));
fi;
end:


#--------------------------------------------------------------
#alingen(n,k,i,j,coeff): inputs integers n,k,i,j,coeff and outputs
#the nth term of the sequence generated by the linear recurrence
#relation a(n)=coeff*a(n-2*i*j)-coeff*a(n-4*i*j)+a(n-6*i*j)

alingen := proc(n,k,i,j,coeff) option remember;

if n<6*i*j+1 then 
	RETURN(agen(n,k,i,j));
else
	RETURN( coeff*(alingen(n-2*i*j,k,i,j,coeff)-alingen(n-4*i*j,k,i,j,coeff))+alingen(n-6*i*j,k,i,j,coeff));
fi;

end:


#--------------------------------------------------------------
#guessLinN(List,N,a): Inputs a sequence, List, an integer N, and
#a symbol a. Outputs a linear recurrence of order N that 
#annihilates List.

guessLinN := proc(List,N,a) local eq,var,c,i,newEq,solutions,r;

if N > nops(List) then
	RETURN(`List not long enough`);
fi;

eq := {};
var := {};
c := 0;
#print(N*floor(evalf(nops(List)/N)));
for i from 1 to nops(List)-N-1 do

	newEq := sum(List[j]*a[j-i], j=i..N+i-1);
	#print(i,newEq);
	c := c+N;
	eq := eq union {newEq=0};
	var := var union {seq(a[k-i], k=i..N+i-1)};
	
od;

#print(eq);

solutions := solve(eq,var);
if solutions <> null then
	r := createRecurrence(solutions,a);
	RETURN(r);
else 
	RETURN(FAIL);
fi;
end:

#--------------------------------------------------------------

createRecurrence := proc(List,a) local s, out;

s := sum(a[j]*f(n+j-nops(List)+1), j=0..nops(List)-1);

out := factor(subs(List,s));

out;

end:

#--------------------------------------------------------------
#GuessLin(List,a): intputs a sequence, List, and a symbol a,
# and outputs a linear recurrence of minimal order that
#annihilates List.

GuessLin := proc(List,a) local i,s;

for i from 1 to floor(nops(List)/2) do

	s := guessLinN(List,i,a);
	if s <> FAIL and s <> 0 then
		RETURN(s);
	fi;

od;

RETURN(FAIL);

end:

#--------------------------------------------------------------
#GuessLinShifts(L,S,a): inputs a sequence, L, a list of shifts, S,
#and a symbol a. Outputs a recurrence [n1,n2,...,n(nops(S))] such that
#add(nj*L[i-s],s in S)=0. 

GuessLinShifts := proc(L,S,a) local n,eq,var,ans;

n := nops(S);
eq := {seq(L[i]=sum(a[j]*L[i-S[j]],j=1..n),i=S[n]+1..nops(L))};
var:=[seq(a[k], k=1..n)];

ans := solve(eq,{op(var)});
if ans = NULL then
	RETURN(FAIL);
fi;

var := subs(ans,var);
[var,S];

end:

#--------------------------------------------------------------
#CheckQuadLin(k,i,j,a): inputs integers k,i,j and a symbol a. 
#the program checks to see if the sequence seq(a(n,k,i,j),n=1..4*6*i*j)
#satisfies the linear recurrence guessed by GuessLinShifts(L,S,a)
#where S=[2*i*j,4*i*j,6*i*j]. returns false if it doesn't work
#out, and the recurrence if it does work.

CheckQuadLin := proc(k,i,j,a) local S,L,Lin,test;

S := [2*i*j,4*i*j,6*i*j];
L := [seq(agen(n,k,i,j),n=1..4*6*i*j)];
Lin := GuessLinShifts(L,S,a);
if Lin <> FAIL then
test:={seq(L[m]-Lin[1][1]*L[m-S[1]]-Lin[1][2]*L[m-S[2]]-Lin[1][3]*L[m-S[3]],m=S[3]+1..nops(L))};
	if test = {0} then
		RETURN(Lin);
	else
		RETURN(Lin,false);
	fi;
fi;
FAIL;

end:

#--------------------------------------------------------------

GuessPolyXYn := proc(X,Y,n,x,a) local P,var,m,eq,var1;

P := a[n+1]^x*add(a[i]*x^i, i=0..n);
var := {seq(a[i], i=0..n+1)};
m := min(nops(X),nops(Y));
eq := {seq(subs(x=X[j],P)=Y[j],j=1..m)};

var1 := {solve(eq,var)};

if nops(var1)=0 then
	RETURN(FAIL);
fi;

P := subs(op(var1),P);

RETURN(P);

end:

#--------------------------------------------------------------

GuessPolyXY := proc(X,Y,x,a) local i,P;

for i from 1 do

	P := GuessPolyXYn(X,Y,i,x,a);
	if P <> FAIL then
		RETURN(P);
	fi;

od;

end:

#--------------------------------------------------------------

#KevenIoneProof(K,n): completes the proof of the linear recurrence
#formula for the case when K is even and i=1.

KevenIoneProof:=proc(K,n) local lu,gu,i,j,z:
print(`Theorem:`):
print(` Let K be a positive integer`):
print(`Define an integer sequence a(n) (n>=0) by the non-linear recurrence`):
print(``):
print(a(n)=(a(n-1)*a(n-2*K+1)+2*a(n-K))/a(n-2*K)):
print(`subject to the initial conditions a(i)=1 for 0<=i<=2K`):
print(`Then lo-and-behold all the terms are integers `):
print(``):
lu:=GlobalPaul(K,n):
gu:=EmilieRecG1(K,z):

print(`Let b(n) be the sequence annihilated by the linear`):
print(`recurrence equation with constant coefficients`):
print(add(coeff(gu,z,i)*b(n-K*2*i+i),i=0..degree(gu,z))=0 ):
print(`With the following initial conditions`):
for j from  1 to nops(lu) do
 print(`For n in the discrete interval`):
 print(lu[j][1] ):
 print(`it equals the polynomial`):
 print(lu[j][2]):
od:
print(`Note that b(n) is well-defined and manifestly an integer`):
print(` (by induction)`):
print(` I claim that b(n)=a(n).`):
print(`All we have to prove is that b(n) satisfies the same non-linear`):
print(`recurrence as a(n), with the same initial conditions`):
print(`The initial conditions are right, of course, Emilie`):
print(`has already proved that if `):
print(`b(n)b(n-2K)-b(n-1)b(n-2K+1)-2b(n-K)=0`):
print(`is true for n<=6K-2, then it is true for ever after.`):
print(`So it only remains to prove this for n<=6K-2`):
print(`But this is a routine verification`):
print(`using the above explicit expressions for b(n), and is indeed true.`):
print(`  QED`):
end:

#--------------------------------------------------------------

#GlobalPaul(K,n): Given symbolic K and n conjecutures
#the first 6K+3 terms of the Hogan sequence
#Try: GlobalPaul(K,n);

GlobalPaul:=proc(K,n) local K1,GU,LU,a,b,pol,i,j,LU1:
GU:=[]:

for K1 from 10 to 20 do
 GU:=[op(GU),[op(1..5,Paul(K1,n))]]:
od:

LU:=[]:

for i from 1 to 5 do
a:=GuessPol( [seq(GU[j][i][1][1],j=1..nops(GU))],K,1):
b:=GuessPol( [seq(GU[j][i][1][2],j=1..nops(GU))],K,1):
pol:=GuessPol( [seq(GU[j][i][2],j=1..nops(GU))],K,1):
if a=FAIL or b=FAIL or pol=FAIL then
 RETURN(FAIL):
fi:
a:=expand(subs(K=K-9,a)):
b:=expand(subs(K=K-9,b)):
pol:=expand(subs(K=K-9,pol)):

LU:=[op(LU),[[a,b],pol]]:
od:

for K1 from 1 to 20 do
LU1:=subs(K=K1,LU):
if
[seq(
seq(subs(n=i,LU1[j][2]),i=LU1[j][1][1]..LU1[j][1][2])
,j=1..5)]
<>[seq(anK(i,K1),i=1..6*K1-2)] then
 print(`It didn't work out at K=`, K1):
  RETURN(LU):
fi:
od:
LU:
end:

#--------------------------------------------------------------

#Paul(K,n): guesses a piece-wise polynomial beginning
#to the Hogan sequence for k=2K, for example,
#try: Paul(5,n);
Paul:=proc(K,n) local L, n1:
L:=[seq(anK(n1,K),n1=1..7*K+20)]:
GuessPWP(L,n):
end:

#--------------------------------------------------------------

#GuessPol(L,n,s0): guesses a polynomial of degree d in n for
# the list L, such that L[i]=P[i] for i>=s0 for example, try: 
#GuessPol([(seq(i,i=1..10),n,1);

GuessPol:=proc(L,n,s0) local d,gu:

for d from 0 to nops(L)-s0-3 do
 gu:=GuessPol1(L,d,n,s0):
 if gu<>FAIL then
    RETURN(gu):
 fi:
od:

FAIL:

end:

#--------------------------------------------------------------

anK:=proc(n,K) local k:option remember:
k:=2*K:

if n<=k then
 1:

else
(anK(n-1,K)*anK(n-k+1,K)+2*anK(n-k/2,K))/anK(n-k,K):
fi:
end:

#--------------------------------------------------------------

#GuessPWP(L,n): guesses a piece-wise polynomial for a sequence L
#as far as it goes
#starting at L[1], it outputs a list whose entries are
#[[a,b],pol] which means that L[i]=pol(i) for a<=i<=b
#For example, try:
#GuessPWP([1$10,2$10],n)
GuessPWP:=proc(L,n) local L1,pol,gu,mu,kodem:
gu:=[]:
L1:=L:
mu:=GrabPolHead(L1,n):
if mu=FAIL then
 RETURN(FAIL):
fi:
kodem:=0:
while mu<>FAIL do
pol:=mu[1]:
gu:=[
      op(gu),
      [ [kodem+1,kodem+mu[2]], expand(subs(n=n-kodem,pol)) ] 
    ]:
kodem:=kodem+mu[2]:
L1:=[op(kodem+1..nops(L),L)]:
mu:=GrabPolHead(L1,n):
od:
gu:
end:

#--------------------------------------------------------------

#GuessPol1(L,d,n,s0): guesses a polynomial of degree d in n for
# the list L, such that L[i]=P[i] for i>=s0 for example, try: 
#GuessPol1([(seq(i,i=1..10),1,n,1);
GuessPol1:=proc(L,d,n,s0) local P,i,a,eq,var:
if d>=nops(L)-s0-2 then
 ERROR(`the list is too small`):
fi:

P:=add(a[i]*n^i,i=0..d):
var:={seq(a[i],i=0..d)}:
eq:={seq(subs(n=i,P)-L[i],i=s0..s0+d+3)}:

var:=solve(eq,var):

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

subs(var,P):

end:

#--------------------------------------------------------------

#GrabPolHead(L,n): grabs the beginning of the list L
#that is a polynomial (if found) together with the number of
#terms it agrees with
GrabPolHead:=proc(L,n) local L1,d,pol,i:
pol:=FAIL:
for d  from 0 to nops(L)-4 while pol=FAIL do
L1:=[op(1..d+4,L)]:
pol:=GuessPol1(L1,d,n,1):
od:
if pol=FAIL then
 RETURN(FAIL):
fi:

for i from 1 to nops(L) do
 if  L[i]<>subs(n=i,pol)  then
   RETURN(pol,i-1):
 fi:
od:
pol,nops(L):
end:

#--------------------------------------------------------------

#EmilieRecG1(K,z) : guesses a linear recurrence equation with
#constant coefficients for the Hogan sequence for symbolic K
#Try: EmilieRecG(K,z);

EmilieRecG1:=proc(K,z) local K1,lu,i,j,N:
lu:=[seq(EmilieRec(K1,N),K1=1..6)]:

add(GuessPol([seq(coeff(lu[i],N,2*i*j-j),i=1..nops(lu))],K,1)*z^j,j=0..3):
end:

#--------------------------------------------------------------

#EmilieRec(K,N) : guesses a linear recurrence equation with
#constant coefficients for the Hogan sequence for K
#Try: EmilieRec(1,N);
EmilieRec:=proc(K,N) local gu,n1:
gu:=[seq(anK(n1,K),n1=1..12*K+30)]:
findCrec(gu,N):
end:

#--------------------------------------------------------------

#findCrec(f,N): Finds a linear recurrence equation with
#constant coefficients for the sequence f
#For example try: findCrec([1$20],N);
findCrec:=proc(f,N) local d,gu,n:
gu:=FAIL:

for d from 0 to nops(f)-6 do
 gu:=findrec1(f,0,d,n,N):
 if gu<>FAIL then
   RETURN(gu):
 fi:
od:
FAIL:
end:

#--------------------------------------------------------------

#findrec1(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrec([seq(i,i=1..10)],0,2,n,N);
findrec1:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a:
option remember:

if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to nops(f)-ORDER do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
  od:
 
   eq:= eq union {eq1}:
od:
 
var1:=solve(eq,var):
 
kv:={}:
 
for i from 1 to nops(var1) do
 mu:=op(i,var1):
 
if op(1,mu)=op(2,mu) then
   kv:= kv union {op(1,mu)}:
 fi:
 
od:

ope:=subs(var1,ope):

if ope=0 then
  RETURN(FAIL):
fi:

ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}:

if nops(ope)>1 then
RETURN(Yafe(ope[1],N)[2]):
elif nops(ope)=1 then
RETURN(Yafe(ope[1],N)[2]):
else
 RETURN(FAIL):
fi:

end:

#--------------------------------------------------------------

Yafe:=proc(ope,N) local i,ope1,coe1,L: 
if ope=0 then
 RETURN(1,0):
fi:
ope1:=expand(ope):
L:=degree(ope1,N):
coe1:=coeff(ope1,N,L):
ope1:=normal(ope1/coe1):
ope1:=normal(ope1):
ope1:=
convert(
[seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`):
factor(coe1),ope1:
end: