// This code correspond to the algorithm described in section 4 of the article
// "How far away must forced letters be so that squares are still avoidable?"
// It should be compiled with c++11 or higher and optimization O3 wich gives:
//    g++ lemma_7.cpp -O3 -lgmp-lgmpxx -o <executable_name>
// where <N> is the size of the alphabet
// that is, for testing the first statement of Theorem 9 one could do:
//      g++ lemma_7.cpp -O3 -DTEST3 -o testing_3.out
//      ./ testing_3.out
#include <iostream>
#include <vector>
#include <array>
#include <algorithm>
#include <bitset>
using namespace std;

#ifdef TEST3
  #define P 61
  #define sA 3
  #define N 320
  array<vector<int>, 28> W;// we keep the mirrors of the elements of W
  array<int, 28> f={4,19,19,19,22,22,22,28,28,28,36,36,36,50,50,50,63,63,63,88,88,88,118,118,118,148,148,148};
  void initW(){
    int t=0;
    W[t++]=vector<int>(9,-1);
    for(int i=19; i<=27;i++)
      for(int a=0; a<3; a++){
        W[t]=vector<int>(i,-1);
        W[t++][0]=a;
      }
  }
#elif defined TEST4
  #define P 18
  #define sA 4
  #define N 96
  array<vector<int>, 21> W;// we keep the mirrors of the elements of W
  array<int, 21> f={2,5,5,5,5,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8};
  void initW(){
    int t=0;
    W[t++] = vector<int>(1,-1);
    for(int a=0; a<sA; a++){
      W[t] = vector<int>(4,-1);
      W[t][1] = a;
      t++;
    }
    for(int a=0; a<sA; a++){
        for(int b=0; b<sA; b++){
        W[t] = vector<int>(6,-1);
        W[t][0]=a;
        W[t][3]=b;
        t++;
      }
    }
  }
#elif defined TEST6
  #define P 12
  #define sA 6
  #define N 96
  array<vector<int>, 43> W;// we keep the mirrors of the elements of W
  array<int, 43> f;
  void initW(){
    int t=0;
    W[t] = vector<int>(1,-1);
    f[t++] = 3;
    for(int a=0; a<sA; a++){
      W[t] = vector<int>(3,-1);
      W[t][1] = a;
      f[t++]=6;
    }
    for(int a=0; a<sA; a++){
        for(int b=0; b<sA; b++){
        W[t] = vector<int>(4,-1);
        W[t][0]=a;
        W[t][2]=b;
      f[t++]=6;
      }
    }
  }
#endif


class optstring:bitset<N>{
public:
  using bitset<N>::operator[];
  // [0]...[7] is the size of the string in binary
  // [8]...[N-1] is the mirror of the string
  // letter i of the mirror image is [8+2i]+2*[8+2i+1] if sA <= 4 or  [8+2i]+2*[8+2i+1]+4*[8+2i+1] if 4 < sA <= 8
  optstring():bitset<N>(){}
  void Cut(){// find g(f(u))
    int n = size();
    for(int i=0; i<n; i++){
      setsize(i+1);
      if(noNeedToComplete() || i+1 >= 2*P-3) return;
    }
  }
  
  void pop_back() {
    setsize(size()-1);
  }
  
  int size() const{// return the size
    int n=0;
    for(int i=0; i<7; i++){
      n += ((*this)[i]?1:0)<<i;
    }
    return n;
  }
  
  void setsize(int n){// set the size to a new value
    for(int i=0; i<7; i++)
      (*this)[i]=(n&(1<<i));
  }
  
#ifndef TEST6
  optstring(optstring& s, int a):bitset<N>(){// a copie of s with a at the end (so a starts the mirror image)
    int n = s.size();
    setsize(n+1);
    (*this)[8] = a%2;
    (*this)[9] = a/2;
    for(int i=8; i<8+2*n; i++)
      (*this)[i+2]=s[i];
  }
  
  int back() const{//return the last letter of the mirror (so the first letter of the word)
    int n=size()-1;
    return ((*this)[n*2+8]?1:0)+((*this)[n*2+9]?2:0);
  }
  
  optstring& operator += (int c){ //catenation of a letter at the end of the mirror image
    int n = size();
    (*this)[n*2+8]=(c%2);
    (*this)[n*2+9]=(c/2);
    setsize(n+1);
    return *this;
  }
  
  optstring& operator++(int){ //add 1 to the last letter
    int n=size()-1;
    (*this)[n*2+9] = (*this)[n*2+9] | (*this)[n*2+8];
    (*this)[n*2+8] = !(*this)[n*2+8];
    return *this;
  }
  
  friend bool operator< (const optstring& x, const optstring& y){
    //any kind of comparator is fine 
    //we need one to manipulate set in time O(log(n))
    int n= x.size();
    if(n<y.size()) return true;
    if(n>y.size()) return false;
    for (int i = 8; i < 8+2*n; i++) {
        if (x[i] && !y[i]) return false;
        if (!x[i] && y[i]) return true;
    }
    return false;
  }
  
  bool hasSquareLast()const{//return true if there is a square that contains the last letter
    int s = size();
    for(int p = 1; p*2<=s; p++){//p = period of the possible square
      bool isrep = true;
      for(int j=0; j<p; j++){
        if(((*this)[(s-1-j)*2+8]!= (*this)[(s-1-j-p)*2+8])||((*this)[(s-1-j)*2+9]!= (*this)[(s-1-j-p)*2+9])){
          isrep = false;
          break;
        }
      }
      if(isrep) return true;
    }
    return false;
  }
  
  bool hasSquareFirst()const{//return true if there is a square that contains the first letter
    int s = size();
    for(int p = 1; p*2<=s; p++){//p = period of the possible square
      bool isrep = true;
      for(int j=0; j<p; j++){
        if(((*this)[(j)*2+8]!= (*this)[(j+p)*2+8])||((*this)[(j)*2+9]!= (*this)[(j+p)*2+9])){
          isrep = false;
          break;
        }
      }
      if(isrep) return true;
    }
    return false;
  }
  
  bool noNeedToComplete(){
    //return true if f(*this)=*this
    //we apply the definition of f from the article
    size_t i=P-1, s=size();
    for(; i>=1; i--){
      bool possiblePeriod = true;
      for(size_t j=0; j+i<s; j++)
        if(((*this)[2*j+8]!= (*this)[2*(j+i)+8])||((*this)[2*j+9]!= (*this)[2*(j+i)+9])){ 
          possiblePeriod = false; 
          break;
        }
      if(possiblePeriod) break;
    }
    if(i == 0) return true;
    else return false;
  }
#else
  optstring(optstring& s, int a):bitset<N>(){// a copie of s with a at the end (so a starts the mirror image)
    int n = s.size();
    setsize(n+1);
    (*this)[8] = a%2;
    (*this)[9] = (a/2)%2;
    (*this)[10] = (a/4);
    for(int i=8; i<8+3*n; i++)
      (*this)[i+3]=s[i];
  }
  
  int back() const{//return the last letter of the mirror (so the first letter of the word)
    int n=size()-1;
    return ((*this)[n*3+8]?1:0)+((*this)[n*3+9]?2:0)+((*this)[n*3+10]?4:0);
  }
  
  optstring& operator += (int c){ //catenation of a letter at the end of the mirror image
    int n = size();
    (*this)[n*3+8]=(c%2);
    (*this)[n*3+9]=(c/2)%2;
    (*this)[n*3+10]=(c/4);
    setsize(n+1);
    return *this;
  }
  
  optstring& operator++(int){ //add 1 to the last letter
    int n=size()-1;
    (*this)[n*3+10] = (*this)[n*3+10] | ((*this)[n*3+8]&(*this)[n*3+9]);
    (*this)[n*3+9] = ((*this)[n*3+9] & (~(*this)[n*3+8])) | ((~(*this)[n*3+9]) & (*this)[n*3+8]);
    (*this)[n*3+8] = ~(*this)[n*3+8];
    return *this;
  }
  
  friend bool operator< (const optstring& x, const optstring& y){
    //any kind of comparator is fine 
    //we need one to manipulate set in time O(log(n))
    int n= x.size();
    if(n<y.size()) return true;
    if(n>y.size()) return false;
    for (int i = 8; i < 8+3*n; i++) {
        if (x[i] && !y[i]) return false;
        if (!x[i] && y[i]) return true;
    }
    return false;
  }
  
  bool hasSquareLast()const{//return true if there is a square that contains the last letter
    int s = size();
    for(int p = 1; p*2<=s; p++){//p = period of the possible square
      bool isrep = true;
      for(int j=0; j<p; j++){
        if(((*this)[(s-1-j)*3+8]!= (*this)[(s-1-j-p)*3+8])||((*this)[(s-1-j)*3+9]!= (*this)[(s-1-j-p)*3+9])||((*this)[(s-1-j)*3+10]!= (*this)[(s-1-j-p)*3+10])){
          isrep = false;
          break;
        }
      }
      if(isrep) return true;
    }
    return false;
  }
  
  bool hasSquareFirst()const{//return true if there is a square that contains the first letter
    int s = size();
    for(int p = 1; p*2<=s; p++){//p = period of the possible square
      bool isrep = true;
      for(int j=0; j<p; j++){
        if(((*this)[(j)*3+8]!= (*this)[(j+p)*3+8])||((*this)[(j)*3+9]!= (*this)[(j+p)*3+9])||((*this)[(j)*3+10]!= (*this)[(j+p)*3+10])){
          isrep = false;
          break;
        }
      }
      if(isrep) return true;
    }
    return false;
  }
  
  bool noNeedToComplete()const{
    //return true if f(*this)=*this
    //we apply the definition of f from the article
    size_t i=P-1, s=size();
    for(; i>=1; i--){
      bool possiblePeriod = true;
      for(int j=0; j+i<s; j++)
        if(((*this)[3*j+8]!= (*this)[3*(j+i)+8])||((*this)[3*j+9]!= (*this)[3*(j+i)+9])||((*this)[3*j+10]!= (*this)[3*(j+i)+10])){ 
          possiblePeriod = false; 
          break;
        }
      if(possiblePeriod) break;
    }
    if(i == 0) return true;
    else return false;
  }
#endif
};

void find(vector<optstring> & listwords){// find the set of non-empty square free words w such that g(f(w))=w
	optstring w; w+=0;
	while(true){
		if(!w.hasSquareLast()){
			if(w.size()<2*P-3 && !w.noNeedToComplete()){
				w += 0;
			}else{
 				listwords.push_back(w);
				while(w.size() > 0 && w.back() >= sA-1) w.pop_back();
				if (w.size() == 0) return;
				w++;
			}
		}else{
      while(w.size() > 0 && w.back()>= sA-1) w.pop_back();
      if (w.size() == 0) return;
      w++;
		}
	}
}

size_t getIndex(vector<optstring>& listwords, optstring& e){// dichotomic search of e in listwords
  size_t lb=0, ub=listwords.size()-1;   
  while(lb<ub){
    size_t mid = (lb+ub)/2;
    if(listwords[mid]<e) lb = mid+1;
    else ub = mid;
  }  
  return lb;
}

void getGraph(vector<optstring>& listwords, vector<array<long long, sA> >& graph){//compute the edges of the graph
	for(size_t i=0; i< listwords.size(); i++){
    if(i%2000000==0) cout<< i<<" vertices done out of "<< listwords.size()<<"  ( "<<100*double(i)/double(listwords.size())<< "% )"<<endl;
		for(int a=0; a < sA; a++){
      optstring wordtest(listwords[i], a);
      if(wordtest.hasSquareFirst()){ graph[i][a]=-1; continue;}
      else{
        wordtest.Cut();
        size_t t = getIndex(listwords, wordtest);
        graph[i][a] = t;
      }
    }
  }
}

void cleanGraph(vector<array<long long,sA> >& graph){//remove all the vertices with P_{i,d,G[X]}(v)<c
  vector<int> worstOutDeg(graph.size(),0);//contains P_{i,d,G[X]}
  vector<int> worstOutDegNew(graph.size(),0);
  vector<bool> valid(graph.size(),true);
  bool todo = true;
  while(todo){//if we removed vertices in the last round, we need to do at least one more round
    todo = false;
    // this is the dynamic algorithm computing P_{i,d,G[X]}
    for(size_t w=0; w<W.size(); w++){
      for(size_t i=0; i<graph.size();i++)
        if(!valid[i]) worstOutDeg[i]=0;
        else worstOutDeg[i]=1;
      for(size_t ff=0; ff<W[w].size(); ff++){// computes P_{ff,d,G[X]} from P_{ff-1,d,G[X]}
        for(size_t i=0; i<graph.size(); i++){
          worstOutDegNew[i]=0;
          if(W[w][ff] !=-1 ){
            if(graph[i][W[w][ff]] != -1 && valid[graph[i][W[w][ff]]]) worstOutDegNew[i] = worstOutDeg[graph[i][W[w][ff]]];//-m[forb];
          }else  
            for(size_t j=0; j<sA; j++) 
              if(graph[i][j] != -1 && valid[graph[i][j]]) worstOutDegNew[i] += worstOutDeg[graph[i][j]];
        }
        worstOutDegNew.swap(worstOutDeg); 
      }  
      for(size_t i=0; i<graph.size();i++) if(valid[i] && worstOutDeg[i]<f[w]) {valid[i] = false; todo = true;} 
    }
  }
  int left=0;
  for(size_t i=0; i<valid.size(); i++) if (valid[i]) left ++;
  cout<<"There are "<< left<<" vertices remaining in the sub-graph satisfying the conditions."<<endl;
}


int main(){
  initW();
  cout<<"Running with:  p= "<<P<<"."<<endl<< "Computation of the vertices of g(f(R_"<<P<<"))"<<endl;
	vector<optstring> listwords;
	find(listwords);
	cout <<"There are "<< listwords.size()<<" vertices in g(f(R_P))."<<endl<<"Sorting the vertices to allow set acces in O(log(n))."<<endl;
  sort(listwords.begin(), listwords.end());
  cout<< "Words sorted."<<endl<<"Starting graph generation."<<endl;
  vector<array<long long,sA>> graph(listwords.size());
  getGraph(listwords, graph);
  vector<optstring>().swap(listwords);
  cout<<"Graph found."<<endl<<"Finding the subgraph."<<endl;
 	cleanGraph(graph);
	return 0;
}