Add a new exact version.

Move the old one in approx.

config.hpp

@@ -6,7 +6,7 @@
 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]
+static const size_t max_len=8; // Must be less or equal tin [6,38]
 
 typedef double Reel;
 typedef __int128 Int;

main.cpp

@@ -36,7 +36,7 @@ void treat(Polygon* P){
   delete P;
 }
 int main(){
-  cout<<std::setprecision(20);
+  cout<<std::setprecision(4);
   disp_info();
   compute_coeffs();
   size_t nb[max_len/2];
@@ -52,10 +52,13 @@ int main(){
   while(gen.next()){
     P=new Polygon;
     gen.set(*P);
-    cilk_spawn treat(P);
+    //    if(P->length==4){
+    //cilk_spawn
+    treat(P);
+      //}
     ++total;
   }
-  cilk_sync;
+  //  cilk_sync;
   Reel d=256;
   for(size_t i=1;2*i<max_len;++i){
     Reel l=2*i+2;

polygon.cpp

@@ -133,10 +133,18 @@ Polygon::fp() const{
       AvxC.get(i,j)=get_coeff(dx,dy);
     }
   }
+  cout<<"Matrix B"<<endl;
+  AvxB.display(ne,ne);
+  cout<<"Matrix C"<<endl;
+  AvxC.display(ne,ne);
   AvxMatrix AvxM;
   AvxM.clear();
   AvxM.from_C_B(AvxC,AvxB,ne);
+  cout<<"Matrix M"<<endl;
+  AvxM.display(ne,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);

exact/Makefile

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

exact/coeffs.cpp

@@ -0,0 +1,60 @@
+#include "coeffs.hpp"
+
+//*****************************
+//* Defintion of externs data *
+//****************************
+
+Reel coeffs[num_coeffs];
+fmpz_poly_t coeffs_poly[num_coeffs];
+Int coeffs_den;
+
+//******************
+//* compute_coeffs *
+//******************
+
+void compute_coeffs(){
+  
+  for(size_t i=0;i<num_coeffs;++i){
+    fmpz_poly_init(coeffs_poly[i]);
+  }
+  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};
+  }
+  coeffs_den=ca[pos(max_ind_coeffs,max_ind_coeffs)].b.denominator();
+  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;
+    }
+  }
+  for(size_t i=0;i<num_coeffs;++i){
+    fmpz_poly_set_coeff_si(coeffs_poly[i],0,
+			   (ca[i].a.numerator()*coeffs_den)/ca[i].a.denominator());
+    fmpz_poly_set_coeff_si(coeffs_poly[i],1,
+			   (ca[i].b.numerator()*coeffs_den)/ca[i].b.denominator());
+  }
+}

exact/coeffs.hpp

@@ -0,0 +1,63 @@
+#ifndef COEFFS_HPP
+#define COEFFS_HPP
+
+#include <cassert>
+#include "config.hpp"
+#include "rationnal.hpp"
+#include "flint/fmpz.h"
+#include "flint/fmpz_poly.h"
+
+
+//**********************
+//* 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];
+extern fmpz_poly_t coeffs_poly[num_coeffs];
+extern Int coeffs_den;
+
+//! 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);
+
+//! Return polynom of coefficient C_{i,j}
+void set_poly_coeff(size_t i,size_t j,fmpz_poly_t p);
+
+//! 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)];
+}
+
+inline void
+set_poly_coeff(size_t i,size_t j,fmpz_poly_t p){
+  size_t ind=i<j?pos(i,j):pos(j,i);
+  fmpz_poly_set(p,coeffs_poly[ind]);
+}
+#endif

exact/config.hpp

@@ -0,0 +1,19 @@
+#ifndef CONFIG_HPP
+#define CONFIG_HPP
+
+#include <iostream>
+
+using namespace std;
+
+//! Maximal self avoiding polygon length to consider
+static const size_t max_len=16; // Must be less or equal tin [6,38]
+
+
+typedef double Reel;
+typedef int64_t Int;
+typedef int64_t int64;
+
+
+
+#endif
+

exact/main.cpp

@@ -0,0 +1,99 @@
+#include <iomanip>
+#include <cilk/cilk.h>
+#include <cilk/cilk_api.h>
+#include "config.hpp"
+#include "coeffs.hpp"
+#include "polygon_generator.hpp"
+#include "results.hpp"
+#include "flint/fmpq_poly.h"
+
+ResultsReducer cilk_result;
+
+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;
+  cout<<" Workers number is "<<__cilkrts_get_nworkers()<<endl;
+}
+
+void treat(Polygon* P){
+  size_t l=P->length;
+  size_t i=l/2-1;
+  ++cilk_result.numbers(i);
+  fmpq_poly_t fp;
+  fmpq_poly_init(fp);
+  P->fp(fp);
+  cilk_result.add_fp(i,fp);
+  delete P;
+}
+
+int main(){
+  disp_info();
+  compute_coeffs();
+  size_t nb[max_len/2];
+  PolygonGenerator gen;
+  size_t total=0;
+  Polygon* P;
+  while(gen.next()){
+    P=new Polygon;
+    gen.set(*P);
+    cilk_spawn treat(P);
+    ++total;
+  }
+  cilk_sync;
+
+  //! Mutiplicative coeffcient for the sum of Fp 
+  fmpq_t c;
+  fmpq_init(c);
+  
+  //! Numerator and denominator of the coefficient c
+  fmpz_t d,n;
+  fmpz_init(n);
+  fmpz_init(d);
+
+  //! Init denominator to 256
+  fmpz_set_si(d,256);
+  for(size_t i=1;2*i<max_len;++i){
+
+    //! Length of the SAP
+    Int l=2*i+2;
+    cout<<endl<<"=== Length "<<l<<" ==="<<endl;
+
+    //! Set numerator of c
+    fmpz_set_si(n,2*l);
+    //! Set c
+    fmpq_set_fmpz_frac(c,n,d);
+
+    //! Retrieve Fp sum in res
+    fmpq_poly_t res;
+    fmpq_poly_init(res);
+    fmpq_poly_set(res,cilk_result.get_fp(i));
+
+    //! Multiply res by c
+    fmpq_poly_scalar_mul_fmpq(res,res,c);
+
+    //! Display res
+    char* str=fmpq_poly_get_str_pretty(res,"x");
+    
+    cout<<" > number : "<<cilk_result.numbers(i)<<endl;
+    cout<<" > coeff  : "<<str<<endl; 
+
+    delete str;
+    
+    //! Update denominator
+    fmpz_mul_si(d,d,16);
+  }
+  cout<<endl<<">>> Total of SAP computed : "<<total<<endl;
+  cout<<endl;
+}

exact/polygon.cpp

@@ -0,0 +1,221 @@
+#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
+//-------------
+
+void
+Polygon::fp(fmpq_poly_t res) 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();
+  fmpz_poly_mat_t B,C,M,B2,X,U,V,R;
+  fmpz_poly_t cc,den,det;
+  fmpz_t z;
+  fmpz_poly_mat_init(B,ne,ne);
+  fmpz_poly_mat_init(C,ne,ne);
+  fmpz_poly_mat_init(M,ne,ne);
+  fmpz_poly_mat_init(B2,ne,ne);
+  fmpz_poly_mat_init(X,ne,1);
+  fmpz_poly_mat_init(U,ne,1);
+  fmpz_poly_mat_init(V,1,ne);
+  fmpz_poly_mat_init(R,1,1);
+  fmpz_poly_mat_zero(B);
+  fmpz_poly_init(cc);
+  fmpz_poly_init(det);
+  fmpz_poly_init(den);
+  fmpz_init(z);
+ 
+  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;
+      fmpz_poly_set_si(fmpz_poly_mat_entry(B,i,j),1);
+      fmpz_poly_set_si(fmpz_poly_mat_entry(B,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);
+      set_poly_coeff(dx,dy,fmpz_poly_mat_entry(C,i,j));
+    }
+  }
+  //! Compute B²
+  fmpz_poly_mat_mul(B2,B,B);
+  //! Compute M=C*B
+  fmpz_poly_mat_mul(M,C,B);
+
+  //! M=4*coeffs_denI+CB
+  //! U[i][0]=1;
+  fmpz_poly_set_si(cc,4*coeffs_den);
+  for(size_t i=0;i<ne;++i){
+    fmpz_poly_add(fmpz_poly_mat_entry(M,i,i),fmpz_poly_mat_entry(M,i,i),cc);
+    fmpz_poly_set_si(fmpz_poly_mat_entry(U,i,0),1);
+  }
+  //! Solve system M*X=U, actually M*X=den*U where den is an output
+  fmpz_poly_mat_solve(X,den,M,U);
+
+  //! Store B^2(i,i) in V(0,i)
+  for(size_t i=0;i<ne;++i){
+    fmpz_poly_set(fmpz_poly_mat_entry(V,0,i),fmpz_poly_mat_entry(B2,i,i));
+  }
+
+  //! Compute determinant of M
+  fmpz_poly_mat_det(det,M);
+
+  //! Determine the sign of den/det
+  int s=fmpz_sgn(fmpz_poly_lead(den))*fmpz_sgn(fmpz_poly_lead(det));
+
+  //! Compute R=V*X, R is matrix 1x1
+  fmpz_poly_mat_mul(R,V,X);
+
+  //! Retrieve the unique coefficient of R
+  fmpq_poly_set_fmpz_poly(res,fmpz_poly_mat_entry(R,0,0));
+
+  //! Compute the multiplicative scalar factor
+  fmpz_set_si(z,coeffs_den*4);
+  fmpz_pow_ui(z,z,ne-1);
+  fmpz_mul_si(z,z,4*s);
+
+  //! Multiply res by z and we are done :)
+  fmpq_poly_scalar_div_fmpz(res,res,z);
+
+  //! Clear all flint object
+  fmpz_poly_mat_clear(B);
+  fmpz_poly_mat_clear(C);
+  fmpz_poly_mat_clear(M);
+  fmpz_poly_mat_clear(B2);
+  fmpz_poly_mat_clear(X);
+  fmpz_poly_mat_clear(U);
+  fmpz_poly_mat_clear(V);
+  fmpz_poly_mat_clear(R);
+  fmpz_poly_clear(cc);
+  fmpz_poly_clear(den);
+  fmpz_poly_clear(det);
+  fmpz_clear(z);
+}
+
+//----------------------
+// 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;
+}

exact/polygon.hpp

@@ -0,0 +1,281 @@
+#ifndef POLYGON_HPP
+#define POLYGON_HPP
+
+#include <cstring>
+#include <cassert>
+#include <stack>
+#include <map>
+
+#include "config.hpp"
+#include "coeffs.hpp"
+#include "polygon_step.hpp"
+
+#include "flint/fmpz_poly.h"
+#include "flint/fmpz_poly_q.h"
+#include "flint/fmpz_poly_mat.h"
+#include "flint/fmpq_poly.h"
+
+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
+  void fp(fmpq_poly_t r) 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

exact/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;
+}

exact/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
+

exact/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

exact/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();
+}

exact/rationnal.hpp

@@ -0,0 +1,232 @@
+#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);
+
+  //! Comparaison operator
+  bool operator==(Int n) const;
+
+  bool operator!=(Int n) const;
+
+  bool operator>(Int n) const;
+
+  bool operator<(Int n) const;
+};
+
+//-----------------
+// 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((num<0)?-num:num,den);
+  if(d==0){
+    cout<<num<<" et "<<den<<endl;
+   
+  }
+  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*n*r.den,s*r.num);
+}
+
+inline bool
+Rationnal::operator==(Int n) const{
+  return num==n*den;
+}
+
+inline bool
+Rationnal::operator!=(Int n) const{
+  return num!=n*den;
+}
+
+inline bool
+Rationnal::operator>(Int n) const{
+  return num>n*den;
+}
+
+inline bool
+Rationnal::operator<(Int n) const{
+  return num<n*den;
+}
+
+#endif
+

exact/results.cpp

@@ -0,0 +1,3 @@
+#include "results.hpp"
+
+esultsReducer cilk_results;

exact/results.hpp

@@ -0,0 +1,98 @@
+#ifndef RESULTS_HPP
+#define RESULTS_HPP
+
+#include <cilk/cilk.h>
+#include <cilk/reducer.h>
+#include "config.hpp"
+#include "flint/fmpq_poly.h"
+
+
+//! Number of results
+static const size_t result_size=max_len/2;
+
+//! Structure storing results
+struct Results{
+  //! Number of SAP
+  size_t numbers[result_size];
+  //! Sum of the Fp
+  fmpq_poly_t sum_fp[result_size];
+  //! Empty constructor
+  Results();
+  //! Reset 
+  void reset();
+};
+
+//!**************************
+//! Class for esults reducer
+//!**************************
+class ResultsReducer{
+private:
+  struct Monoid:cilk::monoid_base<Results>{
+    static void reduce(Results* left, Results* Right);
+  };
+  cilk::reducer<Monoid> imp_;
+public:
+  ResultsReducer();
+  //! Obtain reference to numbers
+  size_t& numbers(size_t i);
+  //! Add an Fp coeff
+  //! \param i lenght of the SAP
+  //! \param p poly of the FP
+  void add_fp(size_t i,fmpq_poly_t p);
+  //! Return pointer to th FP poly
+  fmpq_poly_struct* get_fp(size_t i);
+  void reset();
+};
+
+extern ResultsReducer cilk_results;
+
+//***********
+//* Inlines *
+//***********
+
+inline
+Results::Results(){
+  reset();
+}
+
+inline void
+Results::reset(){
+  for(size_t i=0;i<result_size;++i){
+    numbers[i]=0;
+    fmpq_poly_init(sum_fp[i]);
+    fmpq_poly_zero(sum_fp[i]);
+  }
+}
+
+inline
+ResultsReducer::ResultsReducer(){
+  reset();
+}
+inline void
+ResultsReducer::Monoid::reduce(Results* left,Results* right){
+  for(size_t i=0;i<result_size;++i){
+    left->numbers[i]+=right->numbers[i];
+    fmpq_poly_add(left->sum_fp[i],left->sum_fp[i],right->sum_fp[i]);
+  }
+}
+
+inline void
+ResultsReducer::reset(){
+  imp_.view().reset();
+}
+
+inline size_t&
+ResultsReducer::numbers(size_t i){
+  return imp_.view().numbers[i];
+}
+
+inline void
+ResultsReducer::add_fp(size_t i,fmpq_poly_t p){
+  return fmpq_poly_add(imp_.view().sum_fp[i],imp_.view().sum_fp[i],p);
+}
+
+inline fmpq_poly_struct*
+ResultsReducer::get_fp(size_t i){
+  return imp_.view().sum_fp[i];
+}
+#endif