######################################################################
## EMILIE: Save this file as  EMILIE to use it,                      #
# stay in the                                                        #
## same directory, get into Maple (by typing: maple <Enter> )        #
## and then type:  read EMILIE: <Enter>                               #
## Then follow the instructions given there                          #
##                                                                   #
## Written by Doron Zeilberger, Rutgers University ,                 #
##  zeilberg at math dot rutgers dot edu.                            # 
######################################################################


print(` Started : March 28, 2006: tested for Maple 10 `):
print(`Version of Sept. 19, 2006:  `):
print():
print(`This is EMILIE, A Maple package that completes the proof`):
print(`of the main theorem in the beautiful on-line paper`):
print(`"New Family of Somos-like Recurrences" by Paul Heideman and Emilie `):
print(` Hogan`):

print(`www.math.wisc.edu/~propp/SSL/heiho.pdf `):
print(`It also conjectures [ EmilieS(K,x,Z);] a general generating function`):

print():
print(`The most current version is available on WWW at:`):
print(` http://www.math.rutgers.edu/~zeilberg/tokhniot/EMILIE .`):
print(`Please report all bugs to: zeilberg at math dot rutgers dot edu .`):
print():
print(`For general help, and a list of the functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name);" `):
print():
 

ezra1:=proc()
if args=NULL then

print(` The supporting and less important procedures are `):
print(`CheckEmilie, Emilie, EmilieS , GrabPolHead `):
print(` findCrec, GuessPol, GuessPWP, GuessR, GuessR1, Paul`):
else
ezra(args):
fi:
end:

ezra:=proc()
if args=NULL then

print(` EMILIE: A Maple package for  `):

print(`The MAIN procedures are: `):
print(` anK, EmiliePaul, EmilieRec, EmilieRecG, GlobalPaul`):

elif nargs=1 and args[1]=anK then
print(` anK(n,K) : The sequence a_n of the Heideman-Hogan paper`):
print(`where k=2*K+1` ):
print(`For example, try: anK(4,3); `):

elif nargs=1 and args[1]=CheckEmilie then
print(` CheckEmilie(K,L) : Checks that the coefficients of the conjectured generating function`):
print(`for Emilie(K,x) satisfy the quadratic recurrence (1) of the paper up to n=L`):
print(`For example, try: CheckEmilie(5,1000);`);

elif nargs=1 and args[1]=Emilie then
print(`Emilie(K,x): Conjectures a generating function, in x, for the Somos-Like sequence`):
print(`satisfying the  recurrence`):
print(` a(n)a(n-2*K-1)=a(n-1)*(a(n-2*K)+a(n-K)+a(n-K-1)`):
print(`subject to the inital conditions a(m)=1 for m=0...2*K`):
print(`Note that in the Heideman-Hogan paper, k=2*K+1`):
print(`For example, try: Emilie(5,x);`):

elif nargs=1 and args[1]=EmiliePaul then
print(`EmiliePaul(K,n): completes the proof of Paul Heideman and`):
print(`Emilie Hogan's beautiful conjecture`):

elif nargs=1 and args[1]=EmilieRec then
print(`EmilieRec(K,N) : guesses a linear recurrence equation with`):
print(`constant coefficients for the Heideman-Hogan sequence for K`):
print(`Try: EmilieRec(1,N);`):

elif nargs=1 and args[1]=EmilieRecG then
print(`EmilieRecG(K,N) : guesses a linear recurrence equation with`):
print(`constant coefficients for the Heideman-Hogan sequence for SYMBOLIC K`):
print(`Try: EmilieRecG(K,N);`):

elif nargs=1 and args[1]=EmilieS then
print(`EmilieS(K,x,Z): The conjectured expression for the Emilie(K,x) in terms of`):
print(` Z=x^K. Here K,x,Z are symbols`):
print(`Note that subs({K=K0,Z=x^K0},EmilieS(K,x,Z)); should be Emilie(K0,x); `):
print(`for any numeric K0`):
print(`For example, try: EmilieS(K,x,Z);`);

elif nargs=1 and args[1]=findCrec then
print(`findCrec(f,N): Finds a linear recurrence equation with`):
print(`constant coefficients for the sequence f`):
print(`For example try: findCrec([1$20],N);`):

elif nargs=1 and args[1]=GlobalPaul then
print(`GlobalPaul(K,n): Given symbolic K and n conjecutures`):
print(`the first 6K terms of the Heideman-Hogan sequence`):
print(`Try: GlobalPaul(K,n);`):

elif nargs=1 and args[1]=GrabPolHead then
print(`GrabPolHead(L,n): grabs the beginning of the list L`):
print(`that is a polynomial (if found) together with the number of`):
print(`terms it agrees with`):
print(`For example, try:`):
print(`GrabPolHead([1$8,seq(i^3,i=1..13)],n);`):

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

elif nargs=1 and args[1]=GuessPWP then
print(`GuessPWP(L,n): guesses a piece-wise polynomial for a sequenec L`):
print(`as far as it goes`):
print(`starting at L[1], it outputs a list whose entries are`):
print(`[[a,b],pol] which means that L[i]=pol(i) for a<=i<=b`):
print(`For example, try:`):
print(`GuessPWP([1$10,2$10],n); `):

elif nargs=1 and args[1]=GuessR then
print(`GuessR(L,x): Given a list of numbers L`):
print(`guesses a rational function whose first`):
print(`nops(L) terms are given by L`):
print(`for example, try: GuessR([1$30],x);`):
print(`and you should get 1/(1-x) `):
print(`since the Maclaurin expansion of 1/(1-x) starts with 1+x+x^2+x^3+ ...`):

elif nargs=1 and args[1]=GuessR1 then
print(`GuessR1(L,p,q,x): Given a list of numbers L`):
print(`guesses a rational function in x, with numerator degree p`):
print(`and  denominator degree q`):
print(`whose first`):
print(`nops(L) terms, in its Maclaurin expansion, are given by `):
print(`the entries of L.`):
print(`For example, try: GuessR1([1$30],0,1,x);`):
print(`and you should get 1/(1-x) `):
print(`since the Maclaurin expansion of 1/(1-x) starts with 1+x+x^2+x^3+ ...`):

elif nargs=1 and args[1]=Paul then
print(`Paul(K,n): guesses a piece-wise polynomial (in n)  beginning`):
print(`to the Heideman-Hogan sequence for k=2K+1, for example,`):
print(`try: Paul(5,n);`):


else
 print(`There is no such thing as`, args):

fi:


end:

ez:=proc(): print(`anK(n,K), an1(n), GuessR1(L,p,q,x), GuessR(L,x)`):
print(`Emilie(K,x), Mekh(x,K) `):
print(`A1(K,x), B1(K,x), C1(K,x), D1(K,x), E1(K,x), CheckE(K,L)  `):
print(`EmilieS(K,x,Z), Emilie1(K,x,Z, EmilieS2(K,x,Z) `):
end:

an1:=proc(n) option remember:
if n<=3 then
 1:

else
(an1(n-1)*an1(n-2)+an1(n-1)+ an1(n-2))/an1(n-3):
fi:
end:


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

if n<=k then
 1:

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

GuessR1:=proc(L,p,q,x) local a,b,T,B,i,var,gu,eq:
T:=add(a[i]*x^i,i=0..p):
B:=add(b[i]*x^i,i=0..q):
var:={seq(a[i],i=0..p),seq(b[i],i=0..q)}:
gu:=expand(add(L[i+1]*x^i,i=0..p+q+11)*B-T):
eq:={seq(coeff(gu,x,i),i=0..p+q+5)}:
var:=solve(eq,var):


if var=NULL or subs(var,B)=0 then
 RETURN(FAIL):
else
gu:=normal(subs(var,T)/subs(var,B)):
RETURN(sort(numer(gu)) /sort(denom(gu)) ):
fi:
end:

#GuessR(L,x): Given a list of numbers L
#guesses a rational function whose first
#nops(L) terms are given by L
#for example, try: GuessR([1$30],x);
GuessR:=proc(L,x) local p,q,S1:
for S1 from 0 to nops(L)-14 do
for p from 0 to S1 do
for q from 0 to S1-p do
 if GuessR1(L,p,q,x)<>FAIL then
    RETURN(GuessR1(L,p,q,x)):
 fi:
od:
od:
od:
FAIL:
end:

Emilie:=proc(K,x) local L, k,n:
k:=2*K+1:
L:=[seq(anK(n,K),n=1..6*k+10)]:
GuessR1(L,3*k-4,3*k-3,x);
end:


w:=proc(K): 2*K^2+8*K+3:end:

Mekh:=proc(K,x): (1-x^(2*K))*(1-w(K)*x^(2*K)+x^(4*K)):end:



A1:=proc(K,x): normal((x^(2*K)-1)/(x-1)):end:

B1:=proc(K,x) local i:
expand(x^(2*K)*add((w(K)-2*i)*x^i,i=0..K)):
end:

C1:=proc(K,x) local i:
expand(x^(3*K+1)*add((w(K)-2*K-2*(i+1)*(i+2))*x^i,i=0..K-2)):
end:

D1:=proc(K,x) local i:
expand(x^(5*K-1)*add((2*K+1+2*(i+1)*(i+2))/x^i,i=0..K-1)):end:

E1:=proc(K,x) local i:
expand(x^(6*K-1)*add((2*i+3)/x^i,i=0..K-1)):end:


EmilieC:=proc(K,x) : sort(A1(K,x)-B1(K,x)-C1(K,x)+D1(K,x)+E1(K,x))/Mekh(K,x):end:

CheckEmilie:=proc(K,L) local gu,x,k,i:
k:=2*K+1:
gu:=Emilie(K,x) : 
gu:=taylor(gu,x=0,L+1):
gu:=[seq(coeff(gu,x,i),i=0..L)]:
gu:=[seq(gu[i]*gu[i-k]-gu[i-1]*gu[i-k+1]-gu[i-K]-gu[i-K-1],i=k+1..L)]:
evalb(gu=[0$(L-k)]):
end:


As:=proc(K,x): normal((x^(2*K)-1)/(x-1)):end:

Bs:=proc(K,x) local i:
normal(sum((w(K)-2*i)*x^(i+2*K),i=0..K)):
end:

Cs:=proc(K,x) local i:
normal(sum((w(K)-2*K-2*(i+1)*(i+2))*x^(3*K+1+i),i=0..K-2)):
end:

Ds:=proc(K,x) local i:
normal(sum((2*K+1+2*(K-i)*(K-i+1))*x^(4*K+i),i=0..K-1)):end:

Es:=proc(K,x) local i:
normal(sum((2*(K-i)+1)*x^(5*K+i),i=0..K-1)):end:

MoneS:=proc(K,x) : normal(As(K,x)-Bs(K,x)-Cs(K,x)+Ds(K,x)+Es(K,x)):end:





EmilieS:=proc(K,x,Z) local T,B: 
T:=
-(1-2*x+4*Z^4-4*Z^2-Z^6+8*Z^4*x^2-24*Z^4*x*K-4*Z^4*x*K^2+2*Z^4*x^2*K^2-8*Z^2*x^2*K+16*Z^4*x^2*K+16*Z^2*x*K+4*Z^2*x*K^2-2*Z^
2*x^2*K^2+x^2-2*Z^2*K^2+2*Z^3*x^2+10*Z^2*x+2*Z^4*K^2+2*Z^3*x-2*Z^5*x^2-6*Z^2*x^2+8*Z^4*K-8*Z^2*K-3*Z^6*x^2-12*Z^4*x-2*Z^5*x
+4*Z^6*x)/(x-1)^3:

B:=(1-Z^2)*(1-w(K)*Z^2+Z^4):
T/B:
end:

EmilieS1:=proc(K,x,Z) local luE,luO,T,B: 
luE:=
-1+2*x-8*Z^4*x^2+2*Z^2*K^2-10*Z^2*x-2*Z^4*K^2+6*Z^2*x^2-8*Z^4*K+8*Z^2*K+3*Z^6*x^2+12*Z^4*x-4*Z^6*x+4*Z^2+Z^6-x^2-4*Z^4+24*Z
^4*x*K+4*Z^4*x*K^2-2*Z^4*x^2*K^2+8*Z^2*x^2*K-16*Z^4*x^2*K-16*Z^2*x*K-4*Z^2*x*K^2+2*Z^2*x^2*K^2:
luO:=-2*Z^2*x^2-2*Z^2*x+2*Z^4*x^2+2*Z^4*x:
T:=-expand(luE+luO*Z):
B:=(1-Z^2)*(1-w(K)*Z^2+Z^4)*(1-x)^3:
T/B:
end:


EmilieS2:=proc(K,x,Z) local luE,luO,B: 
luE:=
-1+2*x-8*Z^4*x^2+2*Z^2*K^2-10*Z^2*x-2*Z^4*K^2+6*Z^2*x^2-8*Z^4*K+8*Z^2*K+3*Z^6*x^2+12*Z^4*x-4*Z^6*x+4*Z^2+Z^6-x^2-4*Z^4+24*Z
^4*x*K+4*Z^4*x*K^2-2*Z^4*x^2*K^2+8*Z^2*x^2*K-16*Z^4*x^2*K-16*Z^2*x*K-4*Z^2*x*K^2+2*Z^2*x^2*K^2:
luO:=-2*Z^2*x^2-2*Z^2*x+2*Z^4*x^2+2*Z^4*x:
luE:=subs(Z=Z^(1/2),luE);
luO:=subs(Z=Z^(1/2),luO);
B:=(1-Z)*(1-w(K)*Z+Z^2)*(1-x)^3:
luE/B,luO/B:
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:


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



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


#GuessPWP(L,n): guesses a piece-wise polynomial for a sequenec 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:



#Paul(K,n): guesses a piece-wise polynomial beginning
#to the Hindeman-Hogan sequence for k=2K+1, 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:

#GlobalPaul(K,n): Given symbolic K and n conjecutures
#the first 6K terms of the Heideman-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+1)] then
 print(`It didn't work out at K=`, K1):
  RETURN(LU):
fi:
od:
LU:
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:
 
#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:
 

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

anKS:=proc(K,L) local n1:[seq(anK(n1,K),n1=1..L)]:end:

#EmilieRec(K,N) : guesses a linear recurrence equation with
#constant coefficients for the Heideman-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:



#EmilieRecG(K,N) : guesses a linear recurrence equation with
#constant coefficients for the Heideman-Hogan sequence for symbolic K
#Try: EmilieRecG(K,N);
EmilieRecG:=proc(K,N) local K1,lu,i,j,z:
lu:=[seq(EmilieRec(K1,N),K1=1..6)]:
add(GuessPol([seq(coeff(lu[i],N,2*i*j),i=1..nops(lu))],K,1)*N^(2*K*j),j=0..3):
end:

#EmilieRecG1(K,z) : guesses a linear recurrence equation with
#constant coefficients for the Heideman-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),i=1..nops(lu))],K,1)*z^j,j=0..3):
end:

#EmiliePaul(K,n): completes the proof of Paul Heideman and Emilie Hogan's
#beautiful conjecture
EmiliePaul:=proc(K,n) local lu,gu,i,j,z:
print(`Theorem:`):
print(`(Cojectured and alomost proved by Paul Heideman and Emilie Hogan,`):
print(`in undergraduate research advised by Jim Propp, at the U. of Wisc.)`):
print(``):
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)+a(n-K)+ a(n-K-1))/a(n-2*K-1)):
print(`subject to the initial conditions a(i)=1 for 0<=i<=2K`):
print(`Then lo-and-behold all the terms are integers `):
print(``):
print(`Completion of the proof by Shalosh B. Ekhad `):
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*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, and Paul and Emilie`):
print(`have already proved that if `):
print(b(n)*b(n-2*K-1)-b(n-1)*b(n-2*K)-b(n-K)-b(n-K-1)=0):
print(`is true for n<=6K, then it is true for ever after.`):
print(`So it only remains to prove this for n<=6K`):
print(`But this is a routine verification`):
print(`using the above explicit expressions for b(n), and is indeed true.`):
print(`  QED`):
end: