#zeta-orbits.maple
interface(quiet=true):
with(linalg):
##############################################################
#
#   This file computes the number of disjoint \zeta in S_n
#   with full support, modulo the trivial identites:
#   conjugation with the longest element, taking inverse, and 
#   cyclic shift.
#   It tabulates their rank in the new order, and k, which 
#   is the cardinality of up(\zeta).
#   
#   Has a flaw, in that it doesn't eliminate all permutations 
#   with more than 2 components which are disjoint, hence the 
#   output must be checked by hand.  For example, these are among
#   the permutations output for n=8:
#
#	n    k      cycles
#	10   4   {[3, 6], [4, 7], [1, 2], [5, 8]}  
#	5    3   {[3, 7], [4, 5, 6], [1, 2, 8]} 
#	8    4   {[1, 3], [2, 4], [5, 7], [6, 8]} 
#
#	Note in the first, the cycle [1, 2], which is disjoint from
#	the others.  Similarly, in the second, there is the cycle
#	[4, 5, 6], also disjoint from the others.
#	The third one has 2 components, {[1, 3], [2, 4]} and 
#	{[5, 7], [6, 8]}  which are disjoint.  If this is ever properly
#	implements in C, then this bug should be fixed.
#
#	Frank Sottile
#	19 May 1997
#
###############################################################
n:=4:

#
#	Finds the inverse of a cycle
#
INV:= proc(C,n)  local Cp,Cpp,m,mm,up,j:;

m:=nops(C):

up:=[]:
for j from 2 to m do  up:=[up[],j]: od:
up:=[up[],1]:
Cp:=[]:

for j from 1 to nops(C) do
Cp:= [Cp[],C[nops(C)+1-j]]:  od:

mm:= Cp[1]:
for j from 2 to m do
 mm:= min(mm,Cp[j]): od:

while not(mm=Cp[1]) do
  Cpp:=eval(Cp):
  Cp:=[]:
  for j from 1 to m do
   Cp:=[Cp[],Cpp[up[j]]]: od:
od:
eval(Cp);

end:

#
#	Conjugates a cycle by longest element
#
BAR:= proc(C,n)  local Cp,Cpp,m,mm,up,j:;

m:=nops(C):

up:=[]:
for j from 2 to m do  up:=[up[],j]: od:
up:=[up[],1]:
Cp:=[]:

for j from 1 to nops(C) do
Cp:= [Cp[],n+1-C[j]]:  od:

mm:= Cp[1]:
for j from 2 to m do
 mm:= min(mm,Cp[j]): od:

while not(mm=Cp[1]) do
  Cpp:=eval(Cp):
  Cp:=[]:
  for j from 1 to m do
   Cp:=[Cp[],Cpp[up[j]]]: od:
od:
eval(Cp);

end:

#
#	Performs a cyclic shift on a given cycle
#
CyclicShift:= proc(C,n,down)  local Cp,Cpp,m,mm,up,j:

m:=nops(C):

up:=[]:
for j from 2 to m do  up:=[up[],j]: od:
up:=[up[],1]:


Cp:=[]:
for j from 1 to m  do
Cp:=[Cp[],down[C[up[j]]]]: od:

mm:= Cp[1]:
for j from 2 to m do
 mm:= min(mm,Cp[j]): od:

while not(mm=Cp[1]) do
  Cpp:=eval(Cp):
  Cp:=[]:
  for j from 1 to m do
   Cp:=[Cp[],Cpp[up[j]]]: od:
od:
eval(Cp);

end:

down:=[n]:
N:={1}:
for j from 2 to n do 
N:=N union {j}: down:=[down[],j-1]: od:


Sn:=[]:
pi:=[]:                            
for j from 1 to n do:       #  Initializes pi = identity
pi := [pi[],j]: od:         

asc := 1:                   #  The ASCents of pi
while (asc > 0) do:

w := eval(pi):              #   the value of pi
Sn:=[Sn[],w]:
asc := 0;
for j from 1 to n-1 do:                   #   Finds the last ascent of pi
if (pi[j]<pi[j+1]) then asc := j fi: od;   

if (asc >0) then                            # If pi has an ascent  at asc...
  dummy := eval(pi[asc]);
  for ii from asc+1 to n do:
  if (pi[ii]>dummy) then bigger:=ii fi: od:
  pi[asc]:= eval(w[bigger]):                # kills the ascent at asc
  w[bigger]:= dummy: 
  for jj from 1 to n-asc do;
  pi[asc+jj]:= eval(w[n+1-jj]); od;         # makes pi increasing afterwards
fi;
od:


#print(Sn);
#print(N);

Cycles:=[]:


for zeta in Sn do
 full:= true:
 j:= 0:
 while (full) and (j<n) do
  j:= j+1:
  if (zeta[j]=j) then full:= false: fi: od:
 
 if (full) then 
 Cycle:={}:
 notmoved:=eval(N):

 while not(notmoved ={}) do
  i:= notmoved[1]:
  notmoved:= notmoved minus {i}:
  cycle:=[i]:
  j:= zeta[i]:
  while not(j=i) do
  notmoved:= notmoved minus {j}:
   cycle:=[cycle[],j]:
   j:= zeta[j]:
  od:
  Cycle:=Cycle union {cycle}:
  od:
Cycles:=[Cycles[],Cycle]:
 fi:

od:

#
#	Doesn't quite eliminate all disconnected permutations
#
#
#

Connected:=[]:

for ii from 1 to  nops(Cycles) do

if (nops(op(ii,Cycles))=1) then 
  Connected:=[Connected[],op(ii,Cycles)]:
else
 Comp:= nops(op(ii,Cycles)):
 for A from 1 to Comp-1 do

  for B from 2 to Comp do
   connected:= false:
   C1:= op(A,op(ii,Cycles)):
   C2:= op(B,op(ii,Cycles)):
   for i from 1 to nops(C1)-1 do
    for j from 1 to nops(C2)-1 do
     a:=op(i,C1):  al:=op(i+1,C1):
     b:=op(j,C2):  be:=op(j+1,C2):
     al:= max(al-a,sign(al-a)*(a-al-n)):
     b := max(b-a,sign(b-a)*(a-b-n)):
     be:= max(be-a,sign(be-a)*(a-be-n)):
     if (((b-al)*(be-al))<0) then connected:=true:   fi:
     od: od:
     if (connected) then 
Connected:=[Connected[],op(ii,Cycles)]:fi:
od: od: fi:
od:

Types:={Connected[1]}:
for ii from 2 to nops(Connected) do 
 C:=op(ii,Connected):
  Cb:={};  CI:={}; Cbi:={};
  for i from 1 to nops(C) do
  Cb := Cb  union {BAR(op(i,C),n)}: 
  CI := CI  union {INV(op(i,C),n)}: 
  Cbi:= Cbi union {BAR(INV(op(i,C),n),n)}: 
  od:
 shifts:={eval(C),eval(Cb),eval(CI),eval(Cbi)}:

 for jj from 1 to n-1 do
  Cp:={};  Cb:={};  CI:={}; Cbi:={};
  for i from 1 to nops(C) do
  Cp:= Cp union {CyclicShift(op(i,C),n,down)}: od:
  C:= eval(Cp): 
  for i from 1 to nops(Cp) do
  Cb := Cb  union {BAR(op(i,Cp),n)}: 
  CI := CI  union {INV(op(i,Cp),n)}: 
  Cbi:= Cbi union {BAR(INV(op(i,Cp),n),n)}: 
  od:
  shifts:=shifts union {eval(Cp),eval(Cb),eval(CI),eval(Cbi)}:

  od:

if ((shifts intersect Types)={}) 
then Types:=Types union {op(ii,Connected)}: fi:
od:

M:= nops(Types):

Pairs:=matrix(1,M,0):

for zeta in Types do



pie:=array[1..n,0];
  for i from 1 to nops(zeta) do
   for j from 1 to nops(op(i,zeta))-1 do
    pie[op(j,op(i,zeta))]:= op(j+1,op(i,zeta)): od:
   pie[op(nops(op(i,zeta)),op(i,zeta))]:= op(1,op(i,zeta)): od:

pi:=[pie[1]]:
for i from 2 to n do
 pi:=[pi[],pie[i]]: od:


rank := 0:

up := {}:   
down := {};
for j from 1 to n do:
 if (j<pi[j]) then up:=up union {j} fi:
 if (j>pi[j]) then down:=down union {j} fi:
od:

for a in up do:   for b in down do:
if (pi[a]>pi[b]) then rank := rank +1  fi:
if (a>b) then rank := rank -1 fi:
od:od:

for a in up do:     for aa in up do:
if (a>aa and pi[a]<pi[aa]) then  rank := rank -1 fi:
od:od:

for a in down do:   for aa in down do:
if (a>aa and pi[a]<pi[aa]) then  rank := rank -1 fi:
od:od:

lprint(rank,min(nops(up),n-nops(up)),zeta);   #,pi);

od:


time();
quit;