

#                      PATTERNS AND FRACTIONS (Maple Version)
#  
#                                      By
# 
#                         Aaron Robertson, Doron Zeilberger
#
#                  Temple University, Philadelphia, PA 19122, USA
#                        [aaron,zeilberg]@math.temple.edu 
#                  http://www.math.temple.edu/~[aaron,zeilberg]/
#
#                                and Herb Wilf
#
#             University of Pennsylvania, Philadelphia, PA, 19104, USA
#                            wilf@math.upenn.edu
#                        http://www.cis.upenn.edu/~wilf/


#    ABSTRACT:   We find the generating function for the number of 
#             	 (132)-avoiding permuatations that have a given
#		 number of (123)-patterns, in the form of a continued
#		 fraction.  We show how to extend this to permutations
#		 that have exactly one (132) pattern.  We find some
#		 properties of the continued fraction, which is similar
#		 to, though more general than, those that were studied
#		 by Ramanujan, and raise some questions about it.


#    In a recent article: `Permutations Containing and Avoiding (123)
# and (132) Patterns', one of us (AR) (see http://www.math.temple.edu/~aaron/)
# obtained expressions for the number of permutations on n objects that have 
# exactly s (132) patterns and r (123) patterns for (0<=r,s<=1). 

#    The number of (132) (resp. (123)) patterns of a permutation
# pi of length n, #132(pi) (resp. #123(pi)), is the number of triples 
# 1<=i1<i2<i3<=n such that pi[i1]<pi[i3]<pi[i2] (resp. pi[i1]<pi[i2]<pi[i3]). 
# The number of (12) patterns of pi, #12(pi), is the number of pairs 1<=i1<i2<=n
# such that pi[i1]<pi[i2]. 

#    Our main result is the following.  Let f_r(n) denote the number of 
# permutations of n letters that have no (132) patterns and exactly r (123)
# patterns.  Then the generating function for the {f_r(n)} is

# THEOREM 1: sum(sum(f_r(n)*z^n*q^r,n=0..infinity),r=0..infinity) = 
#
#		 		      1
#				______________
#				
#				         z
#				1 -  ____________
#				    
#				            z
#				     1 - ____________
#				     
#				               z*q
#				         1 - ___________   
#				         
#				                    z*q^3
#				             1 - ___________
#				             
#				                        z*q^6
#				                 1 - ___________
#				                 
#				                         . . .  															
#
# where the n^th numerator is z*q^(binomial(n,2)).  It is remarkable that such
# a continued fraction encodes information about (132)-avoiding permutations. 

#    		             1. THE PATTERNS
 
#    Define the weight of a permutation pi, Weight(pi), to be: 
# Weight(pi) := z^(nops(pi))*q^(#123(pi))*t^(#12(pi)).
# Let P(q,z,t) be the sum of the weights of all (132)-avoiding permutations.
# We will now derive a functional equation for P(q,z,t).

#    If pi is a (132)-avoiding permutation on {seq(i,i=1..n)}, (n>0) and 
# the largest element, n, is at the k^th position, i.e. pi[k]=n, by letting
# pi1:=[op(1..k-1,pi)] and pi2:=[op(k+1..n,pi)], we have that every element in
# convert(pi1,set) must be larger than every element of convert(pi2,set),
# or else a (132) would be formed, with the n serving as the `3' of the (132).
# Hence, convert(pi1,set)={seq(i,i=n-k+1..n-1)}, and convert(pi2,set)=
# {seq(i,i=1..n-k)}. Furthermore, pi1 and pi2 are (132)-avoiding on their
# own right. Conversely, if pi1 and pi2 are (132)-avoiding permutations on
# {seq(i,i=n-k+1..n-1)} and {seq(i,i=1..n-k)} respectively, (for some k,
# 1<=k<=n) then [op(pi1),n,op(pi2)] is a NON-EMPTY (132)-avoiding permutation.

# We have:

#        #123(pi) = #123(pi1) + #123(pi2) + #12(pi1), 

# since a (123) pattern in pi=[op(pi1),n,op(pi2)] may either be totally immersed 
# in the pi1 part, or wholly immersed in the pi2 part, or may be due to the n 
# serving as the `3' of the (123), the number of which is the number of (12) 
# patterns in pi1.

# We also have 

#        #12(pi) = #12(pi1) + #12(pi2) + nops(pi1), 
 
# and, of course

#        nops(pi) = nops(pi1) + nops(pi2) + 1

# Hence:
# Weight(pi)(q,z,t):= q^(#123(pi))*z^(nops(pi))*t^(#12(pi))
# 		    = q^(#123(pi1)+#123(pi2)+#12(pi1))*z^(nops(pi1)+nops(pi2)+1)*
# 		       	t^(#12(pi1)+ #12(pi2)+nops(pi1))
#		    = z*(q^(#123(pi1))*(q*t)^(#12(pi1))*(z*t)^(nops(pi1))*
#                      	q^(#123(pi2))*t^(#12(pi2))*z^(nops(pi2))
# 		    = z*Weight(pi1)(q,z*t,t*q)*Weight(pi2)(q,z,t). 

# So, if nops(pi)>0, pi is (132)-avoiding, and pi=[op(pi1),max(op(pi)),op(pi2)],
# then pi1, pi2 are smaller (132)-avoiding permutations, and:

#	    Weight(pi)(q,z,t) = z*Weight(pi1)(q,z*t,t*q)*Weight(pi2)(q,z,t)

# Now sum over ALL possible (132)-avoiding permutations pi, to get 

#                   P(q,z,t) = 1 + z*P(q,z*t,t*q)*P(q,z,t)

# (the 1 corresponds to the empty permutation: [])

#    Now let Q(q,z,t) be the sum of all the weights of all permutations with 
# exactly ONE (132) pattern.  By adapting the argument from Miklos Bona's paper 
# (Discrete Math  181 (1998), 267-274) we easily see that Q(q,z,t) satisfies:
 
#    Q(q,z,t) = z*P(q,z*t,q*t)*Q(q,z,t) + z*Q(q,z*t,q*t)*P(q,z,t) +
#                  t^2*z^2*P(q,z*t,q*t)*(P(q,z,t)-1) .

#    This holds since our sole (132) pattern can appear in the elements (1) 
# before n, (2) after n, or (3) with n as the `3' in the (132) pattern.  
# z*P(q,z*t,q*t)*Q(q,z,t) corresponds to (1),  z*Q(q,z*t,q*t)*P(q,z,t) 
# corresponds to (2), and t^2*z^2*P(q,z*t,q*t)*(P(q,z,t)-1) corresponds to 
# (3). We see that case (3) follows since  pi=[op(pi1),n-k,n,op(pi2)], where
# convert(pi1,set)={seq(i,i=n-k+2..n-1)} and convert(pi2,set)=
# {seq(i,i=1..n-k-1)} union {n-k+1}, AND k cannot be equal to n.

#                            2. THE FRACTIONS

#    We now study the generating function P(q,z,t) further, and find that it is
# a pretty continued fraction, and discover a fairly explicit form for its 
# numerator and denominator.  First note that

#		 P(q,z,t) = 1/(1 - z*P(q,z*t,t*q)),  (**)

# so by iteration we have that

#		 		      1
#		 P(q,z,t) = 	________________
#				
#				           z
#				1 -  ______________
#				    
#				              z*t
#				     1 - ________________
#				     
#				                 z*t^2*q
#				         1 -  ________________   
#				         
#				                     z*t^3*q^3
#				              1 -  ________________
#				             
#				                         z*t^4*q^6
#				                   1 -  _____________
#				                 
#				                             . . .  
	
#    By letting P(q,z,t) = A(q,z,t)/B(q,z,t), substitution into (**) shows that
# A(q,z,t) = B(q,z*t,q*t), and B satisfies the functional equation
# B(q,z,t) = B(q,z*t,q*t) - z*B(q,t^2*z,t^2*q).

#    As far as the computations have gone, the series B(q,z,t) appears to have only
# coefficients of 0,1,-1.

#    To find out more about its form we write

#		 B(q,z,t) = sum(phi_m(q,t)z^m,m=0..infinity).

# Then phi_0 = 1, and

#		 phi_m(q,t) = t^m*phi_m(q,q*t)-t^(2*m-2)*q^(m-1)*phi_{m-1}(q,q*t^2),

# for m=1,2,3, ... .  It is easy to see by induction that

#		 phi_m(q,t) = - sum(t^(j*m-2)*q^(m*binomial(j,2)-2*j+3))*
#				    phi_{m-1}(q,t*q^j), j=2..infinity).
				    
# For example, we have

#		 phi_1(q,t) = - sum(t^j*q^(binomial(j,2)),j=0..infinity),

# and

#   phi_2(q,t) = sum(sum(t^(l+2*j-4)*
#	         q^((1/2)*l^2+j^2+l*j-(5/2)*j+6),l=2..infinity),j=2..infinity)

#    In general, the exponent of t in phi_m(q,t) will be a linear form in the
# m summation indices, plus a constant, and the exponent of q will be an affine
# form in these indices, i.e., a quadratic form plus a linear form plus a
# constant.  Let's find all of these forms explicity.

#    Hence suppose in general that

#  phi_m(q,t) = (-1)^m*sum(t^A_m ^. j + b_m)*q^((j,Q^m j)+C_m ^. j + d_m),j),
				   
# in which j is the m-vector of summation indices, Q_m is a real symmetric 
# m x m matrix to be determined, A_m, C_m are m-vectors, and b_m, d_m are
# scalars.  Inducively we find that

#		 A_m = [seq(r,r=1..m)]
#		 b_m = -2m
#		 C_m = [seq((-5/2)*r,r=1..m)]
#		 d_m = 3m

# The m x m matrix Q_m is {min(r,s)/2}_{r,s=1}^{m}.  Thus; we have the following 
# formula for B.

# THEOREM 2:  The denominator B(q,z,t) fo the grand generating function
# P(q,z,t) is explicitly given by

#	B(q,z,t) = 1 + sum((-z*q^3/t^2)^m SUM(t^(sum(rj_r,r=1..m) *
#		           q^(1/2)*(sum(sum(min(r,s)j_r*j_s,s=1..m),r=1..m) - 
#		           5*sum(rj_r,r=1..m))),m=1..infinity),
		          
# where SUM is over j_1,j_2,...,j_m >= 2.

#                       3. THE SERIES COMPUTATIONS

#    If f_r(n) denotes the number of permutations of n letters that contain no
# pattern (132) and have exactly r (123)'s, we write 
# AR(r,z):=sum(f_r(n)z^n,n=1..infinity).  Then AR(r,z) is the coefficient of q^r
# in the series development of P(q,z,1).  That is, we have

#		 		      1
#				_____________  = sum(AR(r,z)*q^r,r=0..infinity)
#				
#				         z
#				1 -  ___________
#				    
#				            z
#				     1 - ___________
#				     
#				               z*q
#				         1 - __________   
#				         
#				                    z*q^3
#				             1 - __________
#				             
#				                        z*q^6
#				                 1 - ___________
#				                 
#				                         . . .  


# From above we see that if we termintate the fraction P(q,z,1) at the numerator
# q^N, say, then we'll know all fo the {AR(r,z)}_{r=0}^{N} exactly.

#    Further, if we know the denominator B(q,z,t) exactly through terms of
# order q^N , then by carrying out the division and keeping the same accuracy, 
# we will, after setting t=1, again obtain all of the generating functions
# {AR(r,z)}_{r=0}^{N} exactly.

#    Finally, to find the denominator B(q,z,t) exactly through terms of order
# q^N, it is sufficient to carry out the iteration that is given by its
# functional equation N times, since further iteration will affect only the
# terms involving powers of q higher than the N^th.

#    In this way we computed the AR(r,z)'s for 0<=r<=15 in a few seconds,
# of which AR(r) for r=0,1,2,3,4,5 are shown below.

#    Let g_r(n) denotes the number of permutations of n letters that contain 
# one (132) pattern and have exactly r (123)'s.  Write 
# Aaron(r,z):=sum(g_r(n)z^n,n=1..infinity).  Then Aaron(r,z) is the coefficient 
# of q^r in the series development of Q(q,z,1).  Since the functional equation 
# for Q can be iterated efficiently (due to the quick calulation of AR(r,z)), 
# we have also computed the Aaron(r,z) for 0<=r<=5 in a few minutes, as
# shown below.
 
# AR(0) = (1-z)/(1-2*z)         
# AR(1) = z^3/(1-2*z)^2  
# AR(2) = z^4*(1-z)/(1-2*z)^3   
# AR(3) = z^5*(z-1)^2/(1-2*z)^4   
# AR(4) = -z^4*(z^5-3*z^4+11*z^3-13*z^2+6*z-1)/(1-2*z)^5
# AR(5) = z^5*(z-1)*(z^5-3*z^4+19*z^3-25*z^2+12*z-2)/(1-2*z)^6
 
# Aaron(0) = z^3/(1-2*z)^2                    
# Aaron(1) = 2*z^5/(1-2*z)^3
# Aaron(2) = -z^4*(z^3-6*z^2+4*z-1)/(1-2*z)^4
# Aaron(3) = 2*z^5*(1-z)*(5*z^2-4*z+1)/(1-2*z)^5
# Aaron(4) = z^6*(z^5+12*z^4-55*z^3+65*z^2-30*z+5)/(1-2*z)^6
# Aaron(5) = -2*z^7*(z^6+6*z^5-40*z^4+80*z^3-69*z^2+27*z-4)/(1-2*z)^7 

##############END OF HUMAN INTERFACE/BEGINNING OF MAPLE PROGRAM##########
 
print(`This is a Maple Package AND article`):
print(`written by Aaron Robertson, Herb Wilf, and Doron Zeilberger.`):
lprint(``):
print(`To find the generating function for the number of permutations with `):
print(`0 (132)-patterns and exactly r (123)-patterns, type:  AR(r); .`):
print(`To find the generating function for the number of permutations with `):
print(`exactly 1 (132)-pattern and exactly r (123)-patterns, type: Aaron(r);.`):
lprint(`` ):
print(`Contains procedures AR(r), Aaron(r), ARSlow(r), and AaronSlow(r).`):
lprint(``):
print(`ARSlow(r) and AaronSlow(r) use the functional equations as well`):
print(`as the partial derivatives to calculate AR and Aaron in a`):
print(`straightforward manner.`):
lprint(``):
print(`WARNING: z is a GLOBAL variable!`):
lprint(``):

# Ders(f,vars,r): given a function f in the list of variables vars, and a
# non-negative integer r finds all partial derivatives of order<=r
Ders:=proc(f,vars,r) local gu,mu,i,j:option remember:
if r<0 then RETURN({})  elif r=0 then RETURN({f}): fi:mu:=Ders(f,vars,r-1):
gu:={f}:for i from 1 to nops(mu) do for j from 1 to nops(vars) do
gu:=gu union {D[j](mu[i])}:od:od:gu:end:
 
# Subs(F,SU): Given a set of expressions F, and a set of
# sets of substations, returns the set of substitutions
Subs:=proc(F,SU) local i,j,gu:gu:={}: for i from 1 to nops(F) do for j from 
1 to nops(SU) do gu:=gu union {F[i](op(SU[j]))}:od:od:gu:end:
 
# Ptor(SFE,SP,vars,SU,r): Given a set of functional equations SFE, phrased
# in terms of the functions in the set of functions  SP in the variables vars,
# and a list of substitutions SU solves for the partial derivatives of the
# functions of P up  to r^th order at the given specializations.
Ptor:=proc(SFE,SP,vars,SU,r) local s,lu, eq,var,i: eq:={}:
for i from 1 to nops(SFE) do eq:=eq union Ders(SFE[i],vars,r): od:
eq:=Subs(eq,SU): var:={}:
for i from 1 to nops(SP) do var:=var union Ders(SP[i],vars,r): od:
var:=Subs(var,SU):solve(eq,var): end:


# Comments for ARSlow(r) and AaronSlow(r):
#    The generating function (in z) for the sequence of cardinalities of 132-
# avoiding permutations with exactly r 123-patterns is
# subs(q=0,diff(P(q,z,1),q$r))/r! , and the analogous quantity for 
# permutations with exactly ONE 132-pattern is
# subs(q=0,diff(Q(q,z,1),q$r))/r! . To find these, we take ALL partial     
# derivatives up to the order r, of the functional equations, then do ALL
# substitutions of [0,z,1],[0,z,0],[0,0,0] for [q,z,t], solve the resulting
# (huge) system of algebraic equations, and then extract the desired quantities.


# ARSlow(r): The (ordinary) generating function in z, of the sequence
# number of 132-avoiding permutations with exactly r 123's
ARSlow:=proc(r) local P,t,q,gu,s,lu: gu:=Ptor({((q,z,t)->P(q,z,t)-1-
z*P(q,z*t,q*t)*P(q,z,t))}, {((q,z,t)->P(q,z,t))},[q,z,t],{[0,z,1],[0,z,0],
[0,0,0]},r):lu:=(q,z,t)->P(q,z,t): for s from 1 to r do lu:=D[1](lu): od:
lu:=lu(0,z,1): RETURN(factor(subs(gu,lu)/r!)): end:
 
# AaronSlow(r): The (ordinary) generating function in z, of the sequence of the 
# number of permutations with exactly ONE 132 and exactly r 123's
AaronSlow:=proc(r) local P,t,q,gu,s,lu,Q: 
gu:=Ptor([((q,z,t)->P(q,z,t)-1-z*P(q,z*t,q*t)*P(q,z,t)),
((q,z,t)->Q(q,z,t)-z*P(q,z*t,q*t)*Q(q,z,t)-z*Q(q,z*t,q*t)*P(q,z,t)
-t^2*z^2*P(q,z*t,q*t)*(P(q,z,t)-1))], 
{((q,z,t)->P(q,z,t)),((q,z,t)->Q(q,z,t)) },[q,z,t],{[0,z,1],[0,z,0],
[0,0,0]},r):lu:=(q,z,t)->Q(q,z,t):for s from 1 to r do 
lu:=D[1](lu): od:lu:=lu(0,z,1): RETURN(factor(subs(gu,lu)/r!)): end:

# numrtr(m):  Returns the numerator of the mth convergent to the continued
# fraction P(q,z,1)
numrtr:=proc(m) option remember;if m=0 then RETURN(1) elif m=1 then RETURN(1-z)
else RETURN(expand(numrtr(m-1)-z*q^(binomial(m,2))*numrtr(m-2))) fi:end:

# dnmntr(m):  Returns the denominator of the mth convergent to the continued
# fraction P(q,z,1)
dnmntr:=proc(m) option remember;
if m=0 then RETURN(1-z) elif m=1 then RETURN(1-2*z)
else RETURN(expand(dnmntr(m-1)-z*q^(binomial(m,2))*dnmntr(m-2))) fi: end:

#AR(r): Returns the function AR(r,z)
AR:=proc(r) local yy,cc;
yy:=series(numrtr(r+1)/dnmntr(r+1),q=0,r+1); cc:=simplify(coeff(yy,q,r));
RETURN(cc); end:

PPP:=proc(r) local yy,cc,G,i;  yy:=series(numrtr(r+1)/dnmntr(r+1),q=0,r+1);
G:=(1-z)/(1-2*z); for i from 1 to r do G:=G+simplify(coeff(yy,q^i))*q^i;od:end:

T:=proc(X,n)  local R: if n=1 then RETURN((X-1)/(z*X)); else
R:=(T(X,n-1)-1)/(z*t^(n-1)*q^(binomial(n-1,2))*T(X,n-1)):fi: simplify(R); end:

#We now iterate the functional equation
#Q(q,z,t)=zQ(q,zt,qt)[P(q,z,t)]^2 + t^2 z [P(q,z,t)-1)]^2,
#which can easily be deduced from the equation given in the
#RWZ paper  then subs in t=1 where TP=P(q,zt,qt)
QQ:=proc(r)  local i,m,Z,PP,tmp,ans,d,tmp2,dd,ans2,W:
if r=0 then ERROR(`r must be positive`): fi: m:=max(3,r): Z[m+1]:=1:
for i from m by -1 to 2  do
W:=subs(t=1,T(PP,i)):
Z[i]:=simplify(z*q^(binomial(i,2))*Z[i+1]*W^2+z*q^(2*i+binomial(i,2))*(W-1)^2);
od:
Z[1]:=simplify(z*Z[2]*subs(t=1,T(PP,1))^2+z*q^2*(subs(t=1,T(PP,1))-1)^2);
Z[0]:=simplify(z*Z[1]*PP^2+z*(PP-1)^2); tmp:=simplify(subs(PP=PPP(r),Z[0]));
d:=denom(tmp):tmp2:=collect(simplify(d*tmp),q):ans2:=factor(coeff(tmp2,q^r));
ans:=denom(simplify(d/ans2));dd:=numer(simplify(d/ans2));[ans,dd]:end:

Aaron:=proc(r)  local i,tmp,ans,d,c,G:option remember:
if r=0 then RETURN(z^3/(1-2*z)^2): fi:tmp:=QQ(r):d:=collect(expand(tmp[2]),q):
c:=tcoeff(d,q):G:=0:for i from 0 to r-1 do G:=G+Aaron(i)*coeff(d,q^(r-i)):od:
ans:=factor(simplify((1/c)*(tmp[1]-G)));end:
 
###################END OF MAPLE PROGRAM/END OF ARTICLE########################


