(* Posets.nb

Original design and implementation (Version 1.9) by 
Curtis Greene and Eugenie Hunsicker, Summer 1990/91.

Updates by Curtis Greene, John Dollhopf, Sam Hsiao, 1991-94:

Comments, questions, bug reports to:  
Curtis Greene
Department of Mathematics
Haverford College
Haverford, PA 19041
cgreene@haverford.edu  
*)

(* Changes 1.9 -> 2.0       (October 1991/ Curtis Greene) :                    
                                            
         *: Framed -> Frame (for V.2.0 of Mathematica)
         *: RandomP[n] -> RandomP[n,p]   
         *: PlotColor -> ColorOutput (for V.2.0 of Mathematica)
         *: SubPoset expanded to allow Boolean  selection of indices.
         *: RestoreSubPosetLabels added
         *: RestoreDualPosetLabels added    *)

(* Changes 2.0->2.0z       (June 1993/ Curtis Greene):   
         *: Fixed bug in StrongS[n]. (Noticed by Francesco Brenti.
         *: Some undocumented commands removed (PlaneP,SpaceP,
             PartForm, SolidForm,RDCode). These are available on request. 
*)

 (* Changes 2.0z->2.2a      (July 1993/ John Dollhopf, Curtis Greene):
         *: Links renamed CoverRelations to avoid conflict with 
               Mathematica  V2.1 Links command.
         *: Build[name] now prints out rank sizes in Verbose mode.
         *: Contractions[graph,n] added.
         *: CharPoly[name,q] added.
         *: JoinSubLattice[name,subset], JoinSubLatticeQ[name,subset] added.
         *: MeetSubLattice[name,subset], MeetSubLatticeQ[name,subset] added.
         *: SubLatticeQ[name,subset] added.   *)

 (* Changes 2.2a->2.2b       (August 1993/ Curtis Greene):
        *: Commands renamed:
        *:       ToLinks -> ToRelation
        *:       JILinks -> JICoverRelations, etc.
        *:       Covers -> Successors, CoveredBy ->Predecessors
        *:       MaxChainsBelow -> MaximalChainsDown
        *:       MaxChainsAbove -> MaximalChainsUp    
        *:  LocateIndices[name,set] added .    *)

 (* Changes 2.2b->2.2c       (October 1993/ Curtis Greene):
        *: ThinLines=True  ->  Thinner=p
        *:        (To fix problems with Mathematica  v.2.2 menu command)   *)

 (* Changes 2.2c->2.2d       (July 1994/ Sam Hsiao, Curtis Greene):
         *:  Fixed bug in Build (for 3-element posets).
         *:  Fixed bug in RankedQ (for non-ranked posets).
         *:             (Required changing Build and TopSort.)
         *:  Build now uses internal commands TSort1, TSort2.
         *:  Use of "automatic" removed from Build.
         *:  New public command TopSort[links,N].
         *:  New public command LongestChain[links,N] (previously called 
HeightFunction).
         *:  Fixed bug in TClosure (to handle duplicates in cover relation).
         *:  Fixed bug in NumP (to eliminate duplicates in cover relation).
         *:  PosetP[name] added.  (Partitions of a poset into chains.)  
         *:  ChainsBetweenGF[name,a,b,q] added, replacing ChainsBelow.  
         *:  ZetaPoly[name,a,b,q] and ZetaPoly[name,q] added, replacing 
Z[name].  
         *:  OmegaGF[name,q] and OmegaBarGF[name,q] fixed.  (Degree of 
numerator wrong.)
         *:  IntervalP[name,bottom,top] added.
         *:  Sublattice[name,points] added.
         *: JoinSubLattice[name,points], MeetSubLattice[name,points] changed. 
Now return
         *:        only a list of indices (input for SubPoset).
         *:  DistributiveLatticeQ[name] added.
         *:  InversePermutation[list] added.
         *:  ListCommands[] and ListCommands[string] added.
         *:  LocateIndices[name,set] replaced by Locate[name,x] and 
LocateSet[name,set].
         *:  Vee[name,x,y],Wedge[name,x,y] added (lattice operations).
         *:  Commands renamed:
         *:       Predecessors -> CoCovers
         *:       Successors ->  Covers
         *:       ChainsBelow ->  ChainsBelowGF
         *:       Contractions -> ContractionLattice
         *:       Thinner -> Thinness
         *:       Contractions -> ContractionLattice
         *:  New abbreviations:
         *:       JSL == JoinSubLattice
         *:       MSL == MeetSubLattice
         *:       SL == SubLattice
         *:  NumPosets[n], NumLattices[n] added.
         *:  Small changes to IsomorphicQ.
               *)


 (* More work needed):
        *: Zap command doesn't do all that is claimed.
        *)

Off[General::spell]
Off[General::spell1]

BeginPackage["Posets`"];


                    (* USAGE *)


(* Usage:: Standard poset definitions *)

Subsets::usage = "Subsets[n_] defines the poset of subsets
of an n-element set, ordered by containment.";

Subwords::usage = "Subwords[w_] defines the poset of 
subwords of a word w.";

SetP::usage = "SetP[n_] defines the poset of partitions of 
an n-element set, ordered by refinement. Notation: partitions
are represented by listing, for each element, the smallest
element in its block.";

NonXP::usage = "NonXP[n_] defines the poset of non-crossing
partitions of an n-element set, ordered by refinement.";

WeakS::usage = "WeakS[n_] defines the weak order on the
symmetric group on n elements.  Coverings correspond to
transpositions of adjacent increasing pairs.";

StrongS::usage = "StrongS[n_] defines strong (Bruhat) order
on the symmetric group of n elements."; 

Div::usage="Div[N] defines the poset of divisors of
the integer N.";

NumP::usage = "NumP[N_] defines the poset of partitions of the 
integer N, ordered by refinement.";

YoungsLattice::usage = "YoungsLattice[lambda] defines 
Young's Lattice for the partition lambda, i.e., the poset of 
partitions whose Ferrers diagram fits into that of lambda."; 

MajP::usage = "MajP[N] defines the poset
of partitions of an integer N, ordered by majorization.";

Chain::usage = "Chain[N] defines a chain with N elements";

Antichain::usage="Antichain[N] defines an antichain with
N elements";

ZigZag::usage="ZigZag[N] defines a zigzag poset with 
N elements";

RandomP::usage="RandomP[N,p] defines a random poset with
N elements, generating related pairs with probability p 
and taking the transitive closure.";

ContractionLattice::usage = "ContractionLattice[{list of edges},N] 
defines the lattice of contractions of a graph G with N
vertices.";

PosetP::usage = "PosetP[name] defines the poset whose 
elements are partitions of P[name] into a chains.
The partitions are ordered by refinement. Build should be 
applied next. Input may be in any of the three standard forms." ;

JofP::usage = "JofP[name] returns the poset of order ideals
of P[name]. Output is in {f,minlist,H} form, and Build should be 
applied next. Input may be in any of the three standard forms." ;

(* Usage:: Commands *)

Build::usage = "Build[{coverfunction,min,height},name],
Build[matrix,name], or Build[{relations,card},name}] 
generate the fundamental objects P, Rank, and CoverRelations 
for the poset 'name'.  If <name> is omitted, a unique name of 
the form Poset1, Poset2, etc is generated.";

P::usage = "P[name] is a list of the elements in poset 'name'.";

Rank::usage = "Rank[name][[x]] is the rank of P[name][[x]].";

CoverRelations::usage = "CoverRelations[name] is a list of 
the cover relations in P[name].  Pair {a,b} is present 
if the bth element of P[name] covers the ath element.";

PGraded::usage = "PGraded[name][[k]] is a list of the elements 
at rank k in P[name].";

Card::usage = "Card[name] is the cardinality of P[name].";

NK::usage = "NK[name][[k]] is the number of elements 
at rank k in P[name].";

H::usage = "H[name] is the height of P[name].";

RankedQ::usage = "RankedQ[name] tests to see if a poset is
unranked, ranked, or strongly ranked (all minimal elements 
have rank zero). If unranked, it defines (or redefines) the
objects Rank, P, PGraded, NK, and H using length of longest 
chain as rank.";

RGF::usage = "RGF[name] is the rank generating function
for P[name], represented as a polynomial in q.";

SortByRanks::usage = "SortByRanks[name] reorders the elements 
at each rank of P[name] in lexicographic order.";

Relabel::usage = "Relabel[name,f] relabels the elements of 
P[name] by applying the function f to each element of P[name].
For example, Relabel[name,Compact] makes labels into strings. 
Relabel[name] relabels P[name] using the numbers 
1,2,...,Card[name] as labels.";

Compact::usage="Compact[w] returns a string corresponding to the 
list w.  Relabel[name,Compact] is used to produce a more readable
display.";

Diagram::usage = "Diagram[name] draws the Hasse diagram of 
P[name].  Diagram[name,points,links,options] shows the diagram 
with special points and links highlighted. Options include 
ShowLabels, ThinLines, Jiggled.";

ShowLabels::usage = "ShowLabels -> False (default) causes 
element labels to be suppressed  in the diagram.";

Thinness::usage = "Thinness -> 1 (default) causes thick 
lines and dots to be used in the diagram. Thinness -> p causes
line and dot thickness to be multiplied by p.";

Jiggled::usage = "Jiggled -> 0 (default) displays a diagram 
with uniform alignment.  Jiggled -> k causes the position of 
each element to be perturbed by an amount proportional to k."; 

MaxAntichain::usage =  "MaxAntichain[name] is a list of the 
indices of elements in a maximum-sized antichain in P[name]";

DilworthCover::usage = "DilworthCover[name] is a list of 
links in a minimum-sized chain covering of P[name]."; 

Fuse::usage = "Fuse[links] assembles links into chains (if 
possible)";

ZetaP::usage = "ZetaP[name] is the incidence matrix of P[name].";

Mu::usage = "Mu[name] is the Mobius matrix of P[name]
(i.e., the inverse of ZetaP).";

Up::usage = "Up[name][[x]] is a list of the elements 
greater than or equal to P[name][[x]].";

Down::usage = "Down[name][[x]] is a list of the 
elements less than or equal to P[name][[x]].";

Covers::usage="Covers[name][[x]] is a list of
the elements covering P[name][[x]].";

CoCovers::usage="CoCovers[name][[x]] is a list
of the elements covered by P[name][[x]].";

UpDegree::usage="UpDegree[name][[x]] is the number of 
elements covering P[name][[x]].";

DownDegree::usage="DownDegree[name][[x]] is the number of 
elements covered by P[name][[x]].";

MaximalChainsDown::usage = "MaximalChainsDown[name][[x]] is the
number of maximal chains descending from P[name][[x]].";

MaximalChainsUp::usage = "MaximalChainsUp[name][[x]] is the
number of maximal chains ascending from P[name][[x]].";

ChainsBetweenGF::usage = "ChainsBetweenGF[name,a,b] is a 
generating function (in q) which enumerates (by length) the
chains from a to b in P[name].";

LongestChain::usage = "LongestChain[links,N] returns a list 
showing the length of the longest chain terminating in each element.
Assumes <links> is an acyclic relation.";

ZetaPoly::usage = "ZetaPoly[name,a,b,n] gives the Zeta polynomial
of the interval [a,b], i.e., the number of multichains
a=x_0 <= x_1 <= ... <= x_n = b with endpoints a and b.  
ZetaPoly[name,n] gives the Zeta polynomial of the poset P[name], 
i.e., the number of multichains x_1 <= x_2 <= ... <= x_(n-1)."; 

CharPoly::usage = "CharPoly[name,q] is the characteristic 
polynomial (= generalized chromatic polynomial) of P[name], 
using the variable q.";

LatticeQ::usage = "LatticeQ[name] returns True if P[name] is 
a lattice, and False otherwise.";

DistributiveLatticeQ::usage = "DistributiveLatticeQ[name] returns 
True if P[name] is a distributive lattice, and False otherwise.";  

Vee::usage="Vee[name,x,y] returns the lattice-join of x and y.
Input and output are labels, not indices.";

Wedge::usage="Wedge[name,x,y] returns the lattice-meet of x and y.
Input and output are labels, not indices.";

LJoin::usage = "LJoin[name] is the join table for P[name], 
i.e., LJoin[name][[a,b]] is the join of P[name][[a]] and 
P[name][[b]].";

LMeet::usage = "LMeet[name] is the meet table for P[name],
i.e., LMeet[name][[a,b]] is the meet of P[name][[a]] and 
P[name][[b]].";

JoinSubLatticeQ::usage = "JoinSubLatticeQ[name,S] returns 
True if the elements in P[name] indexed by S form a join
sublattice.";

MeetSubLatticeQ::usage = "MeetSubLatticeQ[name,S] returns 
True if the elements in P[name] indexed by S form a meet
sublattice.";

SubLatticeQ::usage = "SubLatticeQ[name,S] returns True if 
the elements in P[name] indexed by S form a sublattice.";

JI::usage = "JI[name] is a list of the join irreducibles of
P[name], assumed to be a lattice.";

MI::usage = "MI[name] is a list of the meet irreducibles of
P[name], assumed to be a lattice.";

ZetaJI::usage = "ZetaJI[name] is the zeta-matrix of the poset
of join-irreducibles of P[name].";

ZetaMI::usage = "ZetaMI[name] is the zeta-matrix of the poset
of meet-irreducibles of P[name].";

JICoverRelations::usage = "JICoverRelations[name] is the list 
of coverings in the poset of join irreducibles of P[name].";

MICoverRelations::usage = "MICoverRelations[name] is the list 
of coverings in the poset of meet irreducibles of P[name]."; 

LinearExtensions::usage="LinearExtensions[name] is a list
of linear extensions of P[name]. (Can take a long time!)";

ToRelation::usage = "ToRelation[matrix] returns the list of
related pairs defined by a 0-1 matrix.";

ToMatrix::usage = "ToMatrix[relation,n] returns the incidence 
matrix (of order n) corresponding to a set of related pairs.";

G::usage= "G[name,q] is the rational generating 
function (in q) whose nth coefficient equals the number of 
P[name]-partitions of n.";

GBar::usage= "GBar[name,q] is the rational generating 
function (in q) whose nth coefficient equals the number of 
strict P[name]-partitions of n.";

Omega::usage= "Omega[name,n] is the order polynomial
of P[name], i.e., the function (of n) which counts order-
preserving maps from P[name] into an n-element chain.";

OmegaBar::usage= "OmegaBar[name,n] is the strict order 
polynomial of P[name], i.e., the function (of n) which counts 
strict order-preserving maps from P[name] into an n-element 
chain";

OmegaGF::usage= "OmegaGF[name,q] is the rational generating
function (of q) whose coefficients are defined by Omega[name,n].";

OmegaBarGF::usage= "OmegaBarGF[name,q] is the rational 
generating function (of q) whose coefficients are defined by 
OmegaBar[name,n].";

TopSort::usage = "TopSort[name] sorts and renumbers the 
elements of P[name] according to height, remaking Rank[name] 
and CoverRelations[name].";

TClosure::usage = "TClosure[relation,n] computes the transitive
closure and transitive reduction of the relation defined by 
{relation,n}. Returns a pair {matrix,relation}, where <matrix> is 
the incidence matrix of the closure and <relation> is a list of
irredundant pairs.";

IsomorphicQ::usage="IsomorphicQ[name1,name2] returns False
if posets P[name1] and P[name2] fail certain isomorphism tests,
and 'Probably' if all tests are passed.";

Locate::usage = "Locate[name,element] returns
the index in P[name] of the given element.";

LocateSet::usage = "LocateSet[name,set] returns
a list of indices corresponding to the positions in P[name] 
of the elements in <set>.";

Verbose::usage = "Verbose=True (default) causes messages to be 
displayed indicating the progress of computations.";

Version::usage= "Returns current version number.";

ListCommands::usage = "ListCommands[] lists all of the commands
defined by the Poset package.  ListCommands[<string>] lists all of 
the commands that contain  <string>.";

NumPosets::usage = "NumPosets[n] returns the number of 
non-isomorphic posets of size n for any n <= 7. This command
assumes that the file allposets1-7 has been read in.";

NumLattices::usage = "NumLattices[n] returns the number of
non-isomorphic lattices of size n for any n <= 9. This command
assumes that the file alllattices1-9 has been read in.";

Zap::usage = "Zap[name] removes all remembered all values
of objects associated with <name>.";

(* Usage:: Operations for combining posets *)

CartesianProduct::usage= "CartesianProduct[name1,name2]
defines the Cartesian product of P[name1] and P[name2]. Output 
is in {relation,card} form, and Build should be applied next.
Input may be in any of the three standard forms.";

DisjointSum::usage= "DisjointSum[name1,name2]
defines the disjoint sum of P[name1] and P[name2]. Output is
in {relation,card} form, and Build should be applied next.
Input may be in any of the three standard forms.";

OrdinalSum::usage= "OrdinalSum[name1,name2]
defines the ordinal sum of P[name1] and P[name2]. Output is
in {relation,card} form, and Build should be applied next.
Input may be in any of the three standard forms.";

SubPoset::usage= "SubPoset[name,points]
defines the subposet of P[name] induced by <points>. Output 
is in {relation,card} form, and Build should be applied next. 
Input may be in any of the three standard forms. <points> may 
be either a list of indices or a Boolean function which defines 
those indices.";

RestoreSubPosetLabels::usage=
"RestoreSubPosetLabels[name,points,subposetname] restores 
original labels in P[name] to subposet";

IntervalP::usage= "IntervalP[name,bottom,top] defines the interval 
subposet [bottom,top] in P[name]. The two endpoints should be 
labels, not indices. Output is a list of indices. Apply SubPoset 
next, then Build.";

JoinSubLattice::usage = "JoinSubLattice[name,S] returns the join 
sublattice of P[name] generated by S. Both input and output
are sets of indices. Apply SubPoset next, then Build.";

MeetSubLattice::usage = "MeetSubLattice[name,S] returns the meet 
sublattice of P[name] generated by S. Both input and output
are sets of indices. Apply SubPoset next, then Build.";

SubLattice::usage = "SubLattice[name,S] returns the  
sublattice of P[name] generated by S. Both input and output
are sets of indices. Apply SubPoset next, then Build.";

Dual::usage= "Dual[name] returns the dual of P[name]. Output is
in {relation,card} form, and Build should be applied next. Input 
may be in any of the three standard forms.";

RestoreDualPosetLabels::usage=
"RestoreDualPosetLabels[name,dualposetname] restores original 
labels in P[name] to dualposet";

CP::usage= "Abbreviation for CartesianProduct.";

DS::usage="Abbreviation for DisjointSum.";

OS::usage="Abbreviation for OrdinalSum.";

SP::usage="Abbreviation for SubPoset.";

JSL::usage="Abbreviation for JoinSubLattice.";

MSL::usage="Abbreviation for MeetSubLattice.";

SL::usage="Abbreviation for SubLattice.";

(* Usage:: Miscellaneous combinatorics *)

InversionCode::usage="InversionCode[sigma] returns the 
inversion coding of sigma (ith element = # inversions (i>x) 
in sigma)."; 

ToPermutation::usage="ToPermutation[alpha] returns the
permutation with inversion code alpha.";

TensorProduct::usage= "TensorProduct[M,N] computes the
tensor product of two matrices M and N.";

InversePermutation::usage= "InversePermutation[p] gives the
inverse permutation of p.";

DES::usage= "DES[sigma] returns the number of descents of 
permutation sigma";

ASC::usage= "ASC[sigma] returns the number of ascents of 
permutation sigma";

MAJ::usage= "MAJ[sigma] returns the major index of
permutation sigma";

AMAJ::usage= "AMAJ[sigma] returns the major (ascent) 
index of permutation sigma";

GFToPoly::usage="GFToPoly[w(q),q,d,n] returns the polynomial 
coefficient f(n) of q^n in the rational generating 
function w(q)/(1-q)^(d+1)";

PolyToGF::usage="PolyToGF[f,x,q] returns the rational 
generating function w(q)/(1-q)^(d+1) whose coefficients are 
given by f(x).";

ZeroDiffs::usage="ZeroDiffs[sequence] returns a list whose kth
element is the (k-1)st difference at 0 of the sequence.";

InfiniteProduct::usage="Computes infinite product representation
of an ordinary generating function.";

GenFun::usage="GenFun[L,q] makes a generating function where 
the coefficient of q^k is the number of times k appears in 
list L.";

                 Begin["`Inside`"];


(* Standard poset definitions *)

Subsets[n_]:= {Substring,{Array[#&,n]},n};
Substring[w_]:= Substring[w]= 
     Union[Table[Drop[w,{j,j}],{j,Length[w]}]]

Subwords[w_]:= {Substring,{w},Length[w]};
Substring[w_]:= Substring[w]= 
     Union[Table[Drop[w,{j,j}],{j,Length[w]}]]

Div[n_] := {Divides,{n},n}  ;                    
Divides[n_] := 
  (n/#)& /@ Select[Divisors[n],PrimeQ];


SetP[n_]:= {SetRefine, {Range[n]},n-1} ;                     
SetRefine[w_] := Block[
   {outlist={},i,j},
    Do[If[w[[i]]==i,
         Do[If[w[[j]]==j,
             outlist = Append[outlist,
             MapAt[w[[i]]&,w,Position[w,w[[j]]] ] ]              
          ],{j,i+1,Length[w]}
         ]
       ],{i,Length[w]}];
    outlist      
  ]                       

NonXP[n_]:= {NonXPRefine, {Array[#&,n]},n-1}  ;                     
NonXPRefine[w_] := Block[
  {outlist={},i,j,trial},
   Do[If[w[[i]]==i,
       Do[If[w[[j]]==j,
        trial = MapAt[w[[i]]&,w,Position[w,w[[j]]] ];
        If[NonCrossing[trial],
             outlist = Append[outlist,trial]
        ]
        ],{j,i+1,Length[w]}
       ]
      ],{i,Length[w]}];
   outlist      
];                      
NonCrossing[{___,x_,___,y_,___,x_,___,y_,___}] := 
					False /; x!=y
NonCrossing[{___}] := True

WeakS[n_] := {WeakOrder,{Array[#&,n]},Binomial[n,2]}  ;                    
WeakOrder[w_] := Block[
  {outlist={},i,j},
  Do[If[w[[i]]<w[[i+1]],
       AppendTo[outlist,
          Insert[Drop[w,{i}],w[[i]],i+1]]],
  {i,Length[w]-1}];
  outlist      
]

StrongS[n_] := {StrongOrder,{Array[#&,n]},Binomial[n,2]}  ;                    
StrongOrder[w_] := Block[
     {outlist={},i,j,trial},
     Do[If[w[[i]]<w[[j]],trial=w;
         trial[[i]]=w[[j]];
         trial[[j]]=w[[i]];
         If[okstrong[w,i,j],AppendTo[outlist,trial]]],
     {i,Length[w]-1},{j,i+1,Length[w]}];
     outlist      
];
okstrong[w_,i_,j_] := Block[
  {between=Take[w,{i+1,j-1}]},
   (Intersection[between,Range[w[[i]],w[[j]]]]=={})
   ];

YoungsLattice[partition_List]:= {Young, {partition}, 
     Apply[Plus,partition]}                      
Young[w_] := Block[
  {outlist={},i,j},
   Do[If[w[[i]]> w[[i+1]],
      AppendTo[outlist,MapAt[w[[i]]-1&,w,{{i}}]]
      ],
   {i,Length[w]-1}];
   If[Last[w]>0,
     AppendTo[outlist,MapAt[Last[w]-1&,w,{{Length[w]}}]]
   ];        
  outlist      
]                     

MajP[n_]:={Majorize,{{n}},n^2}
Clear[Lineword];
Majorize[w_] :=  Block[
   {i,j,outlist,temp,pp,q},
   outlist = {};
   Do[temp = Append[w,0];
    If[temp[[j+1]]<temp[[j]],
       If[temp[[j+1]]<temp[[j]]-1,
           temp[[j]]=temp[[j]]-1;
           temp[[j+1]]=temp[[j+1]]+1;
           temp=If[Last[temp]==0,
                Drop[temp,-1],temp];
           AppendTo[outlist,temp],
         (*else*)
           pp = Position[temp,temp[[j]]-2];
           If[pp!={},
             q=pp[[1,1]];
             temp[[j]]=temp[[j]]-1;
             temp[[q]]=temp[[q]]+1;
             temp=If[Last[temp]==0,
                  Drop[temp,-1],temp
             ];
             AppendTo[outlist,temp]
           ]   
      ]
    ],
  {j,Length[w]}];
  Union[outlist]      
 ]                

Chain[n_] := {Table[{i,i+1},{i,n-1}],n}  ;                 

Antichain[n_] := {{},n}  ;                 

ZigZag[n_] := 
 {Table[If[OddQ[i],{i,i+1},{i+1,i}],{i,n-1}],n}  ;                 

RandomP[n_,p_] := 
Table[If[i<j,If[Random[]<p,1,0],0],
    {i,n},{j,n}];               

NumP[n_] := {NumberRefine,{{n}},n-1}  ;                 
NumberRefine[w_] := Block[
  {outlist={},i,j,k,delete,result},
   Do[If[w[[i]]> 1,
         k=w[[i]];
         delete = Drop[w,{i,i}];
         For[j=1,j<=k-j,j++,
            result=Sort[Join[delete,{j,k-j}] ];
            outlist = Append[outlist,result]
          ]
       ],
   {i,Length[w]}];        
Union[outlist]       
]           

ContractionLattice[edges_,n_] := {contract[edges,#]&,{Range[n]},n-1};
contract[edges_,g_] := Block[{outlist={},i,j,k,low,high},
  Do[
  {low,high} = Sort[{g[[edges[[k,1]]]],g[[edges[[k,2]]]]}]; 
		  outlist = Append[outlist,
		     MapAt[low&,g,Position[g,high]]],
   {k,Length[edges]}
 ];
 outlist
]

PosetP[name_] := 
 Block[{},
  (*internal function: can two chains merge? *)
  canmergeQ[namea_,c1_,c2_] := 
   Block[{i,j,check=True},
     For[i=1, check && i<=Length[c1], i++, 
      For[j=1, check && j<=Length[c2], j++, 
          check=MemberQ[Down[namea][[Max[{c1[[i]],c2[[j]]}]]],
                         Min[{c1[[i]],c2[[j]]}]]
         ]
        ];
     check
    ];                      
  (*covering function*)
  PosRefine[w_] := 
   Block[{outlist={},check=True,i,j},
    Do[If[w[[i]]==i,
          Do[If[w[[j]]==j && canmergeQ[name,Flatten[Position[w,i]],
                                   Flatten[Position[w,j]]],
                AppendTo[outlist,MapAt[w[[i]]&,w,Position[w,w[[j]]]]]
               ],{j,i+1,Length[w]}
            ]
         ],{i,Length[w]}
      ];
     outlist      
    ];
  Return[{PosRefine,{Range[Card[name]]},Card[name]-1}]
 ]                  

JofP[namea_] := Block[{in},     
     If[Head[namea]===Symbol,in=Transpose[ZetaP[namea]];
          up=Up[namea]
     ];
     If[Head[namea]===List,
          If[MatrixQ[namea],
               in=Transpose[TClosure[namea][[1]]];
               up=Map[Flatten[Position[in[[i]],
                     1]],Range[Length[in]]],
               in=Transpose[TClosure[ToMatrix[namea[[1]],
                    namea[[2]] ]][[1]]];
               up=Map[Flatten[Position[in[[i]],
                     1]],Range[Length[in]]]
          ]
     ];
     PosetElements = Table[j, {j, Length[in]}]; 
     NextIdeals[w_] := Block[{i, PrincipalOrderIdeal, 
         NewPosetElements, output},
         PrincipalOrderIdeal[x_] := Complement[up[[x]], w]; 
         NewPosetElements = Complement[PosetElements, w]; 
         output = {}; 
         Do[  If[PrincipalOrderIdeal[
                   NewPosetElements[[i]]] ==
         	           {NewPosetElements[[i]]}, 
                   AppendTo[output,
                       Union[Append[w,
                       NewPosetElements[[i]] ]]
                   ]
              ], 
              {i, Length[NewPosetElements]}
          ]; 
          output
      ]; 

      Return[{NextIdeals,{{}},
          Card[namea]}]	
]	

(* Constructing the fundamental objects *)

Build[{scan_,minlist_,hbound_Integer},name_:automatic] :=
   Block[{plist=minlist,links,scanned,coversi,i,r,k,pos,
        rank},
        If[name==automatic,name=Unique["Poset"]];
        Zap[name];
        Print["Building poset ",name,"  ..."];
        links={};
        rank = Map[0&,plist];
        scanned = Map[False&,plist];
        i=1;r=0;
        If[Verbose,Print["working on iteration 0"]];
        While[(i<=Length[plist])&&
                  (!scanned[[i]])&&(rank[[i]]<hbound),        
             If[Verbose && (r!=rank[[-1]]), 
               Print["number generated at iteration ",r," = "
                  ,Length[Position[rank,r]]];
               r=rank[[-1]];
               Print["working on iteration ",r]];
             coversi = scan[plist[[i]]];
             scanned[[i]] = true;
             Do[
                  If[  !MemberQ[plist,coversi[[k]]],
                       AppendTo[plist,coversi[[k]]];
                       AppendTo[rank,rank[[i]]+1];
                       AppendTo[scanned,False];
                       AppendTo[links,{i,Length[plist]}],
                   (*else*)
                       AppendTo[links,
                            {i,Position[plist,
                             coversi[[k]]][[1,1]]}
                       ]
                   ],
                   {k,Length[coversi]}
               ];
               i=i+1;
        ];
     Rank[name]=rank;
     CoverRelations[name]=links;
     P[name]=plist;
     If[!Apply[And,Map[(#[[1]] < #[[2]])&, links]], 
         TSort1[name]];
     If[Verbose,Print["number generated at iteration ",rank[[-1]]," = ",
       Length[Position[rank,rank[[-1]]]]]];
     Print["Done"]
];

Build[M_List,name_:automatic]:=
Block[{n=Length[M]},
     If[name==automatic,name=Unique["Poset"]];
     Print["Building poset ",name,"  ..."];
     Zap[name];
     Card[name]=n;
     CoverRelations[name]=ToRelation[M];
     TSort2[name]; 
     {ZetaP[name],CoverRelations[name]}=TClosure[CoverRelations[name],n];
     Print["Done"]
]

Build[{links_,card_Integer},name_:automatic]:=
Block[{n=card},
     If[name==automatic,name=Unique["Poset"]];
     Print["Building poset ",name,"  ..."];
     Zap[name];
     Card[name]=n;
     CoverRelations[name]=links;
     TSort2[name];
     {ZetaP[name],CoverRelations[name]}=TClosure[CoverRelations[name],n];
     Print["Done"]
]

ToRelation[M_] := Position[M,1];

ToMatrix[L_,n_] := Block[{out=Array[0&,{n,n}]},
   Scan[(out[[#[[1]],#[[2]]]]=1)&,L];out];

LongestChain[links_,n_] := Block[
{out=Array[{}&,n],in=hts=Array[0&,n],min={},dlinks,
labeled=0,h=0,k},
update[k_]:= Scan[(in[[#]]--)&,out[[k]]];
dlinks=Select[links,First[#]!=Last[#]&];
Scan[in[[Last[#]]]++&,dlinks];
Scan[AppendTo[out[[First[#]]],Last[#]]&,dlinks];
While[labeled<n,
  Scan[(in[[#]]=-1)&,min];
  min=Flatten[Position[in,0]];
  Scan[
    (hts[[#]]=h;update[#]
     )&,min];
  labeled=labeled+Length[min];
  h=h+1;
  ];
hts
]

TClosure[links_List,n_]:= Block[{covers,a,k,j,
     outmatrix,reducedlinks,redundant={}},
     If[Verbose,
Print["Computing transitive closure, transitive reduction ..."];
        ];
     outmatrix=IdentityMatrix[n];
     covers=Array[{}&,n];
     Scan[(AppendTo[covers[[First[#]]],Last[#]])&,
           links];
     covers=Map[Union,covers];      
     
     doit[x_]:= Block[{},
       Scan[
         (If[outmatrix[[a,#]]==1,
             AppendTo[redundant,{a,#}],
             outmatrix[[a,#]]=1];)&,
               Flatten[Position[outmatrix[[x]],1]] 
       ]
     ];
     For[a=n,a>=1,a--,
          Scan[(doit[#])&,covers[[a]] ]
     ];
     reducedlinks=Complement[links,redundant];
     {outmatrix,reducedlinks}
]        

TopSort[links_,n_] := Block[{height},
     height=LongestChain[links,n];
     Return[Map[#[[2]]&,
          Sort[Table[{height[[i]],i},{i,n}]]]]  
    ];       

(* Internal procedure used by first Build command *)
TSort1[name_] := Block[{height,n=Card[name],permute,labelp},
     height=LongestChain[CoverRelations[name],n];
     labelp=Map[#[[2]]&,
          Sort[Table[{height[[i]],i},{i,n}]]];  
     Rank[name]=Sort[height];
     permute=Sort[Table[{labelp[[i]],i},{i,n}]];
     CoverRelations[name]=
          Map[{permute[[First[#],2]],permute[[Last[#],2]]}&,
               CoverRelations[name]
     ];
     P[name]=Map[P[name][[#]]&,labelp];
];       

(* Internal procedure used by last two Build Commands *)
TSort2[name_] := Block[{height,n=Card[name],permute},
     If[Verbose,
     Print["Renumbering poset ",name];
     Print["Height := length of longest chain."];
     Print["Making P, Rank, CoverRelations ... "]];
     height=LongestChain[CoverRelations[name],n];
     P[name]=Map[#[[2]]&,
          Sort[Table[{height[[i]],i},{i,n}]]];  
     Rank[name]=Sort[height];
     permute=Sort[Table[{P[name][[i]],i},{i,n}]];
     CoverRelations[name]=
          Map[{permute[[First[#],2]],permute[[Last[#],2]]}&,
               CoverRelations[name]
     ];
];       

PGraded[name_] := PGraded[name] = Block[{},
    If[Verbose,Print["Computing PGraded ... "]];
    Table[
       P[name][[Flatten[Position[Rank[name],k]]]],
       {k,0,Max[Rank[name]]}] 
    ]
  
Card[name_] := Card[name] = Length[P[name]];

NK[name_] := NK[name] = Map[Length,PGraded[name]];       

H[name_] := H[name] = Max[Rank[name]];

RGF[name_]:=RGF[name]=Apply[Plus,Global`q^Rank[name]];
    


Relabel[name_,f_:Automatic] := Block[{},
  If[Verbose,Print["Relabeling elements of ",name,"."]];
  (*If[Verbose,Print["Saving old labels as OldP[",name,"]."]];*)
  (*OldP[name]=P[name];*)
  If[f===Automatic,
    P[name]=Range[Card[name]],
  (*else*)
    P[name]=Map[f,P[name]]];
  ];     

Compact[lab_] := Apply[StringJoin,Map[ToString,lab]];

SortByRanks[name_] := Block[{newp},
  If[Verbose,Print["Remaking P, CoverRelations, PGraded ... "]];
  PGraded[name] = Map[Sort,PGraded[name]];
  newp = Flatten[PGraded[name],1];
  CoverRelations[name] =Map[
    Position[newp,P[name][[#]]][[1,1]]&,CoverRelations[name],{2}];
  P[name] = newp;
  ]
  

(* Testing for rank function *)

RankedQ[name_]:= RankedQ[name] = Block[{vars,eqns,maxels,
	moreeqns,S,k,l,m,r,arank,rank=Rank[name],
	links=CoverRelations[name]},
	If[Verbose,Print["Testing if poset is ranked ..."]];
	maxrank=Max[rank];
	maxels=Flatten[Position[rank,maxrank]];
	vars=Table[r[l],{l,Card[name]}];
	eqns=Table[r[links[[m,1]]]==r[links[[m,2]]]-1,
	      {m,Length[links]}
	];
	moreeqns=Table[r[maxels[[k]]]==maxrank,{k,Length[maxels]}];
	eqns=Flatten[AppendTo[eqns,moreeqns]];
	S=Solve[eqns,vars];
	If[S=={},
	    Print["The poset is not ranked."];
	    If[Verbose,Print["Remaking Rank,PGraded,NK..."]];
		Off[Unset::norep];
		PGraded[name]=.;NK[name]=.;Rank[name]=.;H[name]=.;
		Rank[name]=LongestChain[CoverRelations[name],Card[name]];
		Return[False],
      (*otherwise*)
		If[Flatten[vars/.S]==rank,
		     Print["The poset is strongly ranked."];
		     Off[Unset::norep],
		     Print["The poset is weakly ranked."];
		     If[Verbose,Print["Remaking Rank, PGraded, NK."]];
		     Off[Unset::norep];
		     PGraded[name]=.;NK[name]=.;Rank[name]=.;RGF[name]=.;
		     H[name]=.;
		     PGraded[name]=Table[Flatten[Position[Flatten[vars/.S],
		          arank-1]],{arank,maxrank+1}];
	         NK[name]=Map[Length,PGraded[name]];
		     Rank[name]=Flatten[vars/.S];
	    ];
	    On[Unset::norep];
	    Return[True]
	]
	(*]*)
]


(* Hasse Diagram *)

Options[Diagram]={ShowLabels -> False, Thinness -> Automatic,
     Jiggled->Automatic}; 

Diagram[name_Symbol,y___List,opts___Rule]:=Block[{spoints,slines,
     r,i,covline,xy,dotobj,lineobj,labelobj,linethick,pointthick,
	 p=P[name],pgraded=PGraded[name],
	 nk=NK[name],cardp=Card[name],
	 links=CoverRelations[name],showlabels,thinner,jiggled},
     
     If[{y}=={},
         spoints={};slines={}, (* else *)
         If[Length[{y}]==1,
             If[!MatrixQ[y],spoints={y}[[1]]; slines={},
                            spoints={};  slines={y}[[1]]], (*else*)
             spoints={y}[[1]];slines={y}[[2]]
          ]
     ];
     
     showlabels = ShowLabels/.{opts} /.Options[Diagram];
	 thinner = Thinness/.{opts} /.Options[Diagram];
	   If[thinner===Automatic,thinner=1];
	 jiggled = Jiggled/.{opts} /.Options[Diagram];
	   If[jiggled===Automatic,jiggled=0];
	   
	 If[Verbose,Print["Computing XY, Dots, Lines, Labels ..."]];
	 
     nk=NK[name];
  	 xy=Flatten[Table[
      	  Table[{i/(nk[[r]]+1)+
      	     (Random[]-.5)/(20 Max[nk])*jiggled,r},
      	       {i,nk[[r]]}
      	  ],
     {r,H[name]+1}],1
     ];
 	 
     dotobj=Map[Point,xy];
	 covline[{$p1_,$p2_}]:= 
	      Line[{xy[[$p1]],xy[[$p2]]}];
	 lineobj=Map[covline,links];
	 labelobj=Table[Text[p[[$b]],xy[[$b]],{-1.6,-1}],
		         {$b,cardp}
     ];
     
     specialps = Map[Point[xy[[#]]]&,spoints];
     specialels = Map[Line[{xy[[#[[1]] ]],xy[[#[[2]] ]]}]&,
          slines];
          
     linethick = (.0028)*thinner;
     pointthick =  (.02)*thinner;
         
	 If[showlabels,
	      Show[Graphics[{Thickness[linethick],PointSize[pointthick],
                            lineobj,
                            dotobj,
                            PointSize[pointthick*2],
                            RGBColor[0,0,1],
                            specialps,
                            Thickness[linethick*4],
                            RGBColor[.7,0,.5],
                            specialels,
                            labelobj}
                    ],
                    PlotRange->{{0,1},{0,H[name]+2}},
                    Frame->True,FrameTicks->None,
                    ColorOutput->Automatic
          ],
          Show[Graphics[{Thickness[linethick],
				      PointSize[pointthick],
                      lineobj,
                      dotobj,
                      PointSize[pointthick*(2)],
                      RGBColor[0,0,1],
                      specialps,
                      Thickness[linethick*(4)],
                      RGBColor[.7,0,.5],
                      specialels}
                ],
                PlotRange->{{0,1},{0,H[name]+2}},
                Frame->True,FrameTicks->None,
                ColorOutput->Automatic
          ]
     ]
]
     

(* ZetaP, Mu, Up, Down, Covers *)

ZetaP[name_]:= ZetaP[name] = Block[{covers,a,k,j,zetap,
     cardp,links},
     If[Verbose,Print["Computing ZetaP ..."]];
     cardp=Card[name];
     links=CoverRelations[name];
     zetap=IdentityMatrix[cardp];
     covers=Array[{}&,cardp];
     Scan[(AppendTo[covers[[First[#]]],Last[#]])&,
           links];
     
     doit[x_]:= Block[{},
          Scan[(zetap[[a,#]]=1)&,
               Flatten[Position[zetap[[x]],1]] 
          ]
     ];
     For[a=cardp,a>=1,a--,
          Scan[(doit[#])&,covers[[a]] ]
     ];
     zetap
]        

Up[name_]:= Up[name] = Block[{j},
     If[Verbose,Print["Making Up ..."]];
     Table[Flatten[Position[ZetaP[name][[j]],1]],
          {j,Card[name]}
     ]
]

Down[name_] := Down[name] = Block[{j,transpzeta=
     Transpose[ZetaP[name]]},
     If[Verbose,Print["Making Down ..."]];
     Table[Flatten[Position[transpzeta[[j]],1]],
          {j,Card[name]}
     ]
]

Covers[name_] := Covers[name] = Block[{out},
  out=Array[{}&,Card[name]];
  Scan[(AppendTo[out[[First[#]]],Last[#]])&,
         CoverRelations[name]];
  out];

CoCovers[name_] := CoCovers[name] = Block[{out},
  out=Array[{}&,Card[name]];
  Scan[(AppendTo[out[[Last[#]]],First[#]])&,
         CoverRelations[name]];
  out];

UpDegree[poset_]:=Table[Count[CoverRelations[poset],{i,_}],
     {i,Card[poset]}]

DownDegree[poset_]:=Table[Count[CoverRelations[poset],{_,i}],
     {i,Card[poset]}]

Mu[name_] := Mu[name] = Inverse[ZetaP[name]]

(* Counting chains, Zeta polynomial, Order Polynomial *)

MaximalChainsDown[name_] := MaximalChainsDown[name] =
Block[{a,num=Table[0,{Card[name]}], minlist, card=Card[name],
       links=CoverRelations[name]},
minlist=Complement[Range[card],Map[#[[2]]&,links]];
Do[If[MemberQ[minlist,a],num[[a]]=1],
  {a,card}];
Do[num[[links[[a,2]]]]=num[[links[[a,2]]]]+num[[links[[a,1]]]],
  {a,Length[links]}];
num
]

MaximalChainsUp[name_] := MaximalChainsUp[name] =
Block[{a,num=Table[0,{Card[name]}],maxlist, card=Card[name],
       links=CoverRelations[name]},
maxlist=Complement[Range[card],Map[#[[1]]&,links]];
Do[If[MemberQ[maxlist,a],num[[a]]=1],
  {a,card}];
For[a=Length[links],a>=1,a--,
num[[links[[a,1]]]]=num[[links[[a,1]]]]+num[[links[[a,2]]]] ];
num
]

ChainsBetweenGF[name_,a_,b_,q_:Global`q] := ChainsBetweenGF[name,a,b,q]=
 Block[{x,above=Up[name][[a]],chains,interval,intsize},
  between[namea_,top_]:= 
   Intersection[above,Down[namea][[top]]];
  interval=between[name,b];
  intsize=Length[interval];
  chains=Map[0&,interval];
  Do[chains[[x]]=q + q *
   Apply[Plus,Map[chains[[Position[interval,#][[1,1]]]]&,
      Select[interval,
             MemberQ[between[name,interval[[x]]],#]&]
            ]
        ],
  {x,2,intsize}];
  If[intsize==0,Return[0],
   (* else *)
   If[intsize==1,Return[1],
    Return[chains[[intsize]] // Expand]
     ]
    ]
 ];


ZetaPoly[name_,a_,b_,n_:Global`n] := ZetaPoly[name,a,b,n] =
 Block[{lengths},
  lengths=CoefficientList[
         Expand[ChainsBetweenGF[name,a,b,q]],q];
  Do[lengths[[i]]={lengths[[i]],i-1},{i,Length[lengths]}];
  Apply[Plus,
   Map[(Binomial[n,#[[2]]] #[[1]])&, lengths]]//Expand
  ]

ZetaPoly[name_,n_:Global`n] := ZetaPoly[name,n] =
 Block[{chains=Array[0&,Card[name]],lengths,numchainsgf},
  Do[chains[[x]]=Expand[q+
    q Sum[chains[[Down[name][[x,i]]]],
        {i,Length[Down[name][[x]]]-1}]],
    {x,Card[name]}];
  numchainsgf=Apply[Plus,Map[Expand[(1/q) #]&,chains]];
  lengths=CoefficientList[numchainsgf,q];
  Do[lengths[[i]]={lengths[[i]],i-1},{i,Length[lengths]}];
  Apply[Plus,
   Map[(Binomial[n-2,#[[2]]] #[[1]])&, lengths]]//Expand
  ];

DES[s_] := Sum[If[s[[i]]>s[[i+1]],1,0],{i,Length[s]-1}];

ASC[s_] := Sum[If[s[[i]]<s[[i+1]],1,0],{i,Length[s]-1}];

MAJ[s_] := Sum[If[s[[i]]>s[[i+1]],i,0],{i,Length[s]-1}];

AMAJ[s_] := Sum[If[s[[i]]<s[[i+1]],i,0],{i,Length[s]-1}];

LinearExtensions[name_] := LinearExtensions[name] = 
  Block[{out={{}},k,scan},
    If[Verbose,
    Print["Computing linear extensions ...."];
    Print["(This can take a while.)"];];
  scan[L_,k_] := Block[{max,i},
    max=1;
    Do[If[MemberQ[Down[name][[k]],L[[i]]],max=i+1],{i,k-1}];
    Table[Insert[L,k,i],{i,max,Length[L]+1}]
    ];
  Do[out=Map[scan[#,k]&,out];
     out=Flatten[out,1], {k,1,Card[name]}];
  out
  ];

G[poset_,q_] := Block[{i},
   (Plus @@ (q^Map[MAJ,LinearExtensions[poset]]))/
   Product[(1-q^i),{i,Card[poset]}]
   ];

GBar[poset_,q_] := Block[{i},
   (Plus @@ (q^Map[AMAJ,LinearExtensions[poset]]))/
   Product[(1-q^i),{i,Card[poset]}]
   ];

Omega[poset_,n_] := Block[{w},
  w=(Plus @@ (q^Map[DES,LinearExtensions[poset]]));
  GFToPoly[q w,q,Card[poset],n]
  ];

OmegaBar[poset_,n_] := Block[{w},
  w=(Plus @@ (q^Map[ASC,LinearExtensions[poset]]));
  GFToPoly[q w,q,Card[poset],n]
  ];

OmegaGF[poset_,q_] := Block[{w},
  (Plus @@ (q^(Map[DES,LinearExtensions[poset]]+1)))/
   (1-q)^(Card[poset]+1)
  ];

OmegaBarGF[poset_,q_] := Block[{w},
  (Plus @@ (q^(Map[ASC,LinearExtensions[poset]]+1)))/
   (1-q)^(Card[poset]+1)
  ];

(* MaxAntichain and Dilworth covering *)

MaxAntichain[name_] := MaxAntichain[name] = Block[
   {m=n=Card[name],Amatch,Bmatch,Blabel,
         queue,a,b,moretodo=True,success,scan,augment},
   If[Verbose,Print["Computing MaxAntichain and DilworthCover ..."]];
   Amatch=Array[0&,m]; Bmatch=Array[0&,n];
   
   scan[a_] := Block[{nbrs=Drop[Up[name][[a]],1]},
    success=False;
    Scan[
     (If[Bmatch[[#]]==0,success=True;b=#;Return[],
      If[Blabel[[#]]==0,
             Blabel[[#]]=a;
             PrependTo[queue,Bmatch[[#]]]
             ]] )&,nbrs];
    ];
   
    augment := Block[{$a=a,$b=b},
      While[Amatch[[$a]]!=0,
        $b=Amatch[[$a]];
        Bmatch[[b]]=$a;Amatch[[$a]]=b;
        b=$b;$a=Blabel[[b]];
        ];
      Bmatch[[$b]]=$a;Amatch[[$a]]=b;
      ];
     
  While[moretodo,
    Blabel=Array[0&,n];      
    queue=Complement[Range[m],Bmatch];   
    success=False;
    While[queue!={}, 
      a=Last[queue];queue=Drop[queue,-1];
      scan[a];
      If[success,augment;queue={}];
       ];
    moretodo=success;
    ]; 
  DilworthCover[name]=
    Cases[Table[{i,Amatch[[i]]},{i,m}],{_,_?Positive}] ;
  Complement[Range[m],Bmatch[[Flatten[Position[Blabel,0]]]],
     Flatten[Position[Blabel,_?Positive]]]   
  ]   

ff[{u___,{a__,x_},v___,{x_,b__},w___}] := 
       {u,{a,x,b},v,w};
ff[{u___,{x_,b__},v___,{a__,x_},w___}] := 
       {u,{a,x,b},v,w};
ff[{u___}] := {u};
Fuse[L_] := FixedPoint[ff,L];

(* Characteristic Polynomial *)

CharPoly[poset_,q_] := CharPoly[poset,q] = Block[{k},
	Sum[Mu[poset][[1,k]]*
	q^(Rank[poset][[Length[P[poset]]]]-
	Rank[poset][[k]]),
	{k,1,Length[P[poset]]}]]

(* Lattice operations *)

LatticeQ[name_]:= LatticeQ[name] =
    Block[{$a,$b,up=Up[name],ubounds},
    If[Verbose,Print["Testing if P is a lattice ..."]];
    If[up[[1]]!=Range[Card[name]],
       Print["Not a lattice: no 0"];Return[False]];
    For[$a=1,$a<=Card[name],$a++,
    For[$b=1,$b<$a,$b++,
    ubounds = Intersection[up[[$a]],up[[$b]]];
    If[ubounds=={},
      Print["Not a lattice: P[",$a,"] and P[",$b,"] have no upper bounds."];
      Return[False]];
    If[ubounds!=up[[First[ubounds]]],
      Print["Not a lattice: P[",$a,"] and P[",$b,"] have no LUB."];
      Return[False]]
    ]];
    True
    ]    

DistributiveLatticeQ[name_] := 
  (RankedQ[name])&&(LatticeQ[name])&&
    (Length[JI[name]]==Length[MI[name]]==Max[Rank[name]])

Vee[name_,x_,y_] :=  P[name][[
   First[Intersection[
      Up[name][[Locate[name,x]]],Up[name][[Locate[name,y]]]
      ]]
      ]];

Wedge[name_,x_,y_] :=  P[name][[
   Last[Intersection[
      Down[name][[Locate[name,x]]],Down[name][[Locate[name,y]]]
      ]]
      ]];

LJoin[name_]:= LJoin[name]=
    Block[{a,b,cup,up=Up[name]},
    If[!LatticeQ[name],Print["Not a lattice."],
    If[Verbose, Print["Making LJoin ..."]];
    cup[a_,b_] :=  First[Intersection[up[[a]],up[[b]]]];
    Table[cup[a,b],
           {a,Card[name]},{b,Card[name]}]
    ]]
    
    
LMeet[name_]:= LMeet[name]=
    Block[{a,b,cap,down=Down[name]},
    If[!LatticeQ[name],Print["Not a lattice."],
    If[Verbose, Print["Making LMeet ..."]];
    cap[a_,b_] :=  Last[Intersection[down[[a]],down[[b]]]];
    Table[cap[a,b],
           {a,Card[name]},{b,Card[name]}]
    ]]

(* Meet- and join-sublattices; sublattices *)

JoinSubLatticeQ[tag_,points_]:= 
	Complement[LJoin[tag][[points,points]]//Flatten,points]=={}

MeetSubLatticeQ[tag_,points_]:= 
	Complement[LMeet[tag][[points,points]]//Flatten,points]=={}

SubLatticeQ[tag_,points_]:= 
	JoinSubLatticeQ[tag,points]&&MeetSubLatticeQ[tag,points]

(* Join- and meet- irreducibles *)

JI[name_]:=JI[name] = Block[{down=Down[name]},
   ji[x_] := Complement[down[[x]],
      down[[down[[x,Length[down[[x]]]-1]]]]
      ];
   Flatten[Position[
       Map[Length,Map[ji,Range[2,Card[name]]]],1]]+1
   ]

ZetaJI[name_] := ZetaJI[name] =
            ZetaP[name][[JI[name],JI[name]]] 

MI[name_]:=MI[name] = Block[{up=Up[name]},
   mi[x_] := Complement[up[[x]], up[[up[[x,2]]]]];
   Flatten[Position[
       Map[Length,Map[mi,Range[1,Card[name]-1]]],1]]
   ]

ZetaMI[name_] := ZetaMI[name] =
            ZetaP[name][[MI[name],MI[name]]] 

JICoverRelations[name_] := JICoverRelations[name] = Block[
     {StrictJiP,TwoStepJiP,CoversJiP,links,ji=JI[name]},
     StrictJiP =
          ZetaJI[name]-IdentityMatrix[Length[ji]];
     TwoStepJiP = StrictJiP.StrictJiP;
     CoversJiP = StrictJiP-Map[If[#>0,1,0]&,TwoStepJiP,{2}];
     links = Map[ji[[#]]&,Position[CoversJiP,1],{2}]
]

MICoverRelations[name_] := MICoverRelations[name] = Block[
     {StrictMiP,TwoStepMiP,CoversMiP,links,mi=MI[name]},
     StrictMiP =
          ZetaMI[name]-IdentityMatrix[Length[mi]];
     TwoStepMiP = StrictMiP.StrictMiP;
     CoversMiP = StrictMiP-Map[If[#>0,1,0]&,
          TwoStepMiP,{2}
     ];
     links = Map[mi[[#]]&,Position[CoversMiP,1],{2}]
]


(* Operations for combining posets *)

TensorProduct[M_,N_] :=  Map[Flatten,
    Flatten[Transpose[
      Outer[Times,M,N],{1,3,2,4}],1]]

CartesianProduct[poset1_,poset2_] := 
  Block[{mat1,mat2,mat3,card1,card2},
  If[Head[poset1]===Symbol,
     {mat1,card1}={ZetaP[poset1],Card[poset1]},(*else*)
     If[MatrixQ[poset1],
       {mat1,card1}={poset1,Length[poset1]},(*else*)
       {mat1,card1}=
           {ToMatrix[poset1[[1]],poset1[[2]]],poset1[[2]]}
     ]
   ];
  If[Head[poset2]===Symbol,
     {mat2,card2}={ZetaP[poset2],Card[poset2]},(*else*)
     If[MatrixQ[poset1],
       {mat2,card2}={poset2,Length[poset2]},(*else*)
       {mat2,card2}=
           {ToMatrix[poset2[[1]],poset2[[2]]],poset2[[2]]}
     ]
   ];   
  card=card1 card2;
  Do[mat1[[i,i]]=1,{i,card1}];
  Do[mat2[[i,i]]=1,{i,card2}];
  mat3=TensorProduct[mat1,mat2];
  {ToRelation[mat3-IdentityMatrix[card]],card}
  ];

DisjointSum[poset1_,poset2_] := Block[
{card1,links1,card2,links2,links,card},
  If[Head[poset1]===Symbol,
     {links1,card1}={CoverRelations[poset1],Card[poset1]},(*else*)
     If[MatrixQ[poset1],
       {links1,card1}=
            {ToRelation[poset1],Length[poset1]},(*else*)
       {links1,card1}={poset1[[1]],poset1[[2]]}
     ]
   ];
  If[Head[poset2]===Symbol,
     {links2,card2}={CoverRelations[poset2],Card[poset2]},(*else*)
     If[MatrixQ[poset2],
       {links2,card2}=
            {ToRelation[poset2],Length[poset2]},(*else*)
       {links2,card2}={poset2[[1]],poset2[[2]]}
     ]
   ];
  card=card1+card2;
  links=Join[links1,links2+card1];
  {links,card}
  ];

OrdinalSum[poset1_,poset2_] := Block[
{card1,links1,card2,links2,links,card},
  If[Head[poset1]===Symbol,
     {links1,card1}={CoverRelations[poset1],Card[poset1]},(*else*)
     If[MatrixQ[poset1],
       {links1,card1}=
            {ToRelation[poset1],Length[poset1]},(*else*)
       {links1,card1}={poset1[[1]],poset1[[2]]}
     ]
   ];
  If[Head[poset2]===Symbol,
     {links2,card2}={CoverRelations[poset2],Card[poset2]},(*else*)
     If[MatrixQ[poset2],
       {links2,card2}=
            {ToRelation[poset2],Length[poset2]},(*else*)
       {links2,card2}={poset2[[1]],poset2[[2]]}
     ]
   ];
  card=card1+card2;
  links=Join[links1,links2+card1,
       Flatten[Outer[List,
         Range[card1],Range[card2]+card1],1]
       ];
  {links,card}
  ];


Dual[poset_] := Block[
{links1,card,links},
  If[Head[poset]===Symbol,
     {links1,card}={CoverRelations[poset],Card[poset]},(*else*)
     If[MatrixQ[poset],
       {links1,card}=
            {ToRelation[poset],Length[poset]},(*else*)
       {links1,card}={poset[[1]],poset[[2]]}
     ]
   ];
  links=Map[Reverse,Map[(card+1-#)&,links1,{2}]];
  {links,card}
  ];
  
RestoreDualPosetLabels[poset_,dualposet_] := Block[{card=Card[poset]},
  P[dualposet] = Map[P[poset][[card+1 - #]]&,P[dualposet]];
  ];

SubPoset[poset_,points_] := Block[
{mat,links,card,pointlist},
  If[Head[poset]===Symbol,
     {mat,card}={ZetaP[poset],Card[poset]},(*else*)
     If[MatrixQ[poset],
       {mat,card}={poset,Length[poset]},(*else*)
       {mat,card}=
           {ToMatrix[poset[[1]],poset[[2]]],poset[[2]]}
     ]
   ];
  If[Head[points]===List,pointlist=points,(*else*)
     pointlist=Position[Map[points,P[poset]],True]//Flatten];
  card=Length[pointlist];
  links=ToRelation[mat[[pointlist,pointlist]]];
  {links,card}  
  ];
  
RestoreSubPosetLabels[poset_,points_,subposet_] := Block[{pointlist},
  If[Head[points]===List,pointlist=points,(*else*)
     pointlist=Position[Map[points,P[poset]],True]//Flatten];
  P[subposet]=Map[P[poset][[pointlist[[#]]]]&,P[subposet]];
  ];

IntervalP[name_,bottom_,top_] := Block[{bpos,tpos},
  {bpos,tpos} = 
    {Position[P[name],bottom],Position[P[name],top]}//Flatten;
  Intersection[Up[name][[bpos]],Down[name][[tpos]]]
  ]

JoinSubLattice[tag_,points_]:= JoinSubLattice[tag,points] =
  	FixedPoint[fjoinsublattice[tag,#]&,points];

fjoinsublattice[tag_,points_]:= 
	LJoin[tag][[points,points]]//Flatten//Union;

MeetSubLattice[tag_,points_]:= MeetSubLattice[tag,points] =
  	FixedPoint[fmeetsublattice[tag,#]&,points];

fmeetsublattice[tag_,points_]:= 
	LMeet[tag][[points,points]]//Flatten//Union;

SubLattice[tag_,points_]:= SubLattice[tag,points] =
  	FixedPoint[fsublattice[tag,#]&,points];

fsublattice[tag_,points_]:= 
	Union[
	   (LMeet[tag][[points,points]]//Flatten),
	   (LJoin[tag][[points,points]]//Flatten)]


CP[poset1_,poset2_] := CartesianProduct[poset1,poset2];

DS[poset1_,poset2_] := DisjointSum[poset1,poset2];

OS[poset1_,poset2_] := OrdinalSum[poset1,poset2];

SP[poset1_,points_] := SubPoset[poset1,points];

JSL[poset1_,points_] := JoinSubLattice[poset1,points];

MSL[poset1_,points_] := MeetSubLattice[poset1,points];

SL[poset1_,points_] := SubLattice[poset1,points];

(* Miscellaneous commands : Posets *)

Verbose = True;

ListCommands[string__] := 
  Block[{result},
    result = Names[StringJoin["Posets`*",string,"*"]];
    If[result == {}, Print["No Match."],
      If[OddQ[Length[result]],AppendTo[result,"---"]];
      result = Partition[result,Quotient[Length[result],2]];
      result = result//Transpose//MatrixForm]
  ];
  
  ListCommands[] := 
  Block[{result,search},
    result = Names["Posets`*"];
    If[OddQ[Length[result]],AppendTo[result,"---"]];
    result = Partition[result,Quotient[Length[result],2]];
    result = result//Transpose//MatrixForm

  ];

Zap[name_] := Block[{},
    Off[Unset::norep];
    P[name] =. ;
    PGraded[name] =. ;
    CoverRelations[name] =. ;
    Rank[name] =. ;
    RankedQ[name] =. ;
    NK[name] =. ;
    RGF[name] =. ;
    H[name] =. ;
    Card[name] =. ;
    ZetaP[name] =. ;
    Mu[name] =. ;
    Up[name] =. ;
    Down[name] =. ;
    Covers[name] =. ;
    CoCovers[name] =. ;
    MaxAntichain[name]=. ;
    DilworthCover[name]=. ;
    JI[name]=.;
    MI[name]=. ;
    JICoverRelations[name]=.;
    MICoverRelations[name]=. ;
    G[name]=. ;
    GBar[name]=.;
    Omega[name]=.;
    OmegaBar[name]=.;
    OmegaGF[name]=.;
    OmegaBarGF[name]=.;
    LinearExtensions[name]=. ;
    MaximalChainsUp[name]=. ;
    MaximalChainsDown[name]=. ;
    ChainsBetweenGF[name]=. ;
    LatticeQ[name]=. ;
    On[Unset::norep];
    ]
    

Locate[name_,element_] := 
    Position[P[name],element][[1,1]] 

LocateSet[name_,set_] := 
     Map[Position[P[name],#]&,set] //Flatten//Sort

IsomorphicQ[poset1_,poset2_]:=Block[{},
     
     If[Length[CoverRelations[poset1]]!=Length[CoverRelations[poset2]],
          Print["Different number of CoverRelations."];
          Return[False]
     ];
     If[Card[poset1]!=Card[poset2],
          Print["Cardinalities not equal."];
          Return[False]
     ];
     If[Sort[UpDegree[poset1]]!=Sort[UpDegree[poset2]],
          Print["UpDegree distribution is different."];
          Return[False]
     ];
     If[Sort[DownDegree[poset1]]!=Sort[DownDegree[poset2]],
          Print["DownDegree distribution is different."];
          Return[False]
     ];
     If[Count[Flatten[ZetaP[poset1]],1]!=
          Count[Flatten[ZetaP[poset2]],1],
          Print["Number of ones in Zeta matrix is different."];
          Return[False]
     ];
     If[Sort[Flatten[Mu[poset1]]]!=
          Sort[Flatten[Mu[poset2]]],
          Print["Some Mobius function values are different."];
          Return[False]
     ];
     If[RankedQ[poset1]!=RankedQ[poset2],
          Print["One is ranked; the other is not."];
          Return[False]
     ];     
     If[Rank[poset1]!=Rank[poset2],
          Print["Rank distribution is different."];
           Return[False]
     ];     
Return[Global`Probably]
]

NumPosets[n_] := Length[Global`AllPosets[n]];
NumLattices[n_] := Length[Global`AllLattices[n]];

(* Miscellaneous Commands: Combinatorics *)

InfiniteProduct[F_,n_] := Block[{p=F-1,log,loglist,k,a,m,d},
    log = Expand[Sum[(-1)^(k-1) p^k /k,{k,1,n}]];
    loglist =CoefficientList[log,q];
    a[m_]:=Sum[d loglist[[d+1]] MoebiusMu[m/d]If[IntegerQ[m/d],1,0],
                         {d,1,n}] /m ;
    Array[a,n]
    ]

InversionCode[s_] := Table[Count[
    Drop[s,Position[s,i][[1,1]]],
    x_Integer?(i>#&)],{i,Length[s]}];

ToPermutation[code_] := Block[{out={},k},
  For[k=1,k<=Length[code],k++,
    out=Insert[out,k,k-code[[k]]]
    ];out]  

InversePermutation[p_]:=
     ({p,Range[Length[p]]}//Transpose//Sort//Transpose)[[2]]

ZeroDiffs[seq_]:=Block[{current=seq,zerodiffs={seq[[1]]}},
  While[Length[current]!=1 && current*0 !=current,
  current=Drop[current,1]- Drop[current,-1];
  AppendTo[zerodiffs,current[[1]] ]];
  zerodiffs
]

PolyToGF[f_,x_,q_] := Block[{values,d,zerodiffs,i},
  d = Exponent[f,x];
  zerodiffs = ZeroDiffs[Table[f,{x,0,d}]];
  Sum[zerodiffs[[i+1]]*q^i * (1-q)^(d-i),
            {i,0,d}]/(1-q)^(d+1)      
  ] 

GFToPoly[w_,q_,d_,n_] := Block[
   {r,a=CoefficientList[w /. q->1-r,r]},
   Sum[a[[i]]Binomial[n+d-(i-1),d-(i-1)],{i,1,Length[a]}] ];

GenFun[L_,q_] := Block[{},
    L.Table[q^(i-1),{i,Length[L]}]];

Version = "2.2d4b";

                       End[];

                   EndPackage[];