First commit

Makefile

@@ -0,0 +1,13 @@
+EXE 	= test
+CPP 	= g++
+CFLAGS 	=  -mavx2 -mfma -g -O3 -DNDEBUG
+LIBS	= 
+
+%.o:%.cpp %.hpp config.hpp
+	$(CPP) $(CFLAGS) -c $< -o $@
+
+$(EXE): avx_matrix.o polygon_step.o polygon_generator.o polygon.o rationnal.o coeffs.o main.cpp
+	$(CPP) $(CFLAGS) $^ -o $@ $(LIBS)
+
+clean:
+	-$(RM) *.o  $(EXE) *~

README.md

@@ -0,0 +1 @@
+Exploration for Self Avoiding Polygons

avx_matrix.cpp

@@ -0,0 +1,142 @@
+#include "avx_matrix.hpp"
+
+//------------------
+// AvxMatrix::clear
+//------------------
+
+void
+AvxMatrix::clear(){
+  for(size_t i=0;i<avx_size;++i) avx[i]=zeros;
+}
+
+//--------------------
+// AvxMatrix::display
+//--------------------
+
+void
+AvxMatrix::display(size_t nl,size_t nc) const{
+  for(size_t i=0;i<nl;++i){
+    for(size_t j=0;j<nc;++j){
+      cout<<get(i,j)<<'\t';
+    }
+    cout<<endl;
+  }
+}
+
+//--------------------------------
+// AvxMatrix::get_diag_square_sym
+//--------------------------------
+
+Reel
+AvxMatrix::get_diag_square_sym(size_t i,size_t n) const{
+  size_t n_avx=(n-1)/4+1;
+  Block b;
+  b.avx=zeros;
+  for(size_t k=0;k<n_avx;++k){
+    __m256d a=avx[i*avx_rank+k];
+    b.avx=_mm256_fmadd_pd(a,a,b.avx);
+  }
+  return b.data[0]+b.data[1]+b.data[2]+b.data[3];
+}
+
+//---------------------
+// AvxMatrix::from_C_B
+//---------------------
+
+void
+AvxMatrix::from_C_B(const AvxMatrix& C,const AvxMatrix& B,size_t n){
+  size_t n_avx=(n-1)/4+1;
+  Block c;
+  for(size_t i=0;i<n;++i){
+    for(size_t j=0;j<n;++j){
+      c.avx=zeros;
+      for(size_t k=0;k<n_avx;++k){
+	__m256d a=C.avx[i*avx_rank+k];
+	__m256d b=B.avx[j*avx_rank+k];
+	c.avx=_mm256_fmadd_pd(a,b,c.avx);
+      }
+      get(i,j)=(0.25)*(c.data[0]+c.data[1]+c.data[2]+c.data[3]);
+    }
+    get(i,i)+=1;
+    get(i,n)=1;
+  }
+}
+
+//-----------------------
+// AvxMatrix::swap_lines
+//-----------------------
+
+void
+AvxMatrix::swap_lines(size_t i,size_t j,size_t navx){
+  size_t ind_i=i*avx_rank;
+  size_t ind_j=j*avx_rank;
+  for(size_t k=0;k<navx;++k){
+    __m256d a=avx[ind_i];
+    avx[ind_i]=avx[ind_j];
+    avx[ind_j]=a;
+    ++ind_i;
+    ++ind_j;
+  }
+}
+
+//---------------------
+// AvxMatrix::mul_line
+//---------------------
+
+void
+AvxMatrix::mul_line(size_t i,Reel a,size_t navx){
+  __m256d b=_mm256_set1_pd(a);
+  size_t ind_i=i*avx_rank;
+  for(size_t k=0;k<navx;++k){
+    avx[ind_i]=_mm256_mul_pd(avx[ind_i],b);
+    ++ind_i;
+  }
+}
+
+//-------------------------
+// AvxMatrix::add_mul_line
+//-------------------------
+
+void
+AvxMatrix::add_mul_line(size_t i,size_t j,Reel a,size_t navx){
+  __m256d b=_mm256_set1_pd(a);
+  size_t ind_i=i*avx_rank;
+  size_t ind_j=j*avx_rank;
+  for(size_t k=0;k<navx;++k){
+    avx[ind_i]=_mm256_fmadd_pd(avx[ind_j],b,avx[ind_i]);
+    ++ind_i;
+    ++ind_j;
+  }
+}
+
+//------------------
+// AvxMatrix::Gauss
+//------------------
+
+Reel
+AvxMatrix::Gauss(size_t nl,size_t nc){
+  size_t navx=(nc-1)/4+1;
+  Reel det=1;
+  size_t np=0; //np=0
+  for(size_t j=0;j<nc;++j){
+    for(size_t p=np;p<nl;++p){
+      Reel c=get(p,j);
+      if(c!=0){
+	det*=c;
+	mul_line(p,1.0/c,navx);
+	for(size_t k=0;k<nl;++k){
+	  if(k!=p){
+	    add_mul_line(k,p,-get(k,j),navx);
+	  }
+	}
+	if(p!=np){
+	  swap_lines(np,p,navx);
+	  det*=-1;
+	}
+	++np;
+	break;
+      } 
+    }
+  }
+  return det;
+}

avx_matrix.hpp

@@ -0,0 +1,126 @@
+#ifndef AVX_MATRIX
+#define AVX_MATRIX
+
+#include "config.hpp"
+#include <immintrin.h>
+
+
+//*************
+//* AvxMatrix *
+//*************
+
+//! A Matrix class that take advantage of avx instructions
+class AvxMatrix{
+public:
+  //! Maximum size of the matrix
+  static const size_t rank=3*max_len+1;
+  
+  //! Corresponding length for m256 array
+  static const size_t avx_rank=(rank-1)/4+1;
+
+  //! Real number of colums; must be a multiple of 4
+  static const size_t columns=4*avx_rank;
+
+  //! Total length of the m256 array
+  static const size_t avx_size=avx_rank*columns;
+
+  //! Total length of the double array
+  static const size_t capacity=16*avx_size;
+
+  //! Array of the matrix coefficients
+  union{
+    __m256d avx[avx_size];
+    Reel data[capacity];
+  };
+
+  //! The empty constructor
+  AvxMatrix();
+
+  //! Coefficient accessor
+  //! \param i line position of the coefficient
+  //! \param j column position of the coefficient
+  Reel get(size_t i,size_t j) const;
+  Reel& get(size_t i,size_t j);
+
+  //! Fill matrix with 0
+  void clear();
+  
+  //! Display command
+  //! \param nl number of lines to display
+  //! \param nc number of columns to display
+  void display(size_t nl,size_t nc) const;
+
+  //! Assume the current matrix M is symmetric and
+  //! return the ith diagonal term of M square
+  Reel get_diag_square_sym(size_t i,size_t n) const;
+
+  //! Set the current Matrix to be I+1/4.C.B where B
+  //! is supposed to be symmetric
+  //! \param C a matrix of size nxn
+  //! \param B a symmetric matrix of size nxn
+  //! \param n size of matrices B and C
+  void from_C_B(const AvxMatrix& C,const AvxMatrix& B,size_t n);
+
+  //! Swap lines i and j of the matrix
+  //! \param i line number
+  //! \param j line number (differenet from i)
+  //! \param navx number of m256 blocs to consier per lines
+  void swap_lines(size_t i,size_t j,size_t navx);
+
+  //! Multiply line i by a
+  //! \param i line number
+  //! \param a non zero scalar
+  //! \param navx number of m256 blocs to consider per lines 
+  void mul_line(size_t i,Reel a,size_t navx);
+
+  //! Add to line i a multiple by a of line j
+  //! \param i line number
+  //! \param j line number
+  //! \param a non zero scalar
+  //! \param navx number of m256 blocs to consider per lines
+  void add_mul_line(size_t i,size_t j,Reel a,size_t navx);
+
+  //! Perform Gauss reduction on a top-left submatrix
+  //! Return the corresponding minor
+  //! \param nl number of lines of the submatrix
+  //! \param nc number of columns of the submatrix
+  Reel Gauss(size_t nl,size_t nc);
+};
+
+//*****************
+//* AVX constants *
+//*****************
+
+//! Array of 4 zeros
+static const __m256d zeros=_mm256_set1_pd(0);
+
+//*************
+//* Avx block * 
+//*************
+
+//! Block structure mixing representing m256 as array of 4 doubles
+union  Block{
+  __m256d avx;
+  Reel data[4];
+};
+
+//******************
+//* Inline methods *
+//******************
+
+inline
+AvxMatrix::AvxMatrix(){
+}
+
+inline Reel&
+AvxMatrix::get(size_t i,size_t j){
+  return data[i*columns+j];
+}
+
+inline Reel
+AvxMatrix::get(size_t i,size_t j) const{
+  return data[i*columns+j];
+}
+
+
+#endif

coeffs.cpp

@@ -0,0 +1,47 @@
+#include "coeffs.hpp"
+
+//*****************************
+//* Defintion of externs data *
+//****************************
+
+Reel coeffs[num_coeffs];
+
+
+//******************
+//* compute_coeffs *
+//******************
+
+void compute_coeffs(){
+  CoeffAnalytic ca[num_coeffs];
+  ca[pos(0,0)]={0,0};  //C_(0,0)=0
+  ca[pos(0,1)]={-1,0}; //C_(1,0)=-1
+  //Compute C_(i,i)
+  Rationnal r=0;
+  for(size_t i=1;i<=max_ind_coeffs;++i){
+    Rationnal t(1,2*i-1);
+    r=r+t;
+    ca[pos(i,i)]={0,-4*r};
+  }
+  for(size_t j=2;j<=max_ind_coeffs;++j){
+    //C(0,j)=4*C(0,j-1)-C(0,j-2)-2*C(1,j-1)
+    ca[pos(0,j)].a=4*ca[pos(0,j-1)].a-ca[pos(0,j-2)].a-2*ca[pos(1,j-1)].a;
+    ca[pos(0,j)].b=4*ca[pos(0,j-1)].b-ca[pos(0,j-2)].b-2*ca[pos(1,j-1)].b;
+    for(size_t i=1;i<=j-2;++i){
+      //C(i,j)=4*C(i,j-1)-C(i,j-2)-C(i-1,j-1)-C(i+1,j-1)
+      ca[pos(i,j)].a=4*ca[pos(i,j-1)].a-ca[pos(i,j-2)].a-ca[pos(i-1,j-1)].a-ca[pos(i+1,j-1)].a;
+      ca[pos(i,j)].b=4*ca[pos(i,j-1)].b-ca[pos(i,j-2)].b-ca[pos(i-1,j-1)].b-ca[pos(i+1,j-1)].b;
+    }
+    //C(j-1,j)=2*C(j-1,j-1)-C(j-2,j-1)
+    ca[pos(j-1,j)].a=2*ca[pos(j-1,j-1)].a-ca[pos(j-2,j-1)].a;
+    ca[pos(j-1,j)].b=2*ca[pos(j-1,j-1)].b-ca[pos(j-2,j-1)].b;
+  }
+  Rationnal fpii(1725033,5419351);
+  Reel epii=2.275957200481571e-15; //epii=1/pi-1/fpii
+  for(size_t j=0;j<=max_ind_coeffs;++j){
+    for(size_t i=0;i<=j;++i){
+      size_t ind=pos(i,j);
+      Rationnal f=ca[ind].a+fpii*ca[ind].b;
+      coeffs[pos(i,j)]=(Reel)f+(Reel)ca[ind].b*epii;
+    }
+  }
+}

coeffs.hpp

@@ -0,0 +1,49 @@
+#ifndef COEFFS_HPP
+#define COEFFS_HPP
+
+#include "config.hpp"
+#include "rationnal.hpp"
+
+//**********************
+//* Coefficients C_m,n *
+//**********************
+
+//! Maximal coefficinet index
+static const size_t max_ind_coeffs=max_len/2+2;
+
+//! Number of coefficients
+static const size_t num_coeffs=((max_ind_coeffs+1)*(max_ind_coeffs+2))/2;
+
+//! Array of coefficients
+extern Reel coeffs[num_coeffs];
+
+//! Return the indice of coefficient C_{i,j}
+size_t pos(size_t i,size_t j);
+
+//! Return an approximation of coefficient C_{i,j}
+Reel get_coeff(size_t i,size_t j);
+
+//! Compute all the coefficients
+void compute_coeffs();
+
+//! Represent an analytic coefficient
+struct CoeffAnalytic{
+  //! The coefficient is a-4b/pi
+  Rationnal a,b;
+};
+
+//********************
+//* Inline functions *
+//********************
+
+inline size_t
+pos(size_t i,size_t j){
+  return (j*(j+1))/2+i;
+}
+
+inline Reel
+get_coeff(size_t i,size_t j){
+  return i<j?coeffs[pos(i,j)]:coeffs[pos(j,i)];
+}
+
+#endif

config.hpp

@@ -0,0 +1,35 @@
+#ifndef CONFIG_HPP
+#define CONFIG_HPP
+
+#include <iostream>
+
+using namespace std;
+
+//! Maximal self avoiding polygon length to consider
+static const size_t max_len=20; // Must be less or equal tin [6,38]
+
+typedef double Reel;
+typedef __int128 Int;
+typedef int64_t int64;
+
+ostream& operator<<(ostream& os,const Int& x);
+
+//********************
+//* Inline functions *
+//********************
+
+inline ostream&
+operator<<(ostream& os,const __int128& x){
+  if(x==0) return os<<"0";
+  string str;
+  __int128 ax=(x<0)?-x:x;
+  while(ax!=0){
+    str=to_string((int)ax%10)+str;
+    ax=ax/10;
+  }
+  if(x<0) os<<'-';
+  return os<<str;
+}
+
+#endif
+

main.cpp

@@ -0,0 +1,86 @@
+#include <iomanip>
+#include "config.hpp"
+#include "coeffs.hpp"
+#include "polygon_generator.hpp"
+#include "avx_matrix.hpp"
+
+static const size_t numbers[]=
+  {
+   1,
+   2,
+   7,
+   28,
+   124,
+   588,
+   2938,
+   15268,
+   81826,
+   449572,
+   2521270,
+   14385376,
+   83290424,
+   488384528,
+   2895432660,
+   17332874364,
+   104653427012,
+   636737003384
+  };
+
+void box(string str){
+  size_t w=str.size();
+  cout<<"\u250c";
+  for(size_t i=0;i<w+2;++i) cout<<"\u2500";
+  cout<<"\u2510"<<endl;
+  cout<<"\u2502 "<<str<<" \u2502"<<endl;
+  cout<<"\u2514";
+  for(size_t i=0;i<w+2;++i) cout<<"\u2500";
+  cout<<"\u2518"<<endl;
+}
+
+void disp_info(){
+  box("Compute Fp's series of Self Avoiding Polygon");
+  cout<<endl;
+  cout<<" Maximal length is "<<max_len<<endl;
+}
+
+int main(){
+  cout<<std::setprecision(20);
+  disp_info();
+  compute_coeffs();
+  size_t nb[max_len/2];
+  Reel res[max_len/2];
+  Reel coeff[max_len/2];
+  for(size_t i=0;2*i<max_len;++i){
+    nb[i]=0;
+    res[i]=0;
+  }
+  PolygonGenerator gen;
+  size_t total=0;
+  Polygon P;
+  while(gen.next()){
+    P=gen.get();
+    size_t l=P.length;
+    ++nb[l/2-1];
+    ++total;
+    Reel fp=P.fp();
+    res[l/2-1]+=fp;
+  }
+  Reel d=256;
+  for(size_t i=1;2*i<max_len;++i){
+    Reel l=2*i+2;
+    cout<<endl<<"=== Length "<<l<<" ==="<<endl;
+    Reel r=res[i];
+    Reel n=2*l;
+    coeff[i]=(n*r)/d;
+    cout<<" > number : "<<nb[i]<<endl;
+    cout<<" >  value : "<<coeff[i]<<endl; 
+    d*=16;
+  }
+  cout<<endl<<">>> Total : "<<total<<endl;
+  cout<<endl<<">>> Coefficients list : "<<endl;
+  for(size_t i=1;2*i<max_len;++i){
+    cout<<coeff[i]<<" ";
+  }
+  cout<<endl;
+  
+}

polygon.cpp

@@ -0,0 +1,165 @@
+#include "polygon.hpp"
+
+
+//---------------
+// Polygon::init
+//---------------
+
+void
+Polygon::init(){
+  for(size_t i=0;i<grid_size;++i){
+    grid[i]=empty;
+    grid_histo[i]=0;
+  }
+  for(size_t i=0;i<grid_width;++i){
+    grid[pos(i,0)]=border;
+    grid[pos(i,grid_height-1)]=border;
+  }
+  for(size_t i=0;i<grid_height;++i){
+    grid[pos(0,i)]=border;
+    grid[pos(grid_width-1,i)]=border;
+  }
+  for(size_t i=0;i<xinit;++i){
+    grid[pos(i,1)]=border;
+  }
+  head=pos(xinit,yinit);
+  grid[head]=1;
+}
+
+//--------------------
+// Polygon::operator=
+//--------------------
+
+Polygon&
+Polygon::operator=(const Polygon& P){
+  head=P.head;
+  length=P.length;
+  memcpy(grid,P.grid,grid_size*sizeof(size_t));
+  memcpy(grid_histo,P.grid_histo,grid_size*sizeof(size_t));
+  memcpy(step_histo,P.step_histo,(max_len+1)*sizeof(size_t));
+  return *this;
+}
+
+//------------------
+// Polygon::display
+//------------------
+
+void
+Polygon::display() const{
+    for(size_t z=0;z<grid_height;++z){
+    size_t y=grid_height-z-1;
+    for(size_t x=0;x<grid_width;++x){
+      size_t v=grid[pos(x,y)];
+      if(v==border) cout<<"\033[47m \033[0m";
+      else if(v==empty or step_histo[v]!=grid_histo[pos(x,y)] or v>length) cout<<" ";
+      else cout<<"\033[42m\033[30m"<<" "<<"\033[0m";
+    }
+    cout<<endl;
+  }
+}
+
+//-------------
+// Polygon::fp
+//-------------
+
+Reel
+Polygon::fp() const{
+  map<size_t,InfoVertex> infos;
+  size_t id=0;
+  int64 x=xinit;
+  int64 y=yinit;
+  size_t c=pos(x,y);
+  for(size_t t=1;t<=length;++t){
+    size_t l=left(c);
+    size_t r=right(c);
+    size_t u=up(c);
+    size_t d=down(c);
+    auto it=infos.find(c);
+    if(it==infos.end()){
+      infos[c]=InfoVertex(id++,x,y,l,r,u,d);
+    }
+    else{
+      InfoVertex& info=it->second;
+      info.v[0]=l;
+      info.v[1]=r;
+      info.v[2]=u;
+      info.v[3]=d;
+      info.deg=4;
+    }
+    it=infos.find(l);
+    if(it==infos.end()) infos[l]=InfoVertex(id++,x-1,y);
+    it=infos.find(r);
+    if(it==infos.end()) infos[r]=InfoVertex(id++,x+1,y);
+    it=infos.find(u);
+    if(it==infos.end()) infos[u]=InfoVertex(id++,x,y+1);
+    it=infos.find(d);
+    if(it==infos.end()) infos[d]=InfoVertex(id++,x,y-1);
+    size_t s=step_histo[t+1];
+    if(grid[l]==t+1 and s==grid_histo[l]){
+      c=l;
+      --x;
+    }
+    else if(grid[r]==t+1 and s==grid_histo[r]){
+      c=r;
+      ++x;
+    }
+    else if(grid[u]==t+1 and s==grid_histo[u]){
+      c=u;
+      ++y;
+    }
+    else if(grid[d]==t+1 and s==grid_histo[d]){
+      c=d;
+      --y;
+    }
+  }
+  size_t ne=infos.size();
+  AvxMatrix AvxB;
+  AvxMatrix AvxC;
+  AvxB.clear();
+  AvxC.clear();
+  for(auto it=infos.begin();it!=infos.end();++it){
+    InfoVertex& info=it->second;
+    size_t i=info.id;
+    for(size_t k=0;k<info.deg;++k){
+      size_t j=infos[info.v[k]].id;
+      AvxB.get(i,j)=1;
+      AvxB.get(j,i)=1;
+    }
+    for(auto it2=infos.begin();it2!=infos.end();++it2){
+      InfoVertex& info2=it2->second;
+      size_t j=info2.id;
+      int64 dx=abs(info.x-info2.x);
+      int64 dy=abs(info.y-info2.y);
+      AvxC.get(i,j)=get_coeff(dx,dy);
+    }
+  }
+  AvxMatrix AvxM;
+  AvxM.clear();
+  AvxM.from_C_B(AvxC,AvxB,ne);
+  Reel det=AvxM.Gauss(ne,ne+1);
+  Reel res=0;
+  for(size_t i=0;i<ne;++i){
+    res+=AvxB.get_diag_square_sym(i,ne)*AvxM.get(i,ne);
+  }
+  res*=(det*0.25);
+  return res;
+}
+
+//----------------------
+// Polygon::next_vertex
+//----------------------
+
+size_t
+Polygon::next_vertex(size_t t,size_t c) const{
+  size_t l=left(c);
+  size_t r=right(c);
+  size_t u=up(c);
+  size_t d=down(c);
+  size_t s=step_histo[t+1];
+  if(grid[l]==t+1 and s==grid_histo[l]) return l;
+  if(grid[r]==t+1 and s==grid_histo[r]) return r;
+  if(grid[u]==t+1 and s==grid_histo[u]) return u;
+  if(grid[d]==t+1 and s==grid_histo[d]) return d;
+  assert(false);
+  return 0;
+}

polygon.hpp

@@ -0,0 +1,277 @@
+#ifndef POLYGON_HPP
+#define POLYGON_HPP
+
+#include <cstring>
+#include <cassert>
+#include <stack>
+#include <map>
+
+#include "config.hpp"
+#include "avx_matrix.hpp"
+#include "coeffs.hpp"
+#include "polygon_step.hpp"
+
+using namespace std;
+
+//***********
+//* Polygon *
+//***********
+
+//! Class for self avoiding polygon in the square lattice
+//! Polygon are drown on a grid and hace length bounded
+//! by max_len. The grid have a border prevent the polygon
+//! to go too far. A typical example of the grid is
+//! ##########
+//! #        #
+//! ####S    #
+//! ##########
+//! where S is the starting point and # is border.
+//! Polygon will grow applyng step. A grid case contains
+//! contains the path distance from the starting point.
+//! We store the unique id of a step in an historic grid.
+//! In order to perfoem backward without clearing the grid,
+//! we associate the step id to the corresponding path
+//! distance. Hence the case (i,j) of the grid contribute 
+//! to the polygon if grid[(i,j)]=t and histo_grid[t]=s
+//! and step_histo[s]=t.
+
+class Polygon{
+public:
+  //! Width of the grid
+  static const size_t grid_width=max_len-1;
+
+  //! Height of the grid
+  static const size_t grid_height=max_len/2+2;
+
+  //! Size of the grid
+  static const size_t grid_size=grid_width*grid_height;
+
+  //! Position of the starting point
+  static const int64 xinit=max_len/2-2;
+  static const int64 yinit=1;
+
+  //! Border value
+  static const size_t border=0;
+
+  //! Empty value
+  static const size_t empty=-1;
+
+  //! The living grid the polygon
+  size_t grid[grid_size];
+
+  //! The history grid of the polygon
+  size_t grid_histo[grid_size];
+
+  //! The history step of the polygon
+  size_t step_histo[max_len+1];
+
+  //! Head poistion on the grid
+  size_t head;
+
+  //! Lenght of the polygon
+  size_t length;
+
+  //! Empty constructor
+  Polygon();
+
+  //! Copy operator
+  Polygon& operator=(const Polygon& P);
+  
+  //! Init data of the polygon
+  void init();
+
+  //! Display the polygon
+  void display() const;
+ 
+  //! Apply a step to the polygon
+  //! \param s step to apply
+  void apply(Step& s);
+
+  //! Test is the case at position ind is admissible
+  //! as the next case of the polygon; this at path
+  //! length t.
+  //! \param t path_length of the next case
+  //! \param ind index of the case to test
+  //! \param s last step applied to the polygon
+  bool test(size_t t,size_t ind,Step* s) const;
+
+  //! Test if the polygon is closed
+  bool is_closed() const;
+
+  //! Return indice of the next vertex of the polygon
+  //! \param t path length of the current vertex
+  //! \param c current vertex indice
+  size_t next_vertex(size_t t,size_t c) const;
+  
+  //! Return the fp coefficient of the polygon
+  Reel fp() const;
+
+  //------------------
+  // Static functions
+  //------------------
+
+  //! Return x coordinate of the case with indice ind.
+  //! \param ind indice of the case
+  static int64 x(size_t ind);
+
+  //! Return y coordinate of the case with indice ind.
+  //! \param ind indice of the case
+  static int64 y(size_t ind);
+
+  //! Return indice of the case on the left of case with
+  //! indice ind.
+  //! \param ind indice of the case
+  static size_t left(size_t ind);
+
+  //! Return indice of the case on the right of case with
+  //! indice ind.
+  //! \param ind indice of the case
+  static size_t right(size_t ind);
+  
+  //! Return indice of the case on the up of case with
+  //! indice ind.
+  //! \param ind indice of the case
+  static size_t up(size_t ind);
+  
+  //! Return indice of the case on the down of case with
+  //! indice ind.
+  //! \param ind indice of the case
+  static size_t down(size_t ind);
+  
+  //! Return indice of the case at position (x,y).
+  //! \param x x-coordinate of the case
+  //! \prama y y-coordiante of the case
+  static size_t pos(size_t x,size_t y);
+};
+
+//**************
+//* InfoVertex *
+//**************
+
+//! This class is used to obtained the adjacence matrix
+//! of the haired polygon.
+
+class InfoVertex{
+public:
+  //! Indice of the vertex
+  size_t id;
+
+  //! Coordinates of the vertex
+  int64 x,y;
+  
+  //! Degree of the vertex
+  size_t deg;
+
+  //! Neighbhorhoods of the vertex
+  size_t v[4];
+
+  //! Empty constructor (not explicety used)
+  InfoVertex();
+
+  //! Constructor for a vertex with not yet detected neighbhors.
+  //! \param id indice of the vertex
+  //! \param x x-coordinate of the vertex
+  //! \param y y-coordinate of the vertex
+  InfoVertex(size_t id,int64 x,int64 y);
+
+  //! Constructor for a vertex with four identified neighbhors.
+  //! \param id indice of the vertex
+  //! \param x x-coordinate of the vertex
+  //! \param y y-coordinate of the vertex
+  //! \param l indice of the vertex on the left
+  //! \param r indice of the vertex on the right
+  //! \param u indice of the vertex on the up
+  //! \param d indice of the vertex on the dwon
+  InfoVertex(size_t id,int64 x,int64 y,size_t l,size_t r,size_t u,size_t d);
+};
+
+//********************
+//* Inline functions *
+//********************
+
+//---------
+// Polygon 
+//---------
+
+inline
+Polygon::Polygon(){
+}
+
+inline bool
+Polygon::test(size_t t,size_t ind,Step* s) const{
+  size_t v=grid[ind];
+  if(v>s->t) return true;
+  return grid_histo[ind]!=s->histo[v];
+}
+
+inline void
+Polygon::apply(Step& step){
+  length=step.t;
+  grid[head=step.newh]=step.t;
+  grid_histo[head]=step.num();
+  memcpy(step_histo,step.histo,(step.t+1)*sizeof(size_t));
+}
+
+inline bool
+Polygon::is_closed() const{
+  return grid[down(head)]==1;
+}
+
+//----------------
+// Static Polygon
+//----------------
+
+inline size_t
+Polygon::left(size_t ind){
+  return ind-1;
+}
+
+inline size_t
+Polygon::right(size_t ind){
+  return ind+1;
+}
+
+inline size_t
+Polygon::up(size_t ind){
+  return ind+grid_width;
+}
+
+inline size_t
+Polygon::down(size_t ind){
+  return ind-grid_width;
+}
+
+inline int64
+Polygon::x(size_t ind){
+  return ind%grid_width;
+}
+
+inline int64
+Polygon::y(size_t ind){
+  return ind/grid_width;
+}
+
+inline size_t 
+Polygon::pos(size_t x,size_t y){
+  return  x+y*grid_width;
+}
+
+//------------
+// InfoVertex
+//------------
+
+inline
+InfoVertex::InfoVertex(){
+}
+
+inline
+InfoVertex::InfoVertex(size_t _id,int64 _x,int64 _y):
+  id(_id),x(_x),y(_y),deg(0){
+}
+
+inline
+InfoVertex::InfoVertex(size_t _id,int64 _x,int64 _y,size_t l,size_t r,size_t u,size_t d):
+  id(_id),x(_x),y(_y),deg(4),v{l,r,u,d}{
+}
+  
+#endif

polygon_generator.cpp

@@ -0,0 +1,32 @@
+#include "polygon_generator.hpp"
+
+//------------------------
+// PolygonGenerator::next
+//------------------------
+
+bool
+PolygonGenerator::next(){
+  size_t nd=0;
+  while(not s.empty()){
+    Step step=s.top();
+    s.pop();
+    poly.apply(step);
+    if(poly.is_closed()) return true;
+    size_t h=poly.head;
+    int64 xh=h%Polygon::grid_width;
+    int64 yh=h/Polygon::grid_width;
+    ++nd;
+    if(abs(xh-Polygon::xinit)+abs(yh-Polygon::yinit)+step.t<=max_len+1){
+      size_t l=Polygon::left(h);
+      size_t r=Polygon::right(h);
+      size_t d=Polygon::down(h);
+      size_t u=Polygon::up(h);
+      size_t nt=step.t+1;
+      if(poly.test(nt,l,&step)) s.emplace(&step,l,++n);
+      if(poly.test(nt,r,&step)) s.emplace(&step,r,++n);
+      if(poly.test(nt,d,&step)) s.emplace(&step,d,++n);
+      if(poly.test(nt,u,&step)) s.emplace(&step,u,++n);
+    }
+  }
+  return false;
+}

polygon_generator.hpp

@@ -0,0 +1,61 @@
+#ifndef POLYGON_GENERATOR_HPP
+#define POLYGON_GENRATOR_HPP
+
+#include "polygon.hpp"
+
+//********************
+//* PolygonGenerator *
+//********************
+
+//! Generates all self avoidinfg polygon upto length max_len.
+//! The class contains a unique instance of a polygon. We use
+//! a stack to store each Step of the exploration.
+
+class PolygonGenerator{
+private:
+  //! A stack of Step
+  stack<Step> s;
+
+  //! The unique polygon
+  Polygon poly;
+
+public:
+  //! The id of the next step
+  size_t n;
+
+  //! Empty constructor
+  PolygonGenerator();
+
+  //! Explore until the next self avoiding polygon
+  //! Return is such a polygon exists and false otherwise.
+  bool next();
+
+  //! Return a reference to the current self avoiding polygon.
+  const Polygon& get() const;
+
+  //! Set the current self avoiding polygon to P
+  void set(Polygon& P) const;
+};
+
+//********************
+//* Inline functions *
+//********************
+
+inline
+PolygonGenerator::PolygonGenerator(){
+  n=1;
+  poly.init();
+  s.emplace(2,poly.head,Polygon::right(poly.head),n);
+}
+
+inline const Polygon&
+PolygonGenerator::get() const{
+  return poly;
+}
+
+inline void
+PolygonGenerator::set(Polygon& P) const{
+  P=poly;
+}
+#endif
+

polygon_step.hpp

@@ -0,0 +1,39 @@
+#ifndef POLYGON_STEP_HPP
+#define POLYGON_STEP_HPP
+
+#include "config.hpp"
+
+
+struct Step{
+  size_t histo[max_len+1];
+  size_t t,oldh,newh;
+  Step(size_t t,size_t oldh,size_t newh,size_t n);
+  Step(Step* s,size_t newh,size_t n);
+  size_t num() const;
+};
+
+
+inline
+Step::Step(Step* s,size_t _newh,size_t n){
+  memcpy(histo,s->histo,(1+s->t)*sizeof(size_t));
+  t=s->t+1;
+  histo[t]=n;
+  oldh=s->newh;
+  //cout<<"Step histo : "<<histo[0]<<endl;
+  newh=_newh;
+}
+
+inline
+Step::Step(size_t _t,size_t _oldh,size_t _newh,size_t n){
+  t=_t;oldh=_oldh;newh=_newh;
+  histo[0]=0;
+  histo[1]=0;
+  histo[_t]=1;
+}
+
+inline size_t
+Step::num() const{
+  return histo[t];
+}
+
+#endif

rationnal.cpp

@@ -0,0 +1,24 @@
+#include "rationnal.hpp"
+
+//-----
+// gcd
+//-----
+Int gcd(Int a,Int b){
+  Int r=a;
+  Int rp=b;
+  while(r!=0){
+    Int q=rp/r;
+    Int t=r;
+    r=rp-(q*r);
+    rp=t;
+  }
+  return rp;
+}
+
+//------------
+// operator<<
+//------------
+ostream& operator<<(ostream& os,const Rationnal& r){
+  if(r.denominator()==1) return os<<r.numerator();
+  return os<<r.numerator()<<'/'<<r.denominator();
+}

rationnal.hpp

@@ -0,0 +1,199 @@
+#ifndef RATIONNAL_HPP
+#define RATIONNAL_HPP
+
+#include <iostream>
+#include <cstdint>
+#include <iostream>
+#include "config.hpp"
+
+using namespace std;
+
+//*************
+//* Rationnal *
+//*************
+
+//! Class for representing rationnals.
+class Rationnal{
+private:
+  //! Numerator
+  Int num;
+  //! Denominator
+  Int den;
+  //! Normalise th fraction
+  void normalize();
+public:
+  //! Construct the null rationnal
+  Rationnal();
+  
+  //! Construct the rationnal n
+  //! \param n integer
+  Rationnal(Int n);
+
+  //! Construct the rationnal n/d
+  //! \param n numerator integer
+  //! \param d denominator integer
+  Rationnal(Int n,Int d);
+
+  //! Copy constructor
+  //! \param r rationnal to copy
+  Rationnal(const Rationnal& r);
+
+  //! Copy operator
+  //! \param r rationnal to copy
+  Rationnal& operator=(const Rationnal& r);
+
+  //! Numerator accessor
+  Int numerator() const;
+
+  //! Denominator accessor
+  Int denominator() const;
+
+  //! Return an approximation of the rationnal
+  explicit operator Reel() const;
+
+  //! Addition operators 
+  Rationnal operator+(const Rationnal& r) const;
+  Rationnal operator+(Int n) const;
+  friend Rationnal operator+(Int n,const Rationnal& r);
+
+  //! Substration operators
+  Rationnal operator-(const Rationnal& r) const;
+  Rationnal operator-(Int n) const;
+  friend Rationnal operator-(Int n,const Rationnal& r);
+
+  //! Multiplication operators
+  Rationnal operator*(const Rationnal& r) const;
+  Rationnal operator*(Int n) const;
+  friend Rationnal operator*(Int n,const Rationnal& r);
+
+  //! Division operators
+  Rationnal operator/(const Rationnal& r) const;
+  Rationnal operator/(Int n) const;
+  friend Rationnal operator/(Int n,const Rationnal& r);
+};
+
+//-----------------
+// Other functions
+//-----------------
+ostream& operator<<(ostream& os,const Rationnal& r);
+
+
+Int gcd(Int a,Int b);
+
+//********************
+//* Inline functions *
+//********************
+
+inline
+Rationnal::Rationnal():num(0),den(1){
+}
+
+inline
+Rationnal::Rationnal(Int n):num(n),den(1){
+}
+
+inline
+Rationnal::Rationnal(Int n,Int d):num(n),den(d){
+  normalize();
+}
+
+inline
+Rationnal::Rationnal(const Rationnal& r):num(r.num),den(r.den){
+}
+
+inline Rationnal&
+Rationnal::operator=(const Rationnal& r){
+  num=r.num;
+  den=r.den;
+  return *this;
+}
+  
+inline Int
+Rationnal::numerator() const{
+  return num;
+}
+
+inline Int
+Rationnal::denominator() const{
+  return den;
+}
+
+inline 
+Rationnal::operator Reel() const{
+  Reel n=num;
+  Reel d=den;
+  return n/d;
+}
+
+inline void
+Rationnal::normalize(){
+  Int d=gcd(abs(num),den);
+  num/=d;
+  den/=d;
+}
+
+inline Rationnal
+Rationnal::operator+(const Rationnal& r) const{
+  return Rationnal(num*r.den+r.num*den,den*r.den);
+}
+
+inline Rationnal
+Rationnal::operator+(Int n) const{
+  return Rationnal(num+n*den,den);
+}
+
+inline Rationnal
+Rationnal::operator-(const Rationnal& r) const{
+  return Rationnal(num*r.den-r.num*den,den*r.den);
+}
+
+inline Rationnal
+Rationnal::operator-(Int n) const{
+  return Rationnal(num-n*den,den);
+}
+
+inline Rationnal
+Rationnal::operator*(const Rationnal& r) const{
+  return Rationnal(num*r.num,den*r.den);
+}
+
+inline Rationnal
+Rationnal::operator*(Int n) const{
+  return Rationnal(num*n,den);
+}
+
+inline Rationnal
+Rationnal::operator/(const Rationnal& r) const{
+  Int s=(r.num<0)?-1:1;
+  return Rationnal(s*num*r.den,s*den*r.num);
+}
+
+inline Rationnal
+Rationnal::operator/(Int n) const{
+  Int s=(n<0)?-1:1;
+  return Rationnal(s*num,s*den*n);
+}
+
+inline Rationnal
+operator+(Int n,const Rationnal& r){
+  return Rationnal(n*r.den+r.num,r.den);
+}
+
+inline Rationnal
+operator-(Int n,const Rationnal& r){
+  return Rationnal(-n*r.den+r.num,r.den);
+}
+
+inline Rationnal
+operator*(Int n,const Rationnal& r){
+  return Rationnal(n*r.num,r.den);
+}
+
+inline Rationnal
+operator/(Int n,const Rationnal& r){
+  Int s=(n<0)?-1:1;
+  return Rationnal(s*r.num,s*r.den*n);
+}
+
+#endif
+