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

using namespace std;
using namespace giac;

FiniteGroup::FiniteGroup(CMfield k) {
  K = k;
  elements.push_back(K.Id);
  dim = 3;
  oldnumber=0;
};

FiniteGroup::FiniteGroup(CMfield k, int d) {
  K = k;
  gen M = K.create_identity(d);
  elements.push_back(M);
  dim = d;
  oldnumber=0;
};

bool FiniteGroup::areMultipleNum(gen M, gen N, gen eps) {
  // not sure why this is here, it is very clumsy code!
  
  gen Mn = K.matevaln(M);
  gen Nn = K.matevaln(N);
  gen A = Mn*K.inverse(Nn);
  //  cout << "A=" << A << endl;
  if (is_greater(abs(A[0][1],&ct),eps,&ct)) {
    //    cout << "01" << abs(A[0,1],&ct) << endl;
    return false;
  }
  else {
    if (is_greater(abs(A[0][2],&ct),eps,&ct)) {
      //      cout << "02" << abs(A[0,2],&ct) << endl;
      return false;
    }
    else {
      if (is_greater(abs(A[1][0],&ct),eps,&ct)) {
	//	cout << "10" << abs(A[1,0],&ct) << endl;
	return false;
      }
      else {
	if (is_greater(abs(A[1][2],&ct),eps,&ct)) {
	  //	  cout << "12" << abs(A[1,2],&ct) << endl;
	  return false;
	}
	else {
	  if (is_greater(abs(A[2][0],&ct),eps,&ct)) {
	    //	    cout << "20" << abs(A[2,0],&ct) << endl;
	    return false;
	  }
	  else {
	    if (is_greater(abs(A[2][1],&ct),eps,&ct)) {
	      //	      cout << "21" << abs(A[2,1],&ct) << endl;
	      return false;
	    }
	    else {
	      if (is_greater(abs(A[0][0]-A[1][1],&ct),eps,&ct)) {
		//		cout << "00-11" << abs(A[0,0]-A[1,1],&ct) << endl;
		return false;
	      }
	      else {
		if (is_greater(abs(A[1][1]-A[2][2],&ct),eps,&ct)) {
		  //		  cout << "11-22" << abs(A[1,1]-A[2,2],&ct) << endl;
		  return false;
		}
		else {
		  return true;
		}
	      }
	    }
	  }			      
	}
      }
    }
  }		      
}

bool FiniteGroup::areMultiple(gen M, gen N) {
  // starts by finding a nonzero element in M, then uses it as a "pivot"
  bool fdnz = false;
  int k=-1;
  int l;
  int dm1 = dim-1;
  while (!fdnz && k!=dm1) {
    k++;
    l=-1; 
    while (!fdnz && l!=dm1) {
      l++;
      fdnz = !operator_equal(M[k][l],K.zeronumber,&ct);
      if (!fdnz && !operator_equal(N[k][l],K.zeronumber,&ct)) {
	return false;
      }
    }
  }
  if (!fdnz) {
    cerr << "Matrix is zero, this should not happen in a group!" << endl;
    exit(1);
  }
  else {
    //    cout << "fdnz" << endl;
    //    cout << "k=" << k << ",  l=" << l << endl;
    gen a = M[k][l];
    gen b = N[k][l];
    bool ismult = true;
    while (ismult && l<dim) {
      gen te = b*M[k][l]-a*N[k][l];
      K.red(te);
      ismult = operator_equal(te,K.zeronumber,&ct);
      l++;
    }
    k++;
    while (ismult && k<dim) {
      l=0;
      while (ismult && l<dim) {
	//	cout << "k=" << k << ",  l=" << l << endl;
	gen te = b*M[k][l]-a*N[k][l];
	K.red(te);
	ismult = operator_equal(te,K.zeronumber,&ct);
	l++;
      }
      k++;
    }
    return ismult;
  }
}

/*
bool FiniteGroup::isNew(gen M) {
  return !isInList(M,elements);
}
*/

bool FiniteGroup::isInList(gen M, vector<gen> list) {
  bool res = false;
  int k=0;
  while (!res && k<list.size()) {
    res = areMultiple(M,list[k]);
    //    res = areMultipleNum(M,list[k],1e-6);
    k++;
  }
  return res;
}

void FiniteGroup::addElement(gen M) {
  if (!isInList(M,elements)) {
    //    cout << "Added new elt: " << M << endl;
    elements.push_back(M);
  }
}

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

void FiniteGroup::expand() {
  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 FiniteGroup::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 FiniteGroup::exploreCosets(vecteur &H, vecteur &transversal) {
  int jmax = 0;
  int omax = 1;
  for (int j=1; j<elements.size(); j++) {
    int oj = K.order(elements[j]);
    if (oj>omax) {
      jmax = j;
      omax = oj;
    }
  }

  //  vecteur H;
  H.push_back(elements[0]);
  gen M = elements[jmax];
  for (int j=1; j<omax; j++) {
      H.push_back(M);
      M = M*elements[jmax];
      K.matred(M);
  }
  // populated the subgroup generated by an element of largest order
  cout << "|H|=" << H.size() << endl;
  
  //  vecteur transversal;
  transversal.push_back(elements[0]);

  bool is_group = false;
  int ct = 0;

  vecteur new_transversal;
  for (int j=1; j<elements.size(); j++) {
    new_transversal.push_back(elements[j]);
  }

  while (!is_group && ct<10) {
    while (new_transversal.size()>0) {
      gen new_el = new_transversal.pop_back();
      bool fd = false;
      int k = 0;
      gen M = K.matinverse2d(new_el);
      while (!fd && k<transversal.size()) {
	gen M1 = M*transversal[k];
	K.matred(M1);
	int l = 0;
	while (!fd && l<omax) {
	  fd = areMultiple(M1,H[l]);
	  l++;
	}
	k++;
      }
      if (!fd) {
	transversal.push_back(new_el);
      }
    }
    cout << "Number of elements in transversal: " << transversal.size() << endl;
    bool te = true;
    int nt = transversal.size();
    for (int j=0; j<nt; j++) {
      if (j!=jmax) {
	gen Aj = transversal[j];
	for (int m=0; m<omax; m++) {
	  gen Am = Aj * H[m];
	  K.matred(Am);
	  for (int k=0; k<nt; k++) {
	    if (k!=jmax) {
	      gen Ak = Am * transversal[k];
	      K.matred(Ak);
	      // check if there exists an l such that gl^-1*A is in H
	      int l = 0;
	      bool tel = false;
	      while (!tel && l<nt) {
		int o = 0;
		bool is_in_H = false;
		while (!is_in_H and o<omax) {
		  gen A = transversal[l]*H[o];
		  K.matred(A);
		  is_in_H = areMultiple(Ak,A);
		  o++;
		}
		tel = is_in_H;
		l++;
	      }
	      if (!tel) {
		  te = false;
		  new_transversal.push_back(Ak);
	      }
	    }
	  }
	}
      }
    }
    is_group = te;
    ct++;
  }
  if (!is_group) {
    cout << "Failed constructing a group (after 10 iterations, perhaps you could try increasing that parameter)" << endl;
    return false;
  }
  else {
    //    cout << "Subgroup of order " << omax << " has " << transversal.size() << " cosets " << endl;
    return true;
  }
}