// Permutation Groups, Mathematics Computing Project, March 2005
// Martin Jaggi m.jaggi@gmx.net
//
// Implementation of the Schreier Sims Algorithm using Jerrums Filter
// using only O( n^2 ) in space
//
#include <iostream>
using namespace std;

const int maxn = 2100;		// Maximum size of omega = {1, ,n}

class Permutation {			// interface for permutations
	public:
		int p[maxn];														// the images of the points 0..maxn-1
		Permutation () { n=maxn; };											// constructors
		Permutation (int m) { n=m; };
		Permutation (int m, char c) { n=m; switch(c) {
			case 'i': for (int i=0;i<n;i++) p[i] = i; break; // identity
			default:  for (int i=0;i<n;i++) p[i] =-1; break; // undefined
		} }
		Permutation operator * (Permutation param) const {					// multiplication
			Permutation result (n);
			for (int i=0;i<n;i++) result.p[i] = param.p[ p[i] ];
			return (result);
		}
		void operator *= (Permutation param) {								// direct multiplication
			for (int i=0;i<n;i++) p[i] = param.p[ p[i] ];
		}						
		Permutation inverse() const {										// inverse
			Permutation result (n);
			for (int i=0;i<n;i++) result.p[ p[i] ] = i;
			return (result);
		}
		bool isdefined() const { return (p[0] > -1); }						// if it is defined
		bool isidentity() const {											// if it is the identity
			for (int i=0;i<n;i++) if (p[i] != i) return false; return true;
		}				
		bool operator == (Permutation param) const {						// comparison
			for (int i=0;i<n;i++) if (param.p[i] != p[i]) return false; return true;
		}				
		void input()        { for (int i=0;i<n;i++) { cin >> p[i]; p[i]--; } }			   // input
		void output() const { for (int i=0;i<n;i++) cout << p[i]+1 << " "; cout << endl; } // output
		void setn(int m) { n=m; }
	private:
		int n;								// size of omega = {1, ,n}
};

/****** globals ******/
int n;										// size of omega = {1, ,n}

int r;										// number of generators
Permutation* g;								// the generators

int cosreps;								// number of    (= size of orbit of alpha)
Permutation* cosrep;						// coset representatives (to store output of SchreierTree)

Permutation undefined (maxn,'u');
Permutation identity (maxn,'i');


/****** ScheierTree ******/
void ScheierTree(int alpha) {				// depth first search to determine the orbit of alpha
	static Permutation gen (n,'i');	// group element moving original alpha to actual alpha
	static int ag;					// image of actual alpha under generator g[j]
	cosrep[alpha] = gen; cosreps++;
	for (int j=0;j<r;j++) { ag = g[j].p[alpha];
		if (!cosrep[ag].isdefined()) {
			gen *= g[j]; ScheierTree(ag); gen *= g[j].inverse();
		}
	}
}
	
		/****** JerrumsFilter ******/	
		// globals
		struct JVertexT {			// a vertex of jerrum's graph
			int neighbors;					// adjacent edges
			int neighbor[maxn-1];			// list of neighbors. neighbors[i] points to the i-th neighbor
			int perm[maxn-1];				// perm[i] points to the permutation jg describing the edge to to the i-th neighbor (which is jg[ perm[i] ])
			int isneighbor[maxn];			// if i is a neighbor, where it is in the list of neighbors
			bool visited;
		};
		JVertexT* JVertex;		// jerrums graph
		Permutation* jg;		// jerrum's generators (edges of the graph)
		int jgs;
		int freejg[maxn];		// empty slots of the above array. freejg[i] points to the i-th free position in our list of generating permutations
		int freejgs;
		int cycle[maxn]; int length=0;				// here we will store our cycle of length length if we find one
		Permutation h (n);
		
		int SmallestPointMovedBy(Permutation h) {		// self explaining
			for (int i=0;i<n;i++) if (h.p[i] != i) return i; return -1;
		}
		bool FindCycle(int v, int wherewecomefrom) {	// finds a cycle in the graph starting from vertex v
			static int u;
			for (int i=0;i<JVertex[v].neighbors;i++) {
				u = JVertex[v].neighbor[i];
				if (u != wherewecomefrom) {
					if (u != cycle[0]) {
						cycle[ length++ ] = u;
						JVertex[u].visited = true;
						if (FindCycle(u,v)) return true;
						length--;
					} else {
						return true;
					}
				}
			} return false;
		}
		void InsertEdge(Permutation g, int i, int ig) { // insert edge (permutation) moving i to ig (i<ig)
			int nn,jgslot;
			jgs++;
			jgslot = freejg[--freejgs];	// the last free generator place
			jg[jgslot] = g;				
							// now update the two vertices i
				nn = JVertex[i ].neighbors++;
			JVertex[i ].neighbor[ nn ] = ig;
			JVertex[i ].perm[ nn ] = jgslot;
			JVertex[i ].isneighbor[ ig ] = nn;		// means "ig is i's neighbor number nn"
							// and ig
				nn = JVertex[ig].neighbors++;
			JVertex[ig].neighbor[ nn ] = i;
			JVertex[ig].perm[ nn ] = jgslot;
			JVertex[ig].isneighbor[ i ] = nn;
		}
		void RemoveEdge(int i, int ig) {			// remove edge from i to ig
							// remove the permutation of this edge
			jgs--; int jgslot = JVertex[i].perm[ JVertex[i].isneighbor[ig] ];			
			jg[jgslot].p[0] = -1; // just a faster writing for jg[jgslot] = undefined; 
			freejg[freejgs++] = jgslot;	// now we gained a free place
			
							// now update the two vertices i
			JVertex[i].neighbors--;
			for (int j=JVertex[i].isneighbor[ig];j<JVertex[i].neighbors;j++) {
				JVertex[i].neighbor[j] = JVertex[i].neighbor[j+1];
				JVertex[i].perm[j]     = JVertex[i].perm[j+1];
				JVertex[i].isneighbor[ JVertex[i].neighbor[j] ] = j;
			} JVertex[i].isneighbor[ig] = -1;
							// and ig
			JVertex[ig].neighbors--;
			for (int j=JVertex[ig].isneighbor[i];j<JVertex[ig].neighbors;j++) {
				JVertex[ig].neighbor[j] = JVertex[ig].neighbor[j+1];
				JVertex[ig].perm[j]     = JVertex[ig].perm[j+1];
				JVertex[ig].isneighbor[ JVertex[ig].neighbor[j] ] = j;
			} JVertex[ig].isneighbor[i] = -1;
		}
		
		void ThrowIn(Permutation g) {		// throws permutation into the graph and does the required householding steps
			static int i,ig,nn,c,cmin,min,j;
			
			if (!g.isidentity()) {
				i = SmallestPointMovedBy(g); ig = g.p[i];
				nn = JVertex[i].isneighbor[ ig ];
			
				if (nn>=0) {	// if ig is already neighbor of i
					if (!(g == jg[ JVertex[i].perm[nn] ])) {	// if we not already have exactly this generator
														// we have a 'trivial' cycle consisting of two edges.
														// replace gi1 by g and insert h = gi1 * g.inverse()
						h = jg[ JVertex[i].perm[nn] ] * g.inverse();
						
						jg[ JVertex[i ].perm[nn] ] = g;	// replace this edge by our new g;
										
						ThrowIn(h);						// and throw in h
					}
				} else {	// we don't have this edge yet
					InsertEdge(g,i,ig);
																							
					for (c=0;c<n;c++) JVertex[c].visited = false;
					JVertex[i].visited=true; cycle[0] = i; length = 1; 
						
					if (FindCycle(i,-1)) { 					
						cycle[length] = cycle[0];
														// find the smallest point (vertex) of the cycle
						min=n;
						for (c=0;c<length;c++) if (cycle[c]<min) { min=cycle[c]; cmin=c; }
						
						h = identity;					// now walk around the cycle starting from the smallest point, and multiply around the cycle to obtain h
						for (c=cmin;c<length;c++) { i=cycle[c]; j=cycle[c+1];
							if (i < j) { nn = JVertex[i].isneighbor[j];  h *= jg[ JVertex[i].perm[nn] ];           }
							else {       nn = JVertex[j].isneighbor[i];  h *= jg[ JVertex[j].perm[nn] ].inverse(); }
						} for (c=0;c<cmin;c++) {    i=cycle[c]; j=cycle[c+1];
							if (i < j) { nn = JVertex[i].isneighbor[j];  h *= jg[ JVertex[i].perm[nn] ];           }
							else {       nn = JVertex[j].isneighbor[i];  h *= jg[ JVertex[j].perm[nn] ].inverse(); }
						}
														// remove the first edge of the cicle
						RemoveEdge(cycle[cmin],cycle[cmin+1]);
														// and throw in h instead
						ThrowIn(h);				
					}
				} // if ig is already neighbor of i
			} // if g is not the identity
		}
			
/****** JerrumsFilter ******/
void JerrumsFilter() {						// filters the generators

	// initialisation of jerrum's graph and jerrum's generators (jg)
	jgs = 0; freejgs = n; for (int j=0;j<n;j++) { jg[j] = undefined; freejg[j] = j; }
	for (int i=0;i<n;i++) JVertex[i].neighbors=0;
	for (int i=0;i<n;i++) for (int j=0;j<n;j++) JVertex[i].isneighbor[j]=-1;
	
	// schreier lemma loop to get the schreier generators			
	for (int i=0;i<n;i++) { if (cosrep[i].isdefined())
		for (int j=0;j<r;j++) {
		
			// throw the (schreier) generators into the graph, one after each other
			ThrowIn(  cosrep[i] * g[j] * cosrep[ g[j].p[i] ].inverse()  );
		}
	}
	
	r = 0;
	for (int j=0;j<n;j++) if (jg[j].isdefined()) g[r++] = jg[j];
}

/****** SchreierSims ******/
void SchreierSims() {					// schreier sims algorithm
	int alpha = 0;
						cout << "THE ORDER OF THE GROUP:\n";
	do {
		cosreps = 0; for (int i=0;i<n;i++) cosrep[i] = undefined;
										// get the coset representatives for G(alpha)
		ScheierTree(alpha);
										// now perform jerrum's filter for the schreier generators
		JerrumsFilter();
			//cout << "Generators found by jerrum's filter:\n"; for (int j=0;j<r;j++) g[j].output();
		
		alpha++;
		
		cout << cosreps << flush; if (r>0) cout<<"*";
	} while (r>0);
	cout << endl;
}

/****** main ******/
int main () {
	
	cout << "n ( Size of Omega = {1..n} ) ? "; cin >> n;
	
			g = new Permutation[n];
			jg = new Permutation[n];
			cosrep = new Permutation[n];
			JVertex = new JVertexT[n];
			
			// set the size of all permutations to n
			for (int j=0;j<n;j++) { g[j].setn(n); jg[j].setn(n); }
			identity.setn(n); undefined.setn(n);
				
	cout << "How many group generators ? "; cin >> r;
	for (int j=0;j<r;j++) g[j].input();
	
	SchreierSims();
	
			delete[] g; delete[] jg; delete[] cosrep; delete[] JVertex;
	return 0;
}