#include "FiniteLinearGroup.h"
#include "FiniteGroup.h"
#include <cmath>
#include <iostream>
#include <fstream>
#include <string>

using namespace std;
using namespace giac;

FiniteLinearGroup::FiniteLinearGroup(CMfield k) {
  K = k;
  elements.push_back(K.Id2);
  oldnumber=0;
};

bool FiniteLinearGroup::isInList(gen M, vector<gen> list) {
  bool res = false;
  int k=0;
  while (!res && k<list.size()) {
    res = operator_equal(M,list[k],&ct);
    k++;
  }
  return res;
}

void FiniteLinearGroup::addElement(gen M) {
  if (!isInList(M,elements)) {
    elements.push_back(M);
  }
}

void FiniteLinearGroup::addElements(vector<gen> list) {
  for (unsigned k=0; k<list.size(); k++) {
    addElement(list[k]);
  }  
}

void FiniteLinearGroup::expand() {
  //  cout << "expand..." << endl;
  vector<gen> newlist;
  for (unsigned k=oldnumber; k<elements.size(); k++) {
    for (unsigned l=1; l<elements.size(); l++) {
      gen M = elements[k]*elements[l];
      K.matred(M);
      if (!isInList(M,elements) && !isInList(M,newlist)) {
	//	cout << "New matrix: " << M << endl;
	newlist.push_back(M);
      }
    }
  }
  for (unsigned l=oldnumber; l<elements.size(); l++) {
    for (unsigned k=1; k<elements.size(); k++) {
      gen M = elements[k]*elements[l];
      K.matred(M);
      if (!isInList(M,elements) && !isInList(M,newlist)) {
	//	cout << "New matrix: " << M << endl;
	newlist.push_back(M);
      }
    }
  }
  oldnumber = elements.size();
  //  cout << "old nb: " << oldnumber << endl;
  addElements(newlist);
  cout << "Current nb of elements: " << elements.size() << endl;
}

void FiniteLinearGroup::expandToGroup() {
  bool done = false;
  int count = 0;
  while (!done && count<20) {
    int ne = elements.size();
    //    cout << "ne=" << ne << endl;
    expand();
    done = (ne==elements.size());
    count++;
  }
  if (done) {
    cout << "Found a group with " << elements.size() << " elements" << endl;
  }
  else {
    cout << "Did not find a group, even after expanding 20 times..." << endl;    
  }
}

bool FiniteLinearGroup::identify_root_of_unity(gen mu, int o, gen &res) {

  if (o==1) {
    res = K.onenumber;
    return true;
  }
  if (o==2) {
    res = K.onenumber/2;
  }
  gen anro = eval(_subst(makesequence(mu,x__IDNT_e,K.genasrootof),&ct),&ct);
  gen te0 = _simplify(_eval(anro,&ct),&ct);
  if (o==4) {
    if (is_greater(_evalf(im(anro,&ct),&ct),K.zeronumber,&ct)) {
      res = K.onenumber/4;
      return true;
    }
    else {
      if (is_greater(K.zeronumber,_evalf(im(anro,&ct),&ct),&ct)) {
	res = -K.onenumber/4;
	return true;
      }
    }
  }
  
  /*  cout << "te0=" << te0 << endl;
  cout << "convert te0 interval: " << convert_interval(te0,K.nd,&ct) << endl;
  cout << "arg/2pi: " << arg(convert_interval(te0,K.nd,&ct),&ct)/(2*gen("pi",&ct)) << endl;*/
	
  gen te = re(arg(convert_interval(eval(te0,1,&ct),K.nd,&ct),&ct)/(2*gen("pi",&ct)),&ct);
  //  cout << "te=" << te << endl;
  //  cout << "o=" << o << endl;
  // need to identify this as a rational of the form k/o, k=0,...,o-1
  int count = 0;
  for (int k=-o+1; k<o ; k++) {
    gen gk = gen(k);
    gen kogi = gk*(K.onenumber/gen(o));
    gen zi = convert_interval(kogi,K.nd,&ct);
    //    cout << "Comparing with " << zi << endl;
    gen test = zi-te;
    if (is_greater(K.zeronumber,_left(test,&ct),&ct) && is_greater(_right(test,&ct),K.zeronumber,&ct)) {
      count++;
      res = kogi;
    }
  }
  if (count!=1) {
    cout << "Trouble identifying argument of a root of unity, count=" << count << endl;
    cout << "  (this will have no consequences, except possibly for longer runtime)" << endl;
    return false;
  }
  else {
    //    cout << "Success, " << res << endl;
    return true;
  }
}

bool FiniteLinearGroup::studyReflectionSubgroup(int &ord_gp){
  int ne = elements.size();
  vector<int> refl_inds;
  vector<int> refl_orders;
  vecteur mults;
  //  cout << "Reflections in stab:" << endl;
  for (int j=0; j<ne; j++) {
    gen A = elements[j]-K.Id2;
    gen d = A[0][0]*A[1][1]-A[0][1]*A[1][0];
    K.red(d);
    if (!K.isSameMatrix(elements[j],K.Id2) && operator_equal(d,K.zeronumber,&ct)) {
      refl_inds.push_back(j);
      int o = K.linearorder2d(elements[j]);
      refl_orders.push_back(o);
      // 1 is an eigenvalue, so the other eigenvalue is the determinant
      gen mu = elements[j][0][0]*elements[j][1][1] - elements[j][1][0]*elements[j][0][1];
      K.red(mu);
      gen rat;
      if (!identify_root_of_unity(mu,o,rat)) {
	mults.push_back(-K.onenumber);
      }
      else {
	mults.push_back(rat);
      }
    }
  }
  int nr = refl_inds.size();
  cout << "Found " << nr << " reflections" << endl;
  cout << "Orders: " << refl_orders << endl;
  cout << "Multipliers, exp(2*pi*i*...): " << mults << endl;
  cout << endl;
  // NOT ENOUGH TO USE 2 GENERATORS IN GENERAL!

  if (nr==elements.size()-1){
    cout << "Generators are all reflections!" << endl;
    cout << endl;
  }
  else {
    cout << "There are some non-reflection generators..." << endl;
    cout << endl;
  }
  
  int gmax = 1, cmax = 1;
  int jmax = 0;
  int kmax = 0;
  int bomax;
  bool found_one_pair = false;
  bool invjmax, invkmax;
  gen Amax, Bmax;
  int ojmax, okmax;
  
  for (int j=0; j<nr; j++) {
    int oj = K.linearorder2d(elements[refl_inds[j]]);
    //    cout << "numer j: " << _numer(mults[j],&ct) << endl;
    bool tjp = operator_equal(_numer(mults[j],&ct),K.onenumber,&ct); 
    bool tjm = operator_equal(_numer(mults[j],&ct),-K.onenumber,&ct); 
    for (int k=j+1; k<nr; k++) {
      int ok = K.linearorder2d(elements[refl_inds[k]]);
      //     cout << "numer k: " << _numer(mults[k],&ct) << endl;
      bool tkp = operator_equal(_numer(mults[k],&ct),K.onenumber,&ct); 
      bool tkm = operator_equal(_numer(mults[k],&ct),-K.onenumber,&ct); 
      if ((tjp || tjm) && (tkp || tkm)) {

	gen A,B;
	if (tjp) {
	  A = elements[refl_inds[j]];
	}
	else {
	  A = K.matinverse2d(elements[refl_inds[j]]);
	}
	if (tkp) {
	  B = elements[refl_inds[k]];
	}
	else {
	  B = K.matinverse2d(elements[refl_inds[k]]);
	}
	
	int bo = K.braidOrder2d(A,B);
	//	cout << "bo=" << bo << endl;
	//	if (bo>1) {
	if (bo>2) {
	  // I think I have bo>2 only to avoid taking reflections with the same mirror...

	  //	  cout << j << "," << k << endl;
	  //	  cout << "reflections of order " << oj << ", " << ok << ", that braid with length " << bo << endl;
	  gen poj = oj;
	  gen pok = ok;
	  gen pbo = bo;
	  gen ms = simplify((K.onenumber/poj + K.onenumber/pok + K.twonumber/pbo - K.onenumber),&ct);
	  if ( operator_equal(ms,K.zeronumber,&ct) || !operator_equal(_numer(ms,&ct),K.onenumber,&ct) ) {
	    cout << "Trouble computing order of reflection group, " << ms << " should be the inverse of an integer" << endl;
	  }
	  else {
	    ms = K.onenumber/ms;
	    gen oc;
	    if (bo%2==0) {
	      oc = (K.twonumber*K.twonumber*ms)/pbo;
	    }
	    else {
	      oc = (K.twonumber*ms)/pbo;
	    }
	    gen og = 8*K.onenumber*ms*ms/pbo;
	    if (!operator_equal(_denom(og,&ct),K.onenumber,&ct) || !operator_equal(_denom(oc,&ct),K.onenumber,&ct)) {
	      cout << "Trouble computing order of reflections group, " << og << ", " << oc << " should be an integer..." << endl;
	    }
	    else {
	      found_one_pair = true;
	      //	      cout << "og=" << og << endl;
	      //	      cout << "oc=" << oc << endl;
	      int g = _numer(og,&ct).val;
	      int c = _numer(oc,&ct).val;
	      //	      cout << "g=" << g << endl;
	      //	      cout << "c=" << c << endl;
	      if (g>gmax) {
		gmax = g;
		cmax = c;
		ojmax = oj;
		okmax = ok;
		jmax = j;
		kmax = k;
		invjmax = tjm;
		invkmax = tkm;
		Amax = A;
		Bmax = B;
		cout << "Found two reflections of order " << oj << ", " << ok << ", that braid with length " << bo << endl;
		bomax = bo;
	      }	    
	    }
	  }
	}
	if (bo==2) {
	  // Need to check if the mirrors are the same
	  //	  cout << "Found commuting reflections, of orders " << oj << "," << ok << endl;

	  gen MA = A-K.Id2;
	  bool fda = false;
	  int k=0;
	  gen mia;
	  while (!fda && k<2) {
	    gen v0 = _row(makesequence(MA,k),&ct);
	    fda = ( !operator_equal(v0[0],K.zeronumber,&ct) || !operator_equal(v0[1],K.zeronumber,&ct) );
	    if (fda) {
	      vecteur te;
	      te.push_back(-v0[1]);
	      te.push_back(v0[0]);
	      mia = gen(te);
	    }
	    k++;
	  }
	  
	  gen MB = B-K.Id2;	  
	  bool fdb = false;
	  k=0;
	  gen mib;
	  while (!fdb && k<2) {
	    gen v0 = _row(makesequence(MB,k),&ct);
	    fdb = ( !operator_equal(v0[0],K.zeronumber,&ct) || !operator_equal(v0[1],K.zeronumber,&ct) );
	    if (fdb) {
	      vecteur te;
	      te.push_back(-v0[1]);
	      te.push_back(v0[0]);
	      mib = gen(te);
	    }
	    k++;
	  }

	  if (fda && fdb) {
	    gen de = mia[0]*mib[1] - mia[1]*mib[0];
	    K.red(de);
	    if (!operator_equal(de,K.zeronumber,&ct)) {
	      //	      cout << " ... mirrors are distinct" << endl;
	      found_one_pair = true;
	      int g = oj*ok;
	      int c = g;
	      if (g>gmax) {
		gmax = g;
		cmax = c;
		jmax = j;
		kmax = k;
		ojmax = oj;
		okmax = ok;
		invjmax = tjm;
		invkmax = tkm;
		Amax = A;
		Bmax = B;
		cout << "Found two reflections of order " << oj << ", " << ok << ", that commute " << endl;
		bomax = bo;
	      }
	    }
	    /*	    else {
	      cout << " ... they have the same mirror" << endl;
	      }*/
	  }
	  /*	  else {
	    cout << "Trouble computing mirrors" << endl;
	    cout << "MA=" << MA << endl;
	    cout << "MB=" << MB << endl;
	    }*/
	}
      }
    }
  }

  if (!found_one_pair) {
    cout << "Did not find any coherent pair of reflections" << endl;
    return false;
  }

  cout << "Found a subgroup generated by two reflections, of order " << gmax << endl;
  cout << "                                  center of order " << cmax << endl;
  //  cout << "  (in principle, could get a larger one with three generators)" << endl;
  cout << endl;
  //  cout << "Indices of gens: " << refl_inds[jmax] << "," << refl_inds[kmax] << endl;

  if (bomax>2) {
    // what if we don't find any reflections satisfying the coherence test? (will crash in that case)
    //  gen A = elements[refl_inds[jmax]];
    //  gen B = elements[refl_inds[kmax]];
    /*  if (take_inverses) {
	A = K.matinverse(A);
	B = K.matinverse(B);
	}*/
    gen C = Amax*Bmax;
    K.matred(C);
    gen Z; // generator of the center of <A,B>
    if (bomax%2==0) {
      Z = _pow(makesequence(C,bomax/2),&ct);
    }
    else {
      Z = _pow(makesequence(C,bomax),&ct);
    }

    /*  cout << "Order A: " << K.linearorder2d(A) << endl;
	cout << "Order B: " << K.linearorder2d(B) << endl;
	cout << "Braid length: " << K.braidOrder2d(A,B) << endl;
	cout << "Order Z: " << K.linearorder2d(Z) << endl;
	cout << "cmax=" << cmax << endl;*/
  
    vecteur z;
    z.push_back(elements[0]);
    gen M = Z;
    for (int j=1; j<cmax; j++) {
      z.push_back(M);
      M = M*Z;
      K.matred(M);
    }
  //  cout << "size z: " << z.size() << endl;
  
    cout << "Now studying the projective action of the stabilizer" << endl;
    cout << endl;
    
    // now want to check if every element in the stab generator is accounted for
    FiniteGroup S = FiniteGroup(K,2);
    //  cout << "Creating 2D Projective Group" << endl;
    S.addElements(elements);
    vecteur h;
    vecteur tsv;
    if (S.exploreCosets(h,tsv)) {
      
      /*    cout << "Gen for H:" << endl;
	    cout << h[1] << endl;
	    cout << "tsv:" << endl;
	    cout << tsv;
	    cout << "Z:" << endl;
	    cout << Z << endl;*/
      
      //    cout << "Will try to construct:" << endl;
      vector<int> inds_left;
      for (int j=1; j<elements.size(); j++) {
	if (j!=refl_inds[jmax] && j!=refl_inds[kmax]) {
	  inds_left.push_back(j);
	  //	cout << j << ": " << elements[j] << endl;
	}
      }
      vecteur li;
      bool is_refl = false;
      int j=0;
      while (!is_refl && j<tsv.size()) {
	int k=0;
	while (!is_refl && k<h.size()) {
	  int l=0;
	  gen M = tsv[j]*h[k];
	  K.matred(M);
	  // need to check if M has the same projective action as some element in inds_left!
	  bool fd0 = false;
	  vector<int> ml;
	  for (int m = 0; m<inds_left.size(); m++) {
	    if (S.areMultiple(M,elements[inds_left[m]])) {
	      //	    cout << M << " is multiple of " << elements[inds_left[m]] << endl;
	      ml.push_back(m);
	    }
	  }
	  if (ml.size()>0) {
	    //	  cout << j << "," << k << " gives the correct projective action for " << ml.size() << " generators" <<  endl;
	    //	  cout << M << endl;
	    for (int u=0; u<ml.size(); u++) {
	      bool fd = false;
	      int l = 0;
	      while (!fd && l<z.size()) {
		gen N = M*z[l];
		K.matred(N);
		//	      cout << "After multiplying by Z^" << l << ":" << endl;
		//	      cout << N << endl;
		fd = K.isSameMatrix(N,elements[inds_left[ml[u]]]);
		l++;
	      }
	      if (!fd) {
		if (nr==elements.size()-1) {
		  cout << "Stabilizer is a reflection group (actually all generators are reflections)" << endl;
		  cout << "  but it can only be generated by more than two elements..." << endl;
		  cout << "I don't know how an efficient method to describe the center" << endl;
		  cout << "For now I need to enumerate all group elements, this can take a long time" << endl;
		  return false;
		}
		else {
		  cout << "The stabilizer is either not well-generated, or not a reflection group..." << endl;
		  cout << "For now I need to enumerate all group elements, this can take a long time" << endl;
		  return false;
		}
	      }
	    }
	    for (int u=ml.size()-1; u>-1; u--) {
	      inds_left.erase(inds_left.begin()+ml[u]);
	    }
	    //	  cout << "inds_left: " << inds_left << endl;
	  }
	  /*	else {
		cout << j << "," << k << " useless" << endl;
		cout << M << endl;
		}*/
	  is_refl = (inds_left.size()==0);
	  k++;
	}
	j++;
      }
      if (is_refl) {
	cout << "Found all generators using projective action and center computed from 2-generator reflection subgroup" << endl;
	cout << "  the reflection group we just found is the full stabilizer!" << endl;
	cout << endl;
	ord_gp = gmax;
	return true;
	//      cout << "Stabilizer is a reflection group generated by two reflections" << endl;
	//     cout << "  it has order " << gmax << ", center of order " << cmax << endl;
      }
      else {
      cout << "Trouble studying projectivized reflection group" << endl;
      cout << endl;
      cout << "Generators we did not find in 2-generator subgroup:" << endl;
      for (int j=0; j<inds_left.size(); j++) {
		cout << elements[inds_left[0]] << endl;
      }		
      if (nr==elements.size()-1) {
	cout << "Group is not well-generated..." << endl;
	cout << "I don't know how an efficient method to describe the center" << endl;
	cout << "For now I need to enumerate all group elements, this take a very long time" << endl;
      }
      else {
	cout << "The stabilizer is either not well-generated, or not a reflection group..." << endl;
	cout << "For now I need to enumerate all group elements, this take a very long time" << endl;
      }
      return false;
      }
    }
    else {
      cout << "Trouble studying projectivized reflection group" << endl;
      return false;
    }
  }
  else {
    // the two reflections commute

    /*    if (verbose) {
      cout << "Will now populate group generated by two commuting reflections" << endl;
      }*/
    
    vecteur A_powers;
    gen Aj = K.Id2;
    for (int j=0; j<ojmax; j++) {
      A_powers.push_back(Aj);
      Aj = Aj*Amax;
      K.matred(Aj);
    }

    vecteur B_powers;
    gen Bk = K.Id2;
    for (int k=0; k<okmax; k++) {
      B_powers.push_back(Bk);
      Bk = Bk*Bmax;
      K.matred(Bk);
    }

    vecteur obvious_list;
    for (int j=0; j<ojmax; j++) {
      for (int k=0; k<okmax; k++) {
	gen M = A_powers[j]*B_powers[k];
	K.matred(M);
	//	cout << "M=" << M << endl;
	obvious_list.push_back(M);
      }
    }
    
    /*    if (verbose) {
      cout << "... done, constructed " << obvious_list.size() << " elts" << endl;
      }*/
    
    bool are_elements_in_list = true;
    for (int j=1; j<elements.size(); j++) {
      bool fd = false;
      int k = 0;
      while (!fd && k<obvious_list.size()) {
	fd = K.isSameMatrix(obvious_list[k],elements[j]);
	k++;
      }
      if (!fd) {
	//	cout << "Elt #" << j << " is NOT in obvious list" << endl;
	cout << elements[j] << endl;
	are_elements_in_list = false;
      }
      /*      else {
	cout << "Elt #" << j << " is in obvious list" << endl;
	}*/
    }

    
    if (are_elements_in_list) {
      // we understand the group, it is generated by our two commuting reflections
	ord_gp = gmax;
	return true;
    }
    else {
      // could do more, maybe we want to give the list of missing generators?
      return false;
    }
    
  }

}