From: Doron Zeilberger zeilberg AT ekhad.rutgers.edu
Date: Mon, 19 Aug 2013 14:23:31 -0400
Subject: source code (Hardin)

ez:=proc(): print(`Alp(k),  Hkx(k,x),  Gkx(k,x)`): 
print(` IsFollower(w,v), TM(k), Sipur(k,x,L)  `):
end:

#Alp(k): The set of words in  {0,1} of length k
#with no 00
Alp:=proc(k) local gu,mu,mu1:
option remember:
if k=0 then
 RETURN({[]}):
fi:
if k=1 then
 RETURN({[0],[1]}):
fi:

mu:=Alp(k-1):
gu:={}:

for mu1 in mu do
 if mu1[nops(mu1)]=0 then
   gu:=gu union {[op(mu1),1]}:
 else
     gu:=gu union {[op(mu1),0],[op(mu1),1]}:
 fi:
od:
gu:
end:

#IsFollower(w,v): Is v a legal follower of w?
IsFollower:=proc(w,v) local i,k:
k:=nops(w):
if nops(v)<>k then
  RETURN(FAIL):
fi:

 for i from 1 to k do
  if w[i]=0 and v[i]=0 then 
   RETURN(false):
  fi:
od:

 for i from 1 to k-1 do
  if w[i]=0 and v[i+1]=0 then 
   RETURN(false):
  fi:
od:
true:
end:



#TM(k): The transfer matrix of k by n Hardin matrices
TM:=proc(k) local i,j,gu,A,A1,w:
gu:=convert(Alp(k),list):
A:=[]:
for i from 1 to nops(gu) do
w:=gu[i]:
 A1:=[]:
  for j from 1 to nops(gu) do
   if IsFollower(w,gu[j]) then
     A1:=[op(A1),1]:
    else
      A1:=[op(A1),0]:
   fi:
  od:
A:=[op(A),A1]:
od:
A:
end:



#Gkx(k,x): the generating function, in x,
#whose coefficient of x^n is the number of
#all k by n Hardin matrices (without restiction
#that the top-left entry is 1). See A226444.
Gkx:=proc(k,x) local gu,z,eq, eq1,var,i,j:
option remember:
gu:=Alp(k):
var:={seq(z[gu[i]],i=1..nops(gu))}:
eq:={}:

for i from 1 to nops(gu) do
 eq1:=z[gu[i]]-x:
 for j from 1 to nops(gu) do
   if IsFollower(gu[i],gu[j]) then
    eq1:=eq1-x*z[gu[j]]:
   fi:
 od:
eq:=eq union {eq1}:
od:


var:=solve(eq,var):

normal(add(subs(var,z[gu[i]]),i=1..nops(gu))):
end:



#Hkx(k,x): the generating function, in x,
#whose coefficient of x^n is the number of
#all k by n Hardin matrices WITH the restiction
#that the top-left entry is 1. See A228285.
Hkx:=proc(k,x) local gu,z,eq, eq1,var,i,j,mu:
option remember:
gu:=Alp(k):
var:={seq(z[gu[i]],i=1..nops(gu))}:
eq:={}:

for i from 1 to nops(gu) do
 eq1:=z[gu[i]]-x:
 for j from 1 to nops(gu) do
   if IsFollower(gu[i],gu[j]) then
    eq1:=eq1-x*z[gu[j]]:
   fi:
 od:
eq:=eq union {eq1}:
od:


var:=solve(eq,var):

mu:=0:

for i from 1 to nops(gu) do
if gu[i][1]=0 then
 mu:=mu+subs(var,z[gu[i]]):
fi:
od:
normal(mu):
end:


#Sipur(K,x,L): inputs positive integers K and x
#and tells the story of the rigorous enumeration
#of k by n Hardin matrices, for k from 1 to K
#it also prints out 
#the first L terms of each sequence
#and the enumerating sequnce 
#for square matrices for sizes from 1 to K,

#Try: Sipur(6,x);
Sipur:=proc(K,x,L) local f, gu,mu,k,i,t0:
t0:=time():
print(`Rigorously derived generating functions for`):
print(` enumerating k by n Hardin matrices for k from 1`):
print(`to `, K):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Definition: A Hardin matrix has entries from {0,1}`):
print(`with no adjacent 1's horizontally, vertically, or`):
print(`in the NE-SW direction, and whose top-left entry is 1`):
print(``):
print(`Definition:For a fixed k, Let f(k)(x) be the `):
print(`generating function`):
print(`whose coefficient of x^n is the number of n by k`):
print(`We have the following (rigorously-derived) facts`):
print(`followed by the first`, L, `terms of the enumerating`):
print(`sequence (for  the sake of Sloane) `):
mu:=[]:
for k from 1 to K do
gu:=Hkx(k,x):
lprint(f[k](x)=gu):
gu:=taylor(gu,x=0,L+1):
print([seq(coeff(gu,x,i),i=1..L)]):
mu:=[op(mu),coeff(gu,x,k)]:
od:

print(`The first`, K, `terms  of the enumerating sequence`):
print(`for square Hardin matrices is`):
print(mu):

print(`This took `, time()-t0, `seconds, to generate `):

end:
#Doron,
#My good friend Ron Hardin has been enumerating
#matrices of many different classes over the past several years.

#His method is to collect a LOT of data, and then to
#look for patterns in the numbers.

#He just posted the following message to the Sequence
#Fans Mailing List.

#(start)
#The orders of the column (or row) linear recurrences of this s#eem to form
#the Fibonacci series.

#
#T(n,k)=Number of nXk binary arrays with top left value 1 and no two ones
#adjacent horizontally, vertically or nw-se diagonally

#Table starts
#..1...1.....2......3........5.........8..........13.........
#..21.............34

#..1...1.....3......6.......13........28..........60.........
#.129............277

#..2...3....13.....35......112.......337........1034.........
#3154...........9637

#..3...6....35....133......587......2448.......10414........
#44024.........186414

#..5..13...112....587.....3631.....21166......126119.......
#745178........4416695

#..8..28...337...2448....21166....172082.....1428523.....
#11771298.......97268701

#.13..60..1034..10414...126119...1428523....16566199....
#190540884.....2197847780

#.21.129..3154..44024...745178..11771298...190540884...
#3057290265....49208639399

#.34.277..9637.186414..4416695..97268701..2197847780..
#49208639399..1105411581741

#.55.595.29431.789100.26150120.802886174.25325358687.
#791176762937.24801939723742

#Column 1 is A000045
#Column 2 is A002478(n-1)

#Empirical for column k:
# k=1: a(n)=a(n-1)+a(n-2)
#k=2: a(n)=a(n-1)+2*a(n-2)+a(n-3)
#k=3: a(n)=a(n-1)+5*a(n-2)+4*a(n-3)-a(n-5)
#k=4: a(n)=a(n-1)+10*a(n-2)+15*a(n-3)+4*a(n-4)-6*a(n-5)-a(n-6)+3*
#a(n-7)-a(n-8)
#k=5: [order 13]
#k=6: [order 21]
#k=7: [order 34]
#k=8: [order 55]
#k=9: [order 89]

#(end)