#OK to post homework
#Blair Seidler, 2021-02-07, Assignment 5

with(combinat):

Help:=proc():
print(`PermuG(n)`): 
print(`Rank(pi)`,`NextPG(pi)`,`PrevPG(pi)`):
print(`MyRandPermA(n)`,`MyRandPermI(n)`):
print(`NoABC(n)`):
end:

# 1.

#PermuG(n): the list of all permutations of {1, ...,n} using Johnson-Trotter
PermuG:=proc(n)  local S,i,pi,pi1,L,j:

option remember:
if n=0 then
 RETURN([[]]):
fi:


S:=PermuG(n-1):

L:=[]:

for i from 1 to nops(S) do
  pi:=S[i]:
  for j from 1 to n do
    if (i mod 2)=0 then
      pi1:=[op(1..j-1,pi),n,op(j..nops(pi),pi)]:
    else
      pi1:=[op(1..n-j,pi),n,op(n-j+1..nops(pi),pi)]:
    fi:
    L:=[op(L),pi1]:
  od:
od:

L:
end:

# 2.

#Rank(pi): Finds the rank of the permutation (algorithm on page 4 of Dennis^2)
Rank:=proc(pi) local n,i,cur,p1,rp1:
option remember:
n:=nops(pi):
if n<=1 then
  RETURN(0):
fi:
for i from 1 to n do
  if pi[i]=n then cur:=i: fi:
od:
p1:=[op(1..cur-1,pi),op(cur+1..n,pi)]:
rp1:=Rank(p1):
if (rp1 mod 2)=1 then
  n*rp1+cur-1:
else:
  n*rp1+n-cur:
fi:
end:

#NextPG(pi): inputs a permutation pi of {1, ...,n} for some n, and outputs the one that comes
#right after it in the above PermuG(n), but of course, without actually constucting this 
#super-exponentially large set.
NextPG:=proc(pi) local cur,r,r1,i,n,P:
n:=nops(pi):
r:=Rank(pi):
if r=n!-1 then
  RETURN(FAIL):
fi:

for i from 1 to n do
  if pi[i]=n then cur:=i: fi:
od:
r1:=Rank([op(1..cur-1,pi),op(cur+1..n,pi)]):

if (r1 mod 2)=0 then
  # Moving right to left
  if cur=1 then
    RETURN([n,op(NextPG([op(2..n,pi)]))]):
  fi:
  RETURN([op(1..cur-2,pi),n,op(cur-1,pi),op(cur+1..n,pi)]):
else
  # Moving left to right
  if cur=n then
    RETURN([op(NextPG([op(1..n-1,pi)])),n]):
  fi:
  RETURN([op(1..cur-1,pi),op(cur+1,pi),n,op(cur+2..n,pi)]):
fi:

end:

#PrevPG(pi): inputs a permutation pi of {1, ...,n} for some n, and outputs the one that comes
#right before it in the above PermuG(n), but of course, without actually constucting this 
#super-exponentially large set.
PrevPG:=proc(pi) local cur,r,r1,i,n,P:
n:=nops(pi):
r:=Rank(pi):
if r=0 then
  RETURN(FAIL):
fi:

for i from 1 to n do
  if pi[i]=n then cur:=i: fi:
od:
r1:=Rank([op(1..cur-1,pi),op(cur+1..n,pi)]):

if (r1 mod 2)=0 then
  # Moving right to left
  if cur=n then
    RETURN([op(PrevPG([op(1..n-1,pi)])),n]):
  fi:
  RETURN([op(1..cur-1,pi),op(cur+1,pi),n,op(cur+2..n,pi)]):
else
  # Moving left to right
  if cur=1 then
    RETURN([n,op(PrevPG([op(2..n,pi)]))]):
  fi:
  RETURN([op(1..cur-2,pi),n,op(cur-1,pi),op(cur+1..n,pi)]):
fi:

end:

# 3.

#MyRandPermA(n): inputs a non-neg. integer n and outputs a random permutation of {1,...,n}
#(uniformly at random) using the pre-pending paradigm, inspired by PermuA(n)
MyRandPermA:=proc(n) local i,r,S,T:
if n=0 then
 RETURN([[]]):
fi:

T:=[]:
S:=[seq(1..n)]:
for i from 1 to n do
  r:=rand(1..n-i+1)():
  T:=[op(T),S[r]]:
  S:=[op(1..r-1,S),op(r+1..n-i+1,S)]:
od:
T:

end:

#MyRandPermI(n): inputs a non-neg. integer n and outputs a random permutation of {1,...,n}
#(uniformly at random) using the insertion paradigm, inspired by PermuI(n)
MyRandPermI:=proc(n) local i,r,T:
if n=0 then
 RETURN([[]]):
fi:

T:=[]:
for i from 1 to n do
  r:=rand(1..i)():
  T:=[op(1..r-1,T),i,op(r..i-1,T)]:
od:
T:
end:

#I have obviously written some terribly inefficient code here. Some timing results:
# n=10000, randperm 0s, MyRandPermA 3.546s, MyRandPermI 2.890s
# n=100000, randperm 0s, MyRandPermA 346s, MyRandPermI 237s
# n=1000000, randperm 0.062s, MyRandPermA/I not happening

# 4.

#NoABC(n): inputs a positive integer n and outputs the set of permutations AVOIDING the pattern 123.
#That is, outputs all permutations which do not have a triple (i,j,k) such that
# i<j<k and pi(i)<pi(j)<pi(k)
NoABC:=proc(n) local T,S,s,i,j,k,f:
#I'm going to try and brute-force this first. It may not run this century on large inputs
S:=PermuA(n):
T:=[]:
for s from 1 to nops(S) do
  f:=1:
  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 S[s][i]<S[s][j] and S[s][j]<S[s][k] then
          f:=0:
	fi:
      od:
    od:
  od:
  if f=1 then
    T:=[op(T),S[s]]:
  fi:
od:
T:
end:

# [seq(nops(NoABC(n)),n=0..8)]=[1, 1, 2, 5, 14, 42, 132, 429, 1430]
# These are the Catalan numbers (OEIS sequence A000108)
# The brute-force method provides me with no insight as to why this is true,
# so I will have to think about it some more.

#### Code included from C5.txt ####
Help5:=proc(): print(`Redu(L), Expa(pi,S), PermuA(n),PermuI(n) `): end:

#Redu(L): inputs a list L of numbers and outputs
#The multi-set permuation of {1, ...,nops({op(L)}) where n is the number of
#number of 
#distinct members of L and the smallest number in {op(L)}
#replaced by 1, the second smallest by 2, etc.
#For example
#Redu([evalf(Pi),evalf(Pi),evalf(exp(1)),evalf(exp(1))]
#should return
#[2,2,1,1]
##Redu([evalf(Pi),evalf(exp(1)),evalf(exp(1)),5]
#=[2,1,1,3]
#
Redu:=proc(L) local S,T,i:
S:=sort(convert({op(L)},list)):

for i from 1 to nops(S) do
 T[S[i]]:=i:
od:

[seq(T[L[i]],i=1..nops(L))]:

end:

#Expa(pi,S): inputs a permutation pi of 1,..., nops(pi) and
#a sorted list of distinct numbers S, replaces 1 by
#the smallest member of S, etc. 
#For example
#Expa([2,1,3],[11,13,15])= [13,11,15]
Expa:=proc(pi,S) local i:
subs({seq(i=S[i],i=1..nops(S))},pi):
end:

#PermuA(n): inputs a non-neg. integer n and outputs the list of
#all n! permutations of {1,...,n} in lex-order (same as permute(n) in combinat)
PermuA:=proc(n) local S,L,j, J,pi:
option remember:
if n=0 then
 RETURN([[]]):
fi:

S:=PermuA(n-1):

L:=[]:


for j from 1 to n do
  #J:=1...j-1,j+1..n
  J:=[seq(1..j-1), seq(j+1..n)]:
   L:= [op(L),seq([ j, op(Expa(S[i],J)) ],i=1..nops(S))] :
od:
L:
end:

#PermuI(n): the list of all permutations of {1, ...,n} using inertion
PermuI:=proc(n)  local S,i,pi,pi1,L,j:

option remember:
if n=0 then
 RETURN([[]]):
fi:


S:=PermuI(n-1):

L:=[]:

for i from 1 to nops(S) do
  pi:=S[i]:
  for j from 1 to n do
    pi1:=[op(1..j-1,pi),n,op(j..nops(pi),pi)]:
    L:=[op(L),pi1]:
  od:
od:

L:
end:

#### End of code included from C5.txt ####