# OK to post

restart:


read("C4.txt"):
read("C5.txt"):

# 1. Done

# 2. Conjecture and prove:
#    factor(WtE(powerset(n), S->nops(S), x));



PS := n -> powerset({seq(i,i=1..n)}):

# Try n=1..6:
[seq( factor(WtE(PS(n), T->nops(T), x)), n=1..6 )];

# Output (should be):
# [(1+x), (1+x)^2, (1+x)^3, (1+x)^4, (1+x)^5, (1+x)^6]

# Conjecture:
# factor(WtE(PS(n), T->nops(T), x)) = (1+x)^n.

# Proof:
# WtE(PS(n),T->nops(T),x) = Sum_{A subseteq [n]} x^{|A|}.
# Group subsets by size k. There are binomial(n,k) subsets of size k, so
#   Sum_{A subseteq [n]} x^{|A|}
# = Sum_{k=0..n} binomial(n,k) x^k
# = (1+x)^n   (binomial theorem).
# QED.


# 3. Write NumPat3(pi,sig) using redu (from C4.txt)


NumPat3 := proc(pi, sig)
local n, i, j, k, co;

# sig must be a permutation of length 3
if not (type(sig, list) and nops(sig)=3 and {op(sig)}={1,2,3}) then
  return(FAIL);
fi;

n := nops(pi);
co := 0;

for i from 1 to n-2 do
  for j from i+1 to n-1 do
    for k from j+1 to n do
      if redu([pi[i],pi[j],pi[k]]) = sig then
        co := co + 1;
      fi;
    od;
  od;
od;

return(co);
end:

# Test example:
NumPat3([2,1,3,4],[1,2,3]);   # should output 2


# 4. Define L and identify the OEIS sequences for fixed j


L := [seq( WtE(permute(n), pi->NumPat3(pi,[1,2,3]), x), n=1..7 )]:

# My output for L (as polynomials):
# n=1: 1
# n=2: 2
# n=3: 5 + x
# n=4: 14 + 6*x + 3*x^2 + x^4
# n=5: 42 + 27*x + 24*x^2 + 7*x^3 + 9*x^4 + 6*x^5 + 4*x^7 + x^10
# n=6: 132 + 110*x + 133*x^2 + 70*x^3 + 74*x^4 + 54*x^5 + 37*x^6 + 32*x^7
#       + 24*x^8 + 12*x^9 + 16*x^10 + 6*x^11 + 6*x^12 + 8*x^13 + 5*x^16 + x^20
# n=7: 429 + 429*x + 635*x^2 + 461*x^3 + 507*x^4 + 395*x^5 + 387*x^6 + 320*x^7
#       + 260*x^8 + 232*x^9 + 191*x^10 + 162*x^11 + 104*x^12 + 130*x^13 + 100*x^14
#       + 24*x^15 + 74*x^16 + 62*x^17 + 18*x^18 + 32*x^19 + 10*x^20 + 30*x^21
#       + 13*x^22 + 8*x^23 + 10*x^25 + 10*x^26 + 6*x^30 + x^35

# For fixed j, the assignment asks for the sequence (in n):
SeqJ := j -> [seq(coeff(L[i],x,j), i=1..nops(L))]:

SeqJ(0);  # [1,2,5,14,42,132,429]
SeqJ(1);  # [0,0,1,6,27,110,429]
SeqJ(2);  # [0,0,0,3,24,133,635]

# OEIS identification:
# The whole table (rows n, columns j) is the triangle A138159.
# In particular, the OEIS page states:
#   j=0 column: A000108 (Catalan numbers) evaluated at n
#   j=1 column: A003517 evaluated at n-1
#   j=2 column: A001089 evaluated at n
#
# Note: A138159 is stated for pattern 321, but #occurrences of 123 and 321
# have the same distribution by a symmetry bijection (reverse+complement).
# Hence our L for 123 matches that triangle.


# 5. Repeat with the other 5 patterns of length 3:
#    [1,3,2], [2,1,3], [2,3,1], [3,1,2], [3,2,1]


# There are only TWO distinct distributions (by symmetries of permutations):

# (A) Patterns 123 and 321 have the SAME distribution.
# Check quickly:
L321 := [seq( WtE(permute(n), pi->NumPat3(pi,[3,2,1]), x), n=1..7 )]:
evalb(L321 = L);   # should be true

# OEIS: triangle A138159 (same as in Q4).

# (B) The other four patterns 132, 213, 231, 312 all have the SAME distribution.
L132 := [seq( WtE(permute(n), pi->NumPat3(pi,[1,3,2]), x), n=1..7 )]:
L213 := [seq( WtE(permute(n), pi->NumPat3(pi,[2,1,3]), x), n=1..7 )]:
L231 := [seq( WtE(permute(n), pi->NumPat3(pi,[2,3,1]), x), n=1..7 )]:
L312 := [seq( WtE(permute(n), pi->NumPat3(pi,[3,1,2]), x), n=1..7 )]:

evalb(L132=L213 and L213=L231 and L231=L312);  # should be true

# First few row polynomials for this class (e.g. 312):
# n=4: 14 + 5*x + 4*x^2 + x^3
# n=5: 42 + 21*x + 23*x^2 + 14*x^3 + 12*x^4 + 5*x^5 + 3*x^6
# n=6: 132 + 84*x + 107*x^2 + 82*x^3 + 96*x^4 + 55*x^5 + 64*x^6 + ...
#
# OEIS: triangle A263771 (pattern 312), and OEIS states it also covers 132, 213, 231.
# The OEIS page also lists column A-numbers:
#   k=0 column: A000108
#   k=1 column: A002054(n-2) for n>=3
#   k=2 column: A082970
#   k=3 column: A082971
#   k=4 column: A138162
#   k=5 column: A138163