/*
  This code proves Theorem 5.3, i.e. the morphic word avoids additive square.
  It also proves Theorem 5.7, i.e. the morphic word avoids abelian squares of period at least 6
  Theorem 5.1 is thus proved as a corollary of Theorem 5.3 or 5.7.

  
  Code up to line 530 defines some linear algebra prerequisites:
   exact arithmetic in Q[sqrt(3)], vectors, matrices, Smith decomposition.
  
   
   Wrote in C++-11.

  Compilation: 
  Compiling with -DADDITIVE produces a program checking Theorem 5.3.
  Remove -DADDITIVE if you want to check Theorem 5.7.
  with exact arithmetics (needs gmpxx)
  g++ -Wall -DQSQRT3 -DADDITIVE -std=c++0x -O3 -o <file> <file>.cpp -lgmpxx -lgmp

  with double:
  g++ -Wall -DADDITIVE -std=c++0x -O3 -o <file> <file>.cpp

*/


#include <array>
#include <stdlib.h>
#include <assert.h> 
#include <iostream>
#include <vector>
#include <set>
#include <string>
#include <map>
#include <list>
#include <fstream>
#include <algorithm> 
#include <iomanip>
#include <sstream>  
#include <math.h>  
using namespace std;

/* If QSQRT3 is set, we make the computation in exact arithmetic (i.e. in Q[sqrt(3)]), otherwise the computation are made with doubles.*/
#ifdef QSQRT3

#include <gmpxx.h>

/* 
   A class for exact arithmetics in Q[sqrt(3)] 
*/
class Qsqrt3{ 
  mpq_class a,b; // a + b * sqrt(3)
	
public:
  Qsqrt3(int aa=0) : a(aa), b(0) {}
  Qsqrt3(const mpq_class &aa) : a(aa), b(0) {}
  Qsqrt3(const mpq_class &aa, const mpq_class &bb) : a(aa), b(bb) {}
  Qsqrt3(const mpq_class &aa, const mpq_class &bb, const mpq_class &cc) : a(aa/cc), b(bb/cc) {}

  explicit operator double() const {
    return a.get_d()+sqrt(3)*b.get_d();
  }

  friend ostream& operator<<(ostream& out, const Qsqrt3& x) {
    stringstream text;
    if(x.b ==0)text <<x.a ;
    else if(x.a ==0)text<< x.b <<"V(3)";
    else text<< x.a <<((x.b>0)?"+":"")<< x.b <<"V(3)";
    out<< text.str();
    return out;
  }
  Qsqrt3& operator+=(const Qsqrt3& rhs) {          
    a+= rhs.a;
    b+= rhs.b;
    return *this; 
  }
  friend Qsqrt3 operator+(Qsqrt3 lhs, const Qsqrt3& rhs) {
    lhs += rhs; 
    return lhs; 
  }	
  Qsqrt3& operator-=(const Qsqrt3& rhs) {          
    a-= rhs.a;
    b-= rhs.b;
    return *this; 
  }	
  Qsqrt3 operator-() const {
    return Qsqrt3(-a,-b);
  }	
  friend Qsqrt3 operator-(Qsqrt3 lhs, const Qsqrt3& rhs) {
    lhs -= rhs; 
    return lhs; 
  }	
  friend Qsqrt3 operator*(Qsqrt3 lhs, const Qsqrt3& rhs) {
    mpq_class temp = lhs.a*rhs.a+3*lhs.b*rhs.b;
    lhs.b = lhs.a*rhs.b+lhs.b*rhs.a; 
    lhs.a = temp;
    return lhs; 
  }	
  friend Qsqrt3 operator*(int lhs,const Qsqrt3& rhs) {	
    return Qsqrt3(lhs*rhs.a,lhs*rhs.b); 
  }	
  bool isnonneg()const{
    if(a<=0 && b<0) return false;
    if(a >=0 && b >= 0) return true;
    else if(a>0 && a*a >= 3 * b*b) return true;
    else if(a<0 && a*a <= 3 * b*b) return true;
    else return false;
  }
  friend Qsqrt3 operator/(Qsqrt3 lhs, const Qsqrt3& rhs) {
    assert( rhs.a!=0 || rhs.b!=0);
    mpq_class temp = lhs.a*rhs.a-3*lhs.b*rhs.b;
    mpq_class temp2 = rhs.a*rhs.a-3*rhs.b*rhs.b;
    lhs.b = lhs.b*rhs.a-lhs.a*rhs.b; 
    lhs.a = temp;
    lhs.a = lhs.a/temp2;
    lhs.b = lhs.b/temp2;
    return lhs; 
  }	
  friend inline bool operator< (const Qsqrt3& lhs, const Qsqrt3& rhs) {
    Qsqrt3 temp = lhs-rhs;
    return !(temp.isnonneg());		
  }
  friend inline bool operator> (const Qsqrt3& lhs, const Qsqrt3& rhs) {return rhs < lhs;}
  friend inline bool operator<=(const Qsqrt3& lhs, const Qsqrt3& rhs) {return !(lhs > rhs);}
  friend inline bool operator>=(const Qsqrt3& lhs, const Qsqrt3& rhs) {return !(lhs < rhs);}
  friend inline bool operator==(const Qsqrt3& lhs, const Qsqrt3& rhs) {return (lhs.a == rhs.a && lhs.b == rhs.b);}
  friend inline bool operator!=(const Qsqrt3& lhs, const Qsqrt3& rhs) {return (lhs.a != rhs.a || lhs.b != rhs.b);}	
  Qsqrt3 abs()const{
    if (! this->isnonneg()) return -*this;
    return *this;
  }	

  //return the integer (only if this really is an integer)
  int getint() const {
    assert(b==0 && a.get_den().get_si() == 1); // the programm abort if it not an integer
    return a.get_num().get_si();		
  }
};

typedef Qsqrt3 F_t;
#define F(a,b,c) F_t(a,b,c) //F(a,b,c) is a shortcut for the number (a+b*sqrt(3))/c

#define Fdouble(a) (double(a))
bool isNull(const Qsqrt3 &a) { return a==0;}
Qsqrt3 Fabs(const Qsqrt3 &a) { return a.abs();}
int Fint(const Qsqrt3 &a) { return a.getint();}

#else //QSQRT3
// otherwise, we make computations with doubles

typedef double F_t;

#define F(a,b,c) ((double(a)+sqrt(3)*(b))/c)
#define Fdouble(a) (double(a))
bool isNull(double a) { return fabs(a)<0.000000001;}
double Fabs(double a) {return fabs(a);}
int Fint(double a) { return int(round(a));}
#endif  //QSQRT3

bool isNull(int a) { return a==0;}
int Fabs(int a) { return abs(a);}

template<typename Tvalues,size_t n,size_t m> class matrix;
template<typename Tvalues,size_t n> class vect;

template<size_t n> void matinv(array<array<int,n>,n> &m);
template<class Tvalues,size_t n> void matinv(array<array<Tvalues,n>,n> &m);

/* 
   A class for n x m matrices of type Tvalues 
*/
template<typename Tvalues,size_t n,size_t m> 
class matrix : public array<array<Tvalues,m>,n> {
public:
  
  size_t cols() const {return m;}
  size_t rows() const {return n;}

  matrix() {
    for(unsigned int i=0; i<n; i++)
      for(unsigned int j=0; j<m; j++)
	(*this)[i][j]=0;
  }

  matrix(const vector<string> &morph) {
    for(unsigned int i=0; i<n; i++)
      for(unsigned int j=0; j<m; j++)
	(*this)[i][j]=0;
    for(unsigned int i=0; i<n; i++)
      for (unsigned int j=0; j< morph[i].size(); j++)
	(*this)[morph[i][j]-'a'][i]++;
  }

  Tvalues& operator()(int i, int j) { return (*this)[i][j];}
  Tvalues operator()(int i, int j)const { return (*this)[i][j];}
  friend ostream& operator<<(ostream& out, const matrix & x) {
    for(unsigned int i=0; i<x.rows(); i++) {
      for(unsigned int j=0; j<x.cols(); j++)
	out << setw(10)<< x[i][j]<< " ";
      out<<endl;
    }
    return out;
  }

  matrix & operator+=(const matrix & rhs) 
  {          
    for(unsigned int i=0; i<n; i++)
      for(unsigned int j=0; j<m; j++)
	(*this)[i][j]+= rhs[i][j];
    return *this; 
  }	
  matrix  operator+( const matrix & rhs)const
  {
    matrix val(*this);
    val += rhs; 
    return val; 
  }	
  
  template<class T2>
  matrix & operator-=(const matrix<T2,n,m>& rhs) 
  {       
    for(unsigned int i=0; i<n; i++)
      for(unsigned int j=0; j<m; j++)
	(*this)[i][j]-= rhs(i,j);
    return *this; 
  }		
      
  template<class T2>
  matrix  operator-( const matrix<T2,n,m>& rhs)const
  {
    matrix val(*this);
    val -= rhs; 
    return val; 
  }	

  template<typename T2,size_t t> 
  matrix<Tvalues,n,t>  operator* (const matrix<T2,m,t>& rhs)const
  {
    matrix<Tvalues,n,t> val;
    for(unsigned int i=0; i<n; i++)
      for(unsigned int j=0; j<rhs.cols(); j++)
	for(unsigned int k=0; k<m; k++) 
	  val(i,j)+=(*this)[i][k]*rhs(k,j);
    return val; 
  }		

  template<typename T2> 
  friend matrix operator*(const Tvalues& lhs, const matrix<T2,n,m>& rhs)
  {
    matrix val;
    for(unsigned int i=0; i<rhs.rows(); i++)
      for(unsigned int j=0; j<rhs.cols(); j++)
	val(i,j)=lhs*rhs(i,j);
    return val; 
  }

  friend matrix operator/(const matrix & lhs, const Tvalues& rhs)
  {
    matrix  val;
    for(unsigned int i=0; i<lhs.n; i++)
      for(unsigned int j=0; j<lhs.m; j++)
	val(i,j)=lhs(i,j)/rhs;
    return val; 
  }	

  matrix<Tvalues,m,n>  transp() {
    matrix<Tvalues,m,n> trans;
    for(unsigned int i=0; i<n; i++)
      for(unsigned int j=0; j<m; j++)
	trans(j,i)=(*this)[i][j];
    return trans;
  }

  void rowswap(int i, int j) { (*this)[i].swap((*this)[j]); }
  void colswap(int i, int j) {
    for(unsigned int k=0; k<n; k++) {
      Tvalues tmp= (*this)[k][i];
      (*this)[k][i]= (*this)[k][j];
      (*this)[k][j]=tmp;
    }
  }
  matrix  inv() const {
    auto r=*this;
    matinv<>(r);
    return r;
  }
  
  matrix  abs()const
  {
    matrix  val;
    for(unsigned int i=0; i<n; i++)
      for(unsigned int j=0; j<m; j++)
	val(i,j)= Fabs((*this)[i][j]);
    return val;
  }

	
  Tvalues squaredNorm() {
    Tvalues val=0;
    for(unsigned int i=0; i<n; i++)
      for(unsigned int j=0; j<m; j++)
	val+= (*this)[i][j]*(*this)[i][j];
    return val;
  }

  // return an std vector of vect corresponding to the column
  vector<vect<Tvalues,n> >getstdvector()const {
    vector<vect<Tvalues,n> > ret(m,vect<Tvalues,n>());
    for(unsigned int i=0; i<m; i++)
      for(unsigned int j=0; j<n;j++)
	ret[i](j)= (*this)[j][i];	
    return ret;
  }


  bool isId() const {
    if(n!=m) return false;
    for(unsigned int i=0;i<n;i++)
      for(unsigned int j=0;j<n;j++)
	if(!::isNull((*this)[i][j]-int(i==j))) return false;
    return true;
  }
  
  bool isNull() const {
    for(unsigned int i=0;i<n;i++)
      for(unsigned int j=0;j<m;j++)
	if(!::isNull((*this)[i][j])) return false;
    return true;
  }
  
};

/* 
   mathinv(mat) computes the inverse of the matrix mat 
*/

/* The inverse of an unimodular integer matrix is an integer matrix
   but the intermediate computation have to be donne in Q. */
template<size_t n> 
void matinv(array<array<int,n>,n> &mat) {
  matrix<F_t,n,n> tmpMat;
  for(unsigned int i=0; i<n; i++)
    for(unsigned int j=0; j<n; j++)
      tmpMat(i,j)=mat[i][j]; 
  matrix<F_t,n,n> tmpinvMat(tmpMat.inv());
  matrix<int,n,n> invMat;
  for(unsigned int i=0; i<n; i++)
    for(unsigned int j=0; j<n; j++)
      mat[i][j]=Fint(tmpinvMat(i,j));
}

template<class Tvalues,size_t n> 
void matinv(array<array<Tvalues,n>,n> &mat) {
  matrix<Tvalues,n,n> invMat;
  for(unsigned int i=0; i<n; i++) invMat(i,i)=1;
  auto mattemp(mat);
  for(unsigned int i=0; i<n; i++) {
    int minind=i;		
    bool isnul = true;
    for(unsigned int j=i; j<n; j++)
      if(mattemp[j][i] != 0) {minind =j; isnul = false;}
    assert(!isnul);
    mattemp[i].swap(mattemp[minind]);
    invMat.rowswap(i,minind);
    Tvalues coef = mattemp[i][i];
    for(unsigned int k=0; k<n; k++) {
      mattemp[i][k] =mattemp[i][k]/coef;
      invMat(i,k) = invMat(i,k)/coef;
    }
    for(unsigned int j=0; j<n; j++) {
      if(j==i) continue;
      coef = mattemp[j][i];
      for(unsigned int k=0; k<n; k++) {
	mattemp[j][k] -= coef*mattemp[i][k];
	invMat(j,k) -= coef*invMat(i,k);
      }
    }
  }
  for(unsigned int i=0;i<n;i++)
    for(unsigned int j=0;j<n;j++)
      mat[i][j]=invMat(i,j);
}

/*  
    A class for vectors of type Tvalues
*/
template<typename Tvalues,size_t n> 
class vect: public array<Tvalues,n> {
public:
  vect() {
    for(unsigned int i=0; i<n; i++) {
      (*this)[i]=0;
    }
  }
  const Tvalues &operator()(int i) const { return (*this)[i];}
  Tvalues &operator()(int i) { return (*this)[i];}
  friend ostream& operator<<(ostream& out, const vect & x) {
    for(unsigned int i=0; i<n; i++)
      out << setw(10)<< x(i)<< " ";
    return out;
  }	
	
//   template<typename T2,size_t m> 
//   friend vect operator*(const matrix<Tvalues,n,m>& lhs, const vect<T2,m>& rhs) {
//     vect<Tvalues,n> val;
//     for(unsigned int i=0; i<n; i++)
//       for(unsigned int k=0; k<m; k++) 
// 	val[i]+=lhs[i][k]*rhs[k];
//     return val; 
//   }

  const vect &operator+=(const vect & rhs) {
    for(unsigned int i=0;i<n;i++)
      (*this)[i]+=rhs[i];
    return *this;
  }
  const vect &operator-=(const vect & rhs) {
    for(unsigned int i=0;i<n;i++)
      (*this)[i]-=rhs[i];
    return *this;
  }
  vect operator+(const vect & rhs)const {
    vect r=*this;
    for(unsigned int i=0;i<n;i++)
      r[i]+=rhs[i];
    return r;
  }

  vect operator-(const vect & rhs)const {
    vect r=*this;
    for(unsigned int i=0;i<n;i++)
      r[i]-=rhs[i];
    return r;
  }

  Tvalues dot(const vect & rhs) {
    Tvalues res(0);
    for(unsigned int i=0; i<n; i++) res += (*this)[i]*rhs(i);
    return res;
  }	
  
  vect abs() const{
    vect r;
    for(unsigned int i=0;i<n;i++)
      r[i]=Fabs((*this)[i]);
    return r;
  }
  friend vect operator*(const Tvalues& lhs, const vect& rhs) {
    vect r;
    for(unsigned int i=0;i<n;i++)
      r[i]=lhs*rhs[i];
    return r;
  }
  friend vect operator/(const vect& rhs,const Tvalues& lhs) {
    vect r;
    for(unsigned int i=0;i<n;i++)
      r[i]=rhs[i]/lhs;
    return r;
  }
  Tvalues squaredNorm() {
    Tvalues val=0;
    for(unsigned int i=0; i<n; i++)
      val+= (*this)[i]*(*this)[i];
    return val;
  }
};

template<typename Tvalues,typename T2,size_t n,size_t m>
vect<Tvalues,n> operator*(const matrix<Tvalues,n,m>& lhs, const vect<T2,m>& rhs) {
  vect<Tvalues,n> val;
  for(unsigned int i=0; i<n; i++)
    for(unsigned int k=0; k<m; k++)
      val[i]+=lhs[i][k]*rhs[k];
  return val;
}

/* 
   smith(M,U,D,T) computes the Smith normal form of the integer matrix M
   output: integer matrices U,D,V such that UDV=M, D is diagonal and 
    U and T are invertible over integers (i.e. their determinant is 1 or -1).
*/
template<typename T,size_t n> 
void smith(const matrix<T,n,n>& M,matrix<T,n,n>& U, matrix<T,n,n>& D, matrix<T,n,n>& V) {
  for(unsigned int i=0;i<n;i++)
    for(unsigned int j=0;j<n;j++)
      U(i,j)=V(i,j)=(i==j);
  D=M;

  unsigned int colnul=0;
  for(unsigned int i=0; i<D.cols()-colnul; i++) {
    unsigned int minind=i;
    bool isnul = true;
    while(isnul && i<D.cols()-colnul) {
      for(unsigned int j=i; j<D.cols(); j++) if(D(j,i) != 0) { minind =j; isnul = false;}
      if(!isnul) break;
      else{
	colnul++;
	D.colswap(i,D.cols()-colnul);
	V.rowswap(i,D.cols()-colnul);
      }	
    }
    if(isnul) break;
    do{
      //line
      for(unsigned int j=i; j<D.cols(); j++) if(abs(D(j,i))< abs(D(minind,i)) && D(j,i) != 0) minind = j;	
      do{
	D.rowswap(i,minind);
	U.colswap(i,minind);
	for(unsigned int j=i+1; j<D.cols(); j++) {
	  int coef = D(j,i)/D(i,i);
	  for(unsigned int k=0; k<D.rows(); k++) D(j,k)-=coef*D(i,k);
	  for(unsigned int k=0; k<D.rows(); k++) U(k,i)+=coef*U(k,j);
	}
	minind=i;
	for(unsigned int j=i; j<D.cols(); j++) if(abs(D(j,i))< abs(D(minind,i)) && D(j,i) != 0) minind = j;
      }while(minind!=i);
			
      //column
      minind=i;
      for(unsigned int j=i; j<D.cols(); j++) if(abs(D(i,j))< abs(D(i,minind)) && D(i,j) != 0) minind = j;
      do{

	D.colswap(i, minind);
	V.rowswap(i, minind);
	for(unsigned int j=i+1; j<D.cols(); j++) {
	  int coef = D(i,j)/D(i,i);
	  for(unsigned int k=0; k<D.rows(); k++) D(k,j)-=coef*D(k,i);
	  for(unsigned int k=0; k<D.rows(); k++) V(i,k)+=coef*V(j,k);
	}
	minind=i;
	for(unsigned int j=i; j<D.cols(); j++) if(abs(D(i,j))<abs(D(i,minind)) && D(i,j) != 0) minind = j;
      } while(minind!=i);	
      for(unsigned int j=i+1; j<D.cols(); j++) if(D(j,i) != 0) minind = j;
    } while(minind!=i);
  }

}

/* Prints a morphism */
void printm(const vector<string>& m) {
  for (unsigned int i=0; i<m.size(); i++)
    cout << (char)(i+'a')<<"->"<<m[i]<<endl;  
}


#define nbletters 6 // we have 6 letters

/*
  A Class for a template (but this is called Pattern because c++ already uses the keyword template)
*/
class Pattern{
public:
  vector<int> c;
  vect<int,nbletters> u;
	
  Pattern():c(3,-1) {  }
  friend ostream& operator<<(ostream& out,Pattern& p) {
    out << p.c[0]<<" (0,0,0,0) " << p.c[1] << " (" << p.u << ") "<< p.c[2];
    return out;
  }
  inline bool operator <(const Pattern &P) const {
    return c<P.c || (c==P.c &&  u<P.u);
  }	
  inline bool operator == (const Pattern &P) const {
    return c==P.c && u==P.u;
  }
};

/*
  A Class to test the condition |r_i*.x|<bound
*/
class testor{
  vect<F_t,nbletters> x;
public:
  F_t bound;
  testor(F_t& boundn, vect<F_t,nbletters>& x): x(x), bound(boundn) {};
	
  template<typename T> 
  inline bool operator () (const vect<T,nbletters>& parikh) const {
    F_t r=0;
    for(int i =0; i<nbletters; i++)	
      r+= parikh(i)*x(i);
    return Fabs(r)<=bound;
  }	
};

/* 
   A class of one "cut". 
   This corresponds to a left vector x, a right vector y and a letter b
   such that h(b)= ucv; parikh(u)=x; parikh(v)=y and this cut belongs to c
   (it will be on the vector at position c, the function init from the main does that).
*/
class cutvec{
public:
  int letter;
  vect<int,nbletters> prevec;
  vect<int,nbletters> lastvec;
	
  cutvec(int n):letter(-1) {}
	
  friend ostream& operator<<(ostream& out, cutvec& p) {
    out << "letter = "<< p.letter<< "( ";
    for (int i=0; i<nbletters; i++)
      out << p.prevec(i)<< ",";
    out<< ") , (";
    for (int i=0; i<nbletters; i++)
      out << p.lastvec(i)<< ",";
    out<< ")"<<endl;
    return out;
  }
};

/* 
   The following function find the pair (s,p) from suff,pref that maximizes the dot product of evec and the parikh vector of sp.
   This allows the computation of the r_i*.
*/
F_t findrstarbegplusend(const vector<string>& morph, vect<F_t,nbletters>& eVec) {
  F_t ret=0;
  for(unsigned int i= 0; i<morph.size(); i++) {
    for(unsigned int ii= 0; ii<morph.size(); ii++) {
      F_t temp = 0;
      for(unsigned int j=0; j<morph[i].size()-1; j++) {
	temp += eVec(morph[i][j]-'a');
	F_t temp2 = 0;
	for(unsigned int k=morph[ii].size()-1; k>0; k--) {
	  temp2 += eVec(morph[ii][k]-'a');
	  if(Fabs(temp+temp2)>ret) ret = Fabs(temp+temp2);
	}
      }
    }
  }
  return ret;
}

/*
  Find the factor that maximizes the dot product of evec and the parikh vector of this factor.
  This allows to check the bounds on the r_i*.
*/
F_t findrstarfactor(const vector<string>& morph, vect<F_t,nbletters>& eVec) {
  F_t ret=0;
  F_t temp;
  string word=morph[0];
  for(int l=1; l<4; l++)
    for(unsigned int i=0; i<morph[word[l]-'a'].size(); i++)
      word.push_back(morph[word[l]-'a'][i]);
  for(unsigned int j=0; j<word.size()-2; j++) {
    temp = 0;
    for(unsigned int k=j; k<word.size()-1; k++) {
      temp += eVec(word[k]-'a');
      if(Fabs(temp)>ret) ret = Fabs(temp);
    }
  }
  return ret;
}

/*
  Initialize listOfCut: this is a list of the ways a letter can be obtained in the image of an other letter.
  This is usefull to compute it once so we do not have to compute it every time we want to have the parent of a template.
*/
void init(vector<cutvec>* listOfCut, const vector<string>& morph, unsigned int nb_letters) {
  cutvec c0h(nb_letters);
  listOfCut[0].push_back(c0h);
  for(unsigned int i=0; i<nb_letters; i++) {	
    cutvec c1(nb_letters);
    c1.letter=i;
    for (unsigned int j=0; j< morph[i].size(); j++)
      c1.lastvec(morph[i][j]-'a')++;
    for (unsigned int j=0; j< morph[i].size(); j++) {
      listOfCut[0].push_back(c1);
      c1.lastvec(morph[i][j]-'a')--;
      c1.prevec(morph[i][j]-'a')++;	
    }
    listOfCut[0].push_back(c1);
  }
	
  for(unsigned int i=0; i<nb_letters; i++) {	
    cutvec p0(nb_letters);
    p0.letter=i;
    for (unsigned int j=1; j< morph[i].size(); j++)
      p0.lastvec(morph[i][j]-'a')++;
    for (unsigned int j=0; j< morph[i].size()-1; j++) {
      listOfCut[morph[i][j]-'a'+1].push_back(p0);
      p0.lastvec(morph[i][j+1]-'a')--;
      p0.prevec(morph[i][j]-'a')++;	
    }
    listOfCut[morph[i][morph[i].size()-1]-'a'+1].push_back(p0);
  }
}


map<vect<int,nbletters>,list<vect<int,nbletters> > > alreadycomputedpreimages; 

set<Pattern> seenPattern; //this set allows us to add each pattern only once

/*
  The following function computes the list of parents of one pattern.
*/
template<int nb2>
void getParents(Pattern p, vector<Pattern>& listOfParents, vector<cutvec>*  pOfLet, vector<testor>& testlist,
		const matrix<int,nbletters,nbletters>& invU, const matrix<int,nbletters,nbletters>& D, const matrix<int,nbletters,nbletters>& invV,
		const vector<vect<int,nbletters> >& baseoftheker, F_t mumin, const matrix<F_t,nb2,nbletters>& proj, const vect<F_t,nb2>& projbound, const vector<vect<int,nbletters> >& kerforimproovement) {
  //pOfLet[-1] for epsilon
  //unsigned int nbletters = p.nbletters;
  Pattern newPat;
  vect<int,nbletters> tempu;
  vect<int,nbletters> uparents;
  vector<int> combin(baseoftheker.size(),0);
  vector<cutvec>::iterator cut[3];
	
  //we try all the preimages possible for a letter
  int nb_iter = pOfLet[p.c[0]].size() * pOfLet[p.c[1]].size() *pOfLet[p.c[2]].size();
  int compt =0;
  for(cut[0] = pOfLet[p.c[0]].begin(); cut[0] != pOfLet[p.c[0]].end(); cut[0]++)
    for(cut[1] = pOfLet[p.c[1]].begin() ; cut[1]!= pOfLet[p.c[1]].end(); cut[1]++)
      for(cut[2] = pOfLet[p.c[2]].begin() ; cut[2]!= pOfLet[p.c[2]].end(); cut[2]++) {
	if(compt%1000==0 && nb_iter>20000) cerr <<compt<< " out of "<<nb_iter <<" iterations. "<<"\r";
	compt ++;
	//We compute any possible vector
	//The matrix is not invertible so this is not completely trivial
	tempu = p.u - cut[1]->lastvec - cut[2]->prevec + cut[0]->lastvec + cut[1]->prevec;
	for(int i =0; i<3; i++)	newPat.c[i]= cut[i]->letter;
	bool ok = true;
	uparents = invU*tempu;
	for(unsigned int i=0; i<nbletters; i++) {
	  if(D(i,i)==0 && uparents(i) != 0) {ok = false; break;} // no preimage
	  else if(D(i,i)==0) continue;
	  else if(uparents(i) % D(i,i)!= 0) {ok = false; break;} // no integer preimage
	  else uparents(i) = uparents(i)/ D(i,i);
	}
	if (!ok) continue;
	tempu = invV*uparents;
	//If there is at least an integer solution we now have a particular solution
	//We try to take it as centered as possible in the ball so that we get better bounds to find all the other solutions
	//We only add and remove linear combination of some integers elements of the kernel so that we still have a particular solution
	F_t val;
	do{
	  val = (proj*tempu).squaredNorm();
	  for (unsigned int i=0; i<kerforimproovement.size(); i++) {
	    vect<F_t,nb2> PK = proj * kerforimproovement[i];	
	    if(PK.squaredNorm()==0) continue;
	    F_t coef = (proj*tempu).dot(PK)/PK.squaredNorm();
	    tempu = tempu - int(round(double(coef)))* kerforimproovement[i];
	  }
	}while(val*F(99,0,100) > (proj*tempu).squaredNorm()); 

	//This computes the bounds for the ball as described in the part "computation of the parents"
	F_t boundSquared =(projbound+(proj*tempu).abs()).squaredNorm()/mumin;
	int  ibound = floor(sqrt(double(boundSquared)));
	if(alreadycomputedpreimages.find(tempu)==alreadycomputedpreimages.end() ) {
	  auto &tt(alreadycomputedpreimages[tempu]);
	  int kk=0,ll=0;
	  for(unsigned int i=0; i<combin.size(); i++) combin[i]=-ibound;
	  while(true) {
	    kk++;
	    ok = true;
	    newPat.u=tempu;
	    for(unsigned int i=0;i<baseoftheker.size();i++)
	      newPat.u+=combin[i]*baseoftheker[i];
	    for(unsigned int i=0; i< testlist.size(); i++)
	      if(!testlist[i](newPat.u) ) { ok = false; break;}
	    if(ok) {
	      ll++;
	      tt.push_back(newPat.u);
	      if(seenPattern.insert(newPat).second)
		listOfParents.push_back(newPat);
	    }
	    unsigned int l=0;
	    while(l < combin.size() && combin[l]>=ibound) {
	      combin[l]=-ibound;
	      l++;
	    }
	    if(l >= combin.size()) break;
	    (combin[l])++;
	  }
	} else {
	  auto &tt(alreadycomputedpreimages[tempu]);
	  for(auto it=tt.begin();it!=tt.end();++it) {
	    newPat.u=*it;
	    if(seenPattern.insert(newPat).second)
	      listOfParents.push_back(newPat);
	  }
	}
				
      }
}

/* Computes the image of 'a' by the morphism morph. */
string image(const vector<string>& morph,const string &a) {
  string word;
  for(unsigned int j=0; j<a.size(); j++)
    word+=morph[a[j]-'a'];
  return word;
}

/* Fills factors with all the factors of size n in h(preimage) */
void computefactorsofsizenfromsizep(const vector<string>& morph, int n, vector<string> preimage, vector<string> & factors, 
				    set<string>& seenwords) {
  for(unsigned int i=0; i<preimage.size(); i++) {
    string word=image(morph,preimage[i]);

    for(unsigned int j=0; j<word.size()-n+1; j++) {
      string tempword(word.begin()+j,word.begin()+j+n);
      if(seenwords.insert(tempword).second) factors.push_back(tempword);
    }
  }	
}

/* Finds the size of the longest abelian square ending at position s and add them to 'squares' */
int maxSizeOfAbelianSquareEndingAts(const string& w, int s,set<string> &squares) {
  int max =0;
  vect<int,3> diff;
  for(int i = 1; i*2<=s; i++) {		
    diff(w[s-i]-'a')+=2;
    diff(w[s-2*i+1]-'a')--;
    diff(w[s-2*i]-'a')--;
    if(diff.squaredNorm() == 0) {
      max = i;
      string sq(w,s-2*i,2*i);
      squares.insert(sq);
    }
  }
  return max;
}

int main(int argc, char *argv[]) {
  
  //Initialization : everything is hard-coded for the morphic word of Theorem 6

  vector<string> firstmorphism(nbletters);
  firstmorphism[0]="ace";
  firstmorphism[1]="adf";
  firstmorphism[2]="bdf";
  firstmorphism[3]="bdc";
  firstmorphism[4]="afe";
  firstmorphism[5]="bce";

  matrix<int,6,6> M(firstmorphism); // The matrix M_h for the firstmorphism

  matrix<F_t,6,6> J; //the Jordan normal form of M
  J(0,1)=1;
  J(3,3)=3;
  J(4,4)=F(0,1,1);
  J(5,5)=F(0,-1,1);
	
  matrix<F_t,6,6>  P; 
  P(0,0)=F(-1,0,2); P(0,1)=0; 		P(0,2)=-1; P(0,3)=1; P(0,4)=F(2,1,1);   P(0,5)=F(2,-1,1);
  P(1,0)=F(1,0,2);  P(1,1)=-1;  P(1,2)=0;  P(1,3)=1; P(1,4)=F(-2,-1,1); P(1,5)=F(-2,1,1);
  P(2,0)=F(-1,0,2); P(2,1)=1; 		P(2,2)=-1; P(2,3)=1; P(2,4)=-1; 	P(2,5)=-1;
  P(3,0)=0; 	    P(3,1)=0; 		P(3,2)=1;  P(3,3)=1; P(3,4)=F(-3,-2,1); P(3,5)=F(-3,2,1);
  P(4,0)=0; 	    P(4,1)=F(1,0,2); 	P(4,2)=1;  P(4,3)=1; P(4,4)=F(3,2,1);   P(4,5)=F(3,-2,1);
  P(5,0)=F(1,0,2);  P(5,1)=F(-1,0,2);   P(5,2)=0;  P(5,3)=1; P(5,4)=1; 		P(5,5)=1;

  matrix<F_t,6,6>  Pinv = P.inv();

  assert((P*Pinv).isId()); // check that the Pinv is the inverse of P
  assert((P*J*Pinv-M).isNull()); // check that P*J*P^-1=M (that is we have a Jordan decomposition of M)

  matrix<int,6,2> baseofKerofM; // A base of the lattice of the integer points of the kernel of M
  baseofKerofM(0,0)= -1; baseofKerofM(1,0)=1; baseofKerofM(2,0)= -1; baseofKerofM(3,0)= 0; baseofKerofM(4,0)= 0; baseofKerofM(5,0)= 1;
  baseofKerofM(0,1)= -1; baseofKerofM(1,1)=0; baseofKerofM(2,1)= -1; baseofKerofM(3,1)= 1; baseofKerofM(4,1)= 1; baseofKerofM(5,1)= 0;	

  
#ifndef ADDITIVE
  vector<string> h(nbletters);
  h[0]="bbbaabaaac";
  h[1]="bccacccbcc";
  h[2]="ccccbbbcbc";
  h[3]="ccccccccaa";
  h[4]="bbbbbcabaa";
  h[5]="aaaaaaabaa";

  matrix<int,6,6> Mh(h); // The matrix M_h for the second morphism
		
  matrix<int,6,3> baseofKerofh;// A base of the lattice of the integer points of the kernel of Mh
  baseofKerofh(0,0)= 4; baseofKerofh(1,0)=-5; baseofKerofh(2,0)= 6;baseofKerofh(3,0)= 0;baseofKerofh(4,0)= -5;baseofKerofh(5,0)= 0;
  baseofKerofh(0,1)= 3; baseofKerofh(1,1)=3; baseofKerofh(2,1)= -4;baseofKerofh(3,1)= 0;baseofKerofh(4,1)= 0;baseofKerofh(5,1)= -2;
  baseofKerofh(0,2)= 0; baseofKerofh(1,2)=-2; baseofKerofh(2,2)= 1;baseofKerofh(3,2)= 1;baseofKerofh(4,2)= 0;baseofKerofh(5,2)= 0;

#else
  
  matrix<int,6,6> Mh;
  Mh(0,0)=1;Mh(1,0)=0;Mh(2,0)=0;Mh(3,0)=0;Mh(4,0)=0;Mh(5,0)=0;
  Mh(0,1)=1;Mh(1,1)=1;Mh(2,1)=1;Mh(3,1)=0;Mh(4,1)=0;Mh(5,1)=0;
  Mh(0,2)=1;Mh(1,2)=2;Mh(2,2)=1;Mh(3,2)=0;Mh(4,2)=0;Mh(5,2)=0;
  Mh(0,3)=1;Mh(1,3)=0;Mh(2,3)=1;Mh(3,3)=0;Mh(4,3)=0;Mh(5,3)=0;
  Mh(0,4)=1;Mh(1,4)=2;Mh(2,4)=0;Mh(3,4)=0;Mh(4,4)=0;Mh(5,4)=0;
  Mh(0,5)=1;Mh(1,5)=1;Mh(2,5)=0;Mh(3,5)=0;Mh(4,5)=0;Mh(5,5)=0;
  matrix<int,6,3> baseofKerofh;// A base of the lattice of the integer points of the kernel of Mh

  baseofKerofh(0,0)= 0; baseofKerofh(1,0)=2; baseofKerofh(2,0)= -1;baseofKerofh(3,0)= -1;baseofKerofh(4,0)= 0;baseofKerofh(5,0)= 0;
  baseofKerofh(0,1)= -1; baseofKerofh(1,1)=0; baseofKerofh(2,1)= 0;baseofKerofh(3,1)= 0;baseofKerofh(4,1)= -1;baseofKerofh(5,1)= 2;
  baseofKerofh(0,2)= 1; baseofKerofh(1,2)=1; baseofKerofh(2,2)= 0;baseofKerofh(3,2)= -1;baseofKerofh(4,2)= 0;baseofKerofh(5,2)= -1;
  
  
#endif
  cout << "*************** Initialization of the first morphism ***************"<<endl;
  cout << "The first morphism is: "<<endl;
  printm(firstmorphism);
  cout << "The matrix J is:"<<endl<< J <<endl<<endl;
  cout << "The matrix P^-1 is:"<<endl<<Pinv<<endl<<endl;
  cout << "The matrix M associated to the morphism (= P * J * P^-1) is: "<<endl << M<<endl;
  cout << "and P * J * P^-1 is: " <<endl<<  P*J*Pinv<<endl;
  cout << "A base of the kernel of M is:"<<endl<<baseofKerofM<<endl;
  cout << "Check M* the base ="<<endl<< M*baseofKerofM<<endl;
  assert((M*baseofKerofM).isNull()); // Check that the product is null

  cout << "A Jordan decomposition of M^2 is P * J^2 * P^-1."<<endl; //(Note: this is not true in general, but is this particular case it is.)
  auto J2=J*J;
  assert((P*J2*Pinv-M*M).isNull());
  vector<string> firstmorphism2(nbletters); // The composition of firstmorphism by itself
  for(unsigned int i=0;i<nbletters;i++)
    firstmorphism2[i]=image(firstmorphism,firstmorphism[i]);

  assert((matrix<int,6,6>(firstmorphism2)-M*M).isNull());

	
#ifndef ADDITIVE
  cout << "*************** Initialization of the second morphism ***************"<<endl;	
  cout << "The second morphism h is: "<<endl;
  printm(h);
  cout << "The matrix of this morphism is:"<<endl<< Mh<<endl;
  cout << "A base of the kernel of h is:"<<endl<< baseofKerofh <<endl;
  cout << "Check Mh* the base ="<<endl<< Mh* baseofKerofh<<endl;
  assert((Mh* baseofKerofh).isNull()); // Check that the product is null	
#else
    cout << "*************** Initialization of the second morphism ***************"<<endl;  
  cout << "The matrix of the morphism is:"<<endl<< Mh<<endl;
  cout << "A base of the kernel of h is:"<<endl<< baseofKerofh <<endl;
  cout << "Check Mh* the base ="<<endl<< Mh* baseofKerofh<<endl;
  assert((Mh* baseofKerofh).isNull()); // Check that the product is null
#endif
  /////////////////////////////////////////////////////////
  //////////// Computing bounds on the r_i* ///////////////
  /////////////////////////////////////////////////////////

  cout << "*************** Computing the r_i* ***************"<<endl;
  vector<testor> listtest; // we keep here the r_i and r_i* of the proposition 6
  F_t last=0;

  /* This matrix contains the subvectors of P on which we can project the kernel.
     This corresponds to the matrix Q of the "Computation of the parents" in the article. */
  matrix<F_t,2, nbletters> constrkerofmorphism;
  vect<F_t,2> boundkerofmorphism; // this contains the corresponding bounds
  int nbconstrkerofmorphism=0;

  /* This matrix contains all the subevtors of P corresponding to eigein-values of norm <1.
     This correspond to the matrix Q from the proposition 10 in the article.
     This is usefull later to compute the parents by the second morphism */
  matrix<F_t,3,nbletters> constrVect;
  vect<F_t,3> constr; // this contains the corresponding bounds
  int nbconstr = 0;
	
  /* We compute the r_i* as described on the article.
     We use the square of the morphism to get better bounds. */
  for(int i=nbletters-1; i>=0; i--) {
    if(Fabs(J(i,i))<1) {// We need to take care of the eigen-values of norm less than 1
      vect<F_t,nbletters> eVec;
      for(unsigned int j=0; j<nbletters; j++)eVec(j)= Pinv(i,j); // this is the vector corresponding to the projection r_i
	
      // This is the ith coordinate of the vector u in the proposition 3.3
      F_t lambda = 2*findrstarbegplusend(firstmorphism2, eVec)/(1-Fabs(J(i,i))); 
      if(i<(int)nbletters-1) lambda += J2(i,i+1)*last/(1-Fabs(J(i,i)));
      last= lambda;
      cout<<"Bound for eigen value: "<<J(i,i)<< " associated with "<< eVec <<" is r"<<i<<"* = "<< lambda<<"."<<endl;
      cout<<"Test on factor give a lower bound of "<< 2*findrstarfactor(firstmorphism, eVec)<<"."<<endl;
      listtest.push_back(testor(lambda, eVec));
      if(J(i,i) ==0 && (i==0 || J(i-1,i)==0)) {
	for(unsigned int j=0; j<nbletters; j++) constrkerofmorphism(nbconstrkerofmorphism,j) = eVec(j);
	boundkerofmorphism(nbconstrkerofmorphism) = lambda;
	nbconstrkerofmorphism++;
      }
      for(unsigned int j=0; j<nbletters; j++) constrVect(nbconstr, j) = eVec(j);
      constr(nbconstr)= lambda;
      nbconstr++;
    }
  }

#ifndef QSQRT3
  /* If the computation are done without exact arithmetics, we take a larger bound. */
  for(vector<testor>::iterator it=listtest.begin();it!=listtest.end();++it) {
    it->bound*=1.00001;
  }
#endif

	
  /////////////////////////////////////////////////////////
  ////////// Computing bounds on the kernel ///////////////
  /////////////////////////////////////////////////////////	
  cout<<endl<< "*************** Computing bounds on the kernel of h ***************"<<endl<<endl;
  /* All the following is also hard-coded.
     The eigenvalues can be checked by any computer algebra system. */
  matrix<F_t,3,3> tempmath = constrVect*baseofKerofh;
  cout << "The matrix corresponding to Q'B in the proof of proposition 10 in the article is: "<<endl<< tempmath<<endl;
  matrix<F_t,3,3> conconadjointh =tempmath*tempmath.transp();
  cout << "(Q'B)x(Q'B)* = "<<endl<< conconadjointh<<endl;
    
#ifndef ADDITIVE
  F_t muminh=F(11,0,10);
  cout << "One can check that the eigenvalues are aproximatively: 165.3 , 29.70 and 1.108 so we can lowerbound the  Umin of the article by "<<muminh<<endl; 
  cout << "The bound for the kernel of h is: "<< sqrt(double(constr.squaredNorm()/muminh))<< "."<<endl<<endl; 
#else
  F_t muminh=F(15,0,100);
  cout << "One can check that the eigenvalues are aproximatively: 29.3 , 3.5 and 0.155 so we can lowerbound the  Umin of the article by "<<muminh<<endl; 
  cout << "The bound for the kernel of h is: "<< sqrt(double(constr.squaredNorm()/muminh))<< "."<<endl<<endl; 
  F_t squaredbound = constr.squaredNorm()/muminh;
#endif
  
  
#ifndef ADDITIVE
  /////////////////////////////////////////////////////////
  ////////////// Computing parents by h ///////////////////
  /////////////////////////////////////////////////////////	
  cout<<endl<< "*************** Computing the patterns parents by h ***************"<<endl<<endl;
	
  //We initialize everithing before the call of getParents
  vector<Pattern> listOfParents; 
  Pattern p;
  listOfParents.push_back(p);
  seenPattern.insert(p);

  vector<cutvec> pOfLeth[nbletters+1];
  init(pOfLeth, h, nbletters);
	
  matrix<int,nbletters, nbletters> Uh,Vh,Diagh;
  smith(Mh, Uh, Diagh, Vh);
  matrix<int,nbletters, nbletters> invVh(Vh.inv()),invUh(Uh.inv());

  assert((Vh*invVh).isId() && (Uh*invUh).isId());
  assert((Uh*Diagh*Vh-Mh).isNull()); // check we have a correct Smith decomposition
 
  vector<vect<int,nbletters> > kerForImproovement(baseofKerofh.getstdvector());
  /* We could reduce the basis of Ker(h) (e.g. by LLL) to improve x_0, and thus the bound and the execution time.
    But this is useless here. */
  
  vector<vect<int,nbletters> > copybaseofKerofh(baseofKerofh.getstdvector());
  // this computes the parents by h that are realizable by the first morphism
  getParents<3>(listOfParents[0], listOfParents, &pOfLeth[1], listtest, invUh, Diagh, invVh, 
 		copybaseofKerofh, muminh, constrVect, constr, kerForImproovement);
  cout<<"There are " <<listOfParents.size()<<" patterns to avoid in the pre-image of h."<<endl<<endl;
  alreadycomputedpreimages.clear();// Reset alreadycomputedpreimages because it is not valid for the first morphism
#else
  cout<<endl<< "*************** Computing the templates corresponding to an additive square by h ***************"<<endl<<endl;
  vector<Pattern> listOfParents; 
  Pattern p;
  vect<int,3> combin;
  int ibound= (int)floor(sqrt((double)(squaredbound)));
   combin(0)=combin(1)=combin(2)=-ibound;
    while(true){
      p.u= baseofKerofh*combin;
      bool ok = true;
      for(unsigned int i=0; i< listtest.size(); i++)
        if(!listtest[i](p.u)){ok = false; break;}
      if(ok) {
        seenPattern.insert(p);
        listOfParents.push_back(p);
        cout <<p<<endl;
      }
      unsigned int l=0;
      while(l < 3 && combin(l)>=ibound) {
        combin(l)=-ibound;
        l++;
      }
      if(l >= 3) break;
      (combin(l))+=1;
    }

#endif  
  /////////////////////////////////////////////////////////
  ////////////////// Computing bounds /////////////////////
  /////////////////////////////////////////////////////////
  
  cout<<endl<< "*************** Computing the bounds on the morphism's kernel ***************"<<endl<<endl;

  /* All the following is also hard-coded. */
  matrix<F_t,2,2> tempmatmorph =(constrkerofmorphism*baseofKerofM);
  cout<<baseofKerofM<<endl;
  cout << constrkerofmorphism<<endl;
  cout<<"The matrix corresponding to QB in the part explaining computation of the parents in the article is: "<<endl<< tempmatmorph <<endl;
  matrix<F_t,2,2> conconadjointmorph =tempmatmorph * tempmatmorph.transp();
	
  cout <<"(Q'B)x(Q'B)* = "<<endl<<conconadjointmorph<<endl;
  F_t muminmorph(1);
  cout<< "One can check that the eigenvalues are 4 and 1 so Umin = "<<muminmorph<<endl; 

  auto projmorphism=constrkerofmorphism;
  auto projboundmorphism=boundkerofmorphism;

  /////////////////////////////////////////////////////////
  ////////////////// Get patterns /////////////////////////
  /////////////////////////////////////////////////////////
  cout<<endl<< "*************** Computing ancestors of the patterns by the first morphism ***************"<<endl<<endl;
	
  //We initialize everithing before calling getParents
  vector<cutvec> pOfLet[nbletters+1];
  init(pOfLet, firstmorphism, nbletters);
	
  matrix<int,nbletters, nbletters> U,V,Diag;
  smith(M, U, Diag, V);
  matrix<int,nbletters, nbletters> invV(V.inv()),invU(U.inv());

  assert((V*invV).isId() && (U*invU).isId());
  assert((U*Diag*V-M).isNull()); // check if we have a correct Smith decomposition
 
  vector<vect<int,nbletters> > kerforimproovementmorph(baseofKerofM.getstdvector());

  vector<vect<int,nbletters> > copybaseofKerofM(baseofKerofM.getstdvector());
  unsigned int ind=0;
  while(ind<listOfParents.size()) {
    getParents<2>(listOfParents[ind], listOfParents, &pOfLet[1], listtest, invU, Diag, invV, 
		  copybaseofKerofM, muminmorph, projmorphism, projboundmorphism, kerforimproovementmorph);
    if(++ind%10000==1)cout<<ind<<" : " <<listOfParents.size()<<endl;
  }

  cout<<"There are "<< listOfParents.size()<<" ancestors."<<endl<<endl;
  alreadycomputedpreimages.clear();

  /////////////////////////////////////////////////////////
  ////////////////////// Get factors //////////////////////
  /////////////////////////////////////////////////////////	
  cout<<endl<< "*************** Finding the factor size to check ***************"<<endl<<endl;
  int Delta=0;
  for(vector<Pattern>::iterator it= listOfParents.begin(), ie= listOfParents.end(); it != ie ; ++it) {
    int normL1=0;
    for(unsigned int i=0; i<it->u.size(); i++) normL1 += abs(it->u(i));
    if(normL1>Delta) Delta = normL1;
  }
  int delta =0;
  for(unsigned int i=0; i<nbletters; i++) if((int)firstmorphism[i].size()> delta) delta = firstmorphism[i].size();
  int lengthToCheck = Delta + 2*delta+2+1;
  cout <<"Using the formula of propositon 8 we can compute Delta = "<< Delta<< " and delta = "<<delta<< 
    " and we know that we have to check factors of length at most "<< lengthToCheck<< "."<<endl;
	
  cout<<endl<< "*************** Finding the small factors ***************"<<endl<<endl;
  vector<vector<string> >factors(lengthToCheck);
  for(unsigned int i=0; i< nbletters; i++) factors[0].push_back(string(1,i+'a'));
	
  set<string>seenFactors;
  computefactorsofsizenfromsizep(firstmorphism, 2, factors[0], factors[1], seenFactors);
  unsigned int templ=0;
  do{
    templ = factors[1].size();
    computefactorsofsizenfromsizep(firstmorphism, 2, factors[1], factors[1], seenFactors);
  }while(templ < factors[1].size());
	
  for(unsigned int i =2; i<factors.size(); i++) {
    computefactorsofsizenfromsizep(firstmorphism, i+1, factors[i-1], factors[i], seenFactors);
    cout << "There are "<< factors[i].size()<<" factors of length "<<i+1 <<"."<<endl;
  }
	
  /////////////////////////////////////////////////////////
  //////////////// Compare factors to patterns ////////////
  /////////////////////////////////////////////////////////	
  cout<<endl<< "*************** Comparing small factors to patterns ***************"<<endl<<endl;
	
  bool abeliansquarefree = true;
  for(vector<Pattern>::iterator it= listOfParents.begin(), ie= listOfParents.end(); it != ie && abeliansquarefree ; ++it) { 
     //if the word realizing the pattern is aubvc
    int diffLength=0; // = |v|-|u|
    for(unsigned int i=0; i<nbletters; i++)	diffLength += it->u(i);
    int lengthLetters = ((it->c[0]!=-1)?1:0)+((it->c[1]!=-1)?1:0)+((it->c[2]!=-1)?1:0); // lengthLetters = |abc| 
    for(unsigned int i=2; i<factors.size() && abeliansquarefree; i++) {
      int l= i+1-diffLength-lengthLetters;
      if(l<0 || l%2!=0) continue;
      l=l/2;// l= |u|
      if( l+diffLength<0) continue; // check|v| >=0
      int p1=l;
      if(it->c[0]!=-1) p1++;
      for(unsigned int j=0; j<factors[i].size() && abeliansquarefree ; j++) {
	// check if the letters of the template are compatible with the factor
	if(it->c[0]!=-1 && it->c[0]!=factors[i][j][0]-'a') continue;
	if(it->c[1]!=-1 && it->c[1]!=factors[i][j][p1]-'a') continue;
	if(it->c[2]!=-1 && it->c[2]!=factors[i][j][i]-'a') continue;
	
	// check if the vector is the same
	unsigned int lc=0;
	if(it->c[0]!=-1) lc++;	
	vect<int,nbletters>xu;
	for(unsigned int k =lc; k<lc+l; k++) xu(factors[i][j][k]-'a')--;
	lc+=l;
				
	if(it->c[1]!=-1) lc++;	
	for(unsigned int k =lc; k<lc+l+diffLength; k++)	xu(factors[i][j][k]-'a')++;
	lc+=l+diffLength;
				
	if(it->c[2]!=-1) lc++;
	assert(lc==i+1);

	if(!(xu == it->u)) continue;

#ifndef ADDITIVE
	cout << "There is an abelian square which appears! "<<endl<<"The template is: "<< *it<<endl;
#else
  cout << "There is an additive square which appears! "<<endl<<"The template is: "<< *it<<endl;
#endif
	cout<<" The factor realizing it is: "<<endl;
	cout<<factors[i][j]<<endl;
	cout<<endl;
	abeliansquarefree = false;
      }
    }
  }
#ifndef ADDITIVE
  if (abeliansquarefree) {
    cout << "The word obtained by applying the second morphism to a fixed-point of the first morphism avoids long abelian-squares."<<endl;

    /* We compute all the small squares */
    int maxL= h[0].size(); // (the morphism is uniform)
    cout<<"By the condition of the Proposition 11 we know that the squares have periode at most "<<maxL<<endl;
    vector<string> possiblesquares;
    set<string> seenFactorsh;
    computefactorsofsizenfromsizep(h, 2*maxL, factors[1], possiblesquares, seenFactorsh);	

    int longestPeriod=0;
    set<string> squares;
    for(unsigned int i=0; i<possiblesquares.size(); i++)
      for(unsigned int s=2; s<=possiblesquares[i].size(); s++)
	longestPeriod = max(longestPeriod, maxSizeOfAbelianSquareEndingAts(possiblesquares[i], s,squares));
    cout<<"The squares appearings are: "<<endl;
    for(set<string>::const_iterator it=squares.begin();it!=squares.end();++it)
      cout<<*it<<" ";
    cout<<endl<<"Number of squares: "<<squares.size()<<"."<<endl;		
    cout<<endl<<"The longest square has period "<<longestPeriod<<"."<<endl;		

  } else
    cout << "The word obtained by applying the second morphism to a fixed-point of the first morphism is NOT abelian-square-free."<<endl;
#else
   if (abeliansquarefree)
    cout << "The word obtained by applying the second morphism to a fixed-point of the first morphism avoids additive-squares."<<endl;
  else
    cout << "The word obtained by applying the second morphism to a fixed-point of the first morphism is NOT additive-square-free."<<endl;
#endif
  return 0;
}