// Permutation Groups, Mathematics Computing Project, March 2005
// Martin Jaggi m.jaggi@gmx.net
//
// Implementation of the Incremental Version of the Schreier Sims Algorithm
//
#include <iostream>
using namespace std;

const int maxn = 100;		// Maximum size of omega = {1, ,n}
const int maxr = 500;		// Maximum number of generators for the group

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; };
		int fbp;	// to store the first basepoint moved by this permutation
	private:
		int n;								// size of omega = {1, ,n}
};

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

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

int cosreps;								// number of    (= size of orbit of alpha)
Permutation cosrep[maxn];					// coset representatives (to store output of SchreierTree)
Permutation csrep[maxn];					// a copy of these
int firstbasepoint = 0;
int basepoints;
int* basepoint = new int[maxn];				// base
int sgss;
Permutation* sgs = new Permutation[maxr];	// strong generating set


void output();

		/************/
		int SmallestPointMovedBy(Permutation h) {		// self explaining
			for (int i=0;i<n;i++) if (h.p[i] != i) return i; return -1;
		}
		int FirstBasePointMovedBy(Permutation h) {		// self explaining
			for (int i=0;i<basepoints;i++) if (h.p[ basepoint[i] ] != basepoint[i]) return i; return -1;
		}


Permutation gen (n,'i');	// group element moving original alpha to actual alpha

/****** SchreierTree ******/
void SchreierTree(int alpha, int i) {		// depth first search to determine the orbit of alpha under G(i)
	static int ag;							// image of actual alpha under generator sgs[j]
	cosrep[alpha] = gen; cosreps++;
	for (int j=0;j<sgss;j++) { if (sgs[j].fbp >= i) {
		ag = sgs[j].p[alpha];
		if (!cosrep[ag].isdefined()) {
			gen *= sgs[j]; SchreierTree(ag,i); gen *= sgs[j].inverse();
		}
	} }
}

/****** Sifting ******/
bool Sifting(Permutation h) {	// returns if g is in G, given by SGS sgs
	int i;
	if (h.isidentity()) return true;
	else {
		i = FirstBasePointMovedBy(h);
		if (i<firstbasepoint) return false; else {
			cosreps = 0; for (int j=0;j<n;j++) cosrep[j] = undefined; gen = identity;
												// construct the schreier tree with root basepoint[i] under G(i)
			SchreierTree(basepoint[i],i);
			
			if (!cosrep[ h.p[ basepoint[i] ] ].isdefined()) return false;	// basepoint[i]^h is not in orbit
			else {
				return Sifting( h * cosrep[ h.p[basepoint[i]] ].inverse() );
			}
		}
	}
}

/****** SchreierSims ******/
void SchreierSims() {		// schreier sims algorithm
	Permutation s, h;
	
	do {		
		bool stabilizer_is_nontrivial = false;
		for (int j=0;j<sgss;j++)
			if (sgs[j].fbp > firstbasepoint) stabilizer_is_nontrivial = true;
		
		/*** recursive call ***/
		firstbasepoint++;
		if (stabilizer_is_nontrivial) SchreierSims();
		firstbasepoint--;
		
		// now test each schreier generator s for G(b1) to see if it is in H = <new sgs>
					
									// get generators for H an coset representatives
		cosreps = 0; for (int j=0;j<n;j++) cosrep[j] = undefined; gen = identity;
		SchreierTree(basepoint[firstbasepoint],firstbasepoint);	// determine orbit of b1 under G(firstbasepoint)
		for (int j=0;j<n;j++) csrep[j] = cosrep[j];
		
		s = undefined;
		for (int i=0;i<n;i++) { if (csrep[i].isdefined())
			for (int j=0;j<sgss;j++) { if (sgs[j].fbp >= firstbasepoint) {	// generators for G(firstbasepoint)
				h = csrep[i] * sgs[j] * csrep[ sgs[j].p[i] ].inverse();		// schreier generator for G(b1)
				if (!Sifting( h )) { s = h; /*cout<<"outside H: ";s.output();*/ }
			} }
		}
			
		if (s.isdefined()) {
			sgs[sgss] = s;	// add s to sgs
										//cout<<"s outside H found: "; cout<<endl;//s.output();			
				// add point moved by s to base if s fixes all base elements
			if (FirstBasePointMovedBy(s) == -1) {
										//cout<<"s fixes all basepoints, added "<<SmallestPointMovedBy(s)+1<<"\n";
				basepoint[basepoints++] = SmallestPointMovedBy(s);
			}
			sgs[sgss].fbp = FirstBasePointMovedBy(s);
			sgss++;
		}
		else { //cout << "finished since H = Gb1 \n";
			return;
		}
	} while (true); /*** goto beginning ***/

}

/****** main ******/
int main () {
	cout << "n ( Size of Omega = {1..n} ) ? "; cin >> n;
			
			// set the size of all permutations to n
			for (int j=0;j<maxr;j++) { sgs[j].setn(n); }
			identity.setn(n); undefined.setn(n);
	
	cout << "How many group generators ? "; cin >> sgss;
	for (int j=0;j<sgss;j++) sgs[j].input();
	
	// compute initial partial base and partial strong generating set
	basepoints = 0; int i; bool is_in_base_yet;
	for (int j=0;j<sgss;j++) {
		if (sgs[j].isidentity()) {		// delete all identities
			sgss--;
			for (int k=0;k<sgss;k++) sgs[k] = sgs[k+1];
		}
		i = SmallestPointMovedBy(sgs[j]); is_in_base_yet = false;
		for (int k=0;k<basepoints;k++) if (basepoint[k] == i) is_in_base_yet = true;
		if (!is_in_base_yet) basepoint[basepoints++] = i;
	}
	for (int j=0;j<sgss;j++) sgs[j].fbp = FirstBasePointMovedBy(sgs[j]);
	
	// start the schreier sims algorithm
	firstbasepoint = 0;
	SchreierSims();
	
	//cout<<"#### output of the incremental schreier sims algorithm: ####\n";
	//output();
	
	// finally calculate the order of the group using the sgs
	
	cout<<"#### the order of the group: ####\n";
	i = basepoint[0];
	for (int i=0;i<basepoints;i++) {
			cosreps = 0; for (int j=0;j<n;j++) cosrep[j] = undefined; gen = identity;
													// construct the schreier tree with root basepoint[i]
			SchreierTree(basepoint[i],i);
			cout << cosreps << "*";
	}
	cout<<endl;
				
	delete[] sgs;
	return 0;
}


/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */


/****** Output base and sgs ******/
void output() {
	for (int j=0;j<firstbasepoint;j++) cout<<"  ";
	cout<<"base ("<<basepoints - firstbasepoint<<") : ";
	for (int j=firstbasepoint;j<basepoints;j++) cout<<basepoint[j]+1<<" "; cout<<endl;
	for (int j=0;j<firstbasepoint;j++) cout<<"  ";
	cout<<"sgs ("<<sgss<<") :\n";
	for (int j=0;j<sgss;j++) {		
		for (int i=0;i<firstbasepoint;i++) cout<<"  ";
		sgs[j].output();
	}
}
/*
int schreiergs;									// number of
Permutation* schreierg = new Permutation[maxr];	// schreier generators (to store output of SchreierLemma)
****** SchreierLemma ******
void SchreierLemma() {					// generators for the stabilizer of alpha
	for (int i=0;i<n;i++) { if (cosrep[i].isdefined())
		for (int j=0;j<r;j++) {
			schreierg[schreiergs++] = cosrep[i] * g[j] * cosrep[ g[j].p[i] ].inverse();
		}
	}
}
*/