Sage support

Makefile

@@ -0,0 +1,21 @@
+CPP		= g++
+CPPFLAGS	= -O3 -mavx2 -mfma -g
+EXE		= fp
+LIB		= -lmpfr
+SAGE		= sage
+PYTHON 		= $(SAGE) -python
+MODULE		= fp.so
+
+all: $(EXE) sage
+
+sage: $(MODULE)
+
+
+$(EXE): src/*.hpp src/*.cpp
+	$(CPP) $(CPPFLAGS) $^ -o $@ $(LIB)
+
+$(MODULE): sage/setup.py sage/fp.pyx
+	$(PYTHON) sage/setup.py build_ext --inplace
+
+clean:
+	$(RM) -r build *.so sage/*.cpp $(EXE)

other/ressources/Untitled.ipynb

@@ -0,0 +1,6 @@
+{
+ "cells": [],
+ "metadata": {},
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

other/ressources/Untitled1.ipynb

@@ -0,0 +1,74 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "4614de4b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import fp"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "b72dbcbf",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "ici\n",
+      "get f\n"
+     ]
+    },
+    {
+     "ename": "RuntimeError",
+     "evalue": "Invalid step : \u001b[35mf\u001b[0md",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+      "\u001b[0;31mRuntimeError\u001b[0m                              Traceback (most recent call last)",
+      "Input \u001b[0;32mIn [2]\u001b[0m, in \u001b[0;36m<cell line: 1>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0m \u001b[43mfp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mPolygon\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mfd\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m)\u001b[49m\n",
+      "File \u001b[0;32m~/Seafile/Recherche/self_avoiding_polygon/new/sage/polygon.pxi:13\u001b[0m, in \u001b[0;36mfp.Polygon.__cinit__\u001b[0;34m()\u001b[0m\n\u001b[1;32m     11\u001b[0m       print(\"ici\")\n\u001b[1;32m     12\u001b[0m #      cdef string cpp_str = py_unicode_string.encode(str)\n\u001b[0;32m---> 13\u001b[0m       self.cpp=cpp_Polygon(str.encode('utf-8'))\n\u001b[1;32m     14\u001b[0m \n\u001b[1;32m     15\u001b[0m   def size(self):\n",
+      "\u001b[0;31mRuntimeError\u001b[0m: Invalid step : \u001b[35mf\u001b[0md"
+     ]
+    }
+   ],
+   "source": [
+    "fp.Polygon(\"fd\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "37d2e2fd",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "SageMath 9.5",
+   "language": "sage",
+   "name": "sagemath"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.10.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

other/ressources/fp.py

@@ -0,0 +1,35 @@
+var('tau')
+
+def C_int(dx,dy):
+    return -1/pi*integral(1/tau*(1-((tau-I)/(tau+I))**(dx-dy)*((tau-1)/(tau+1))**(dx+dy)),tau,0,infinity)
+    
+def C_diag(m):
+   return -2/pi*(harmonic_number((2*m-1)/2)+log(4))
+
+var('c1','c2','c3','c4','c5')
+var('x')
+l=[0,c1,c2,c3,c4,c5]
+l=[0,-4*x,-4*x*(1+1/3),-4*x*(1+1/3+1/5),-4*x*(1+1/3+1/5+1/7)]
+
+def c(i,j):
+    if i>j:
+        return c(j,i)
+    if (i,j)==(0,1):
+        return -1
+    if i==j:
+        return l[i]
+    if i==0:
+        return 4*c(0,j-1)-c(0,j-2)-2*c(1,j-1)
+    if i==j-1:
+        return 2*c(j-1,j-1)-2*c(j-2,j)
+    return 4*c(i,j-1)-c(i,j-2)-c(i-1,j-1)-c(i+1,j-1)
+
+var('k')
+
+def d(i,j):
+    if i>j:
+        return d(j,i)
+    r=1
+    for k in range(j):
+        r=r*(2*k+1)
+    return r*c(i,j)

other/ressources/tune.py

@@ -0,0 +1,15 @@
+import fp
+C=fp.Coefficients()
+
+def erreur():
+    L=flatten([[C.get(i,j) for i in range(j+1)] for j in range(1001)])
+    R=flatten(load('65536.txt.sobj'))
+    E=set([abs(L[i]-R[i]) for i in range(len(L))])
+    return max(E)
+
+def surface(imax=1000):
+    P=[]
+    for j in range(imax+1):
+        for i in range(j+1):
+            P.append((i,j,exp(-C.get(i,j))))
+    show(list_plot3d(P,interpolation_type='linear'))

other/trash/float128/.ipynb_checkpoints/Untitled1-checkpoint.ipynb

@@ -0,0 +1,6 @@
+{
+ "cells": [],
+ "metadata": {},
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

other/trash/float128/.ipynb_checkpoints/Untitled2-checkpoint.ipynb

@@ -0,0 +1,6 @@
+{
+ "cells": [],
+ "metadata": {},
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

other/trash/float128/Makefile

@@ -0,0 +1,13 @@
+CPP		= g++
+CPPFLAGS	= -O3 -mavx2 -mfma -g
+EXE		= fp
+LIB		= -lmpfr -lquadmath
+
+all: $(EXE) 
+
+
+$(EXE): src/*.hpp src/*.cpp
+	$(CPP) $(CPPFLAGS) $^ -o $@ $(LIB)
+
+clean:
+	$(RM) -r $(EXE)

other/trash/float128/Makefile~

@@ -0,0 +1,21 @@
+CPP		= g++
+CPPFLAGS	= -O3 -mavx2 -mfma -g
+EXE		= fp
+LIB		= -lmpfr -lquadmath
+SAGE		= sage
+PYTHON 		= $(SAGE) -python
+MODULE		= fp.so
+
+all: $(EXE) #sage
+
+sage: $(MODULE)
+
+
+$(EXE): src/*.hpp src/*.cpp
+	$(CPP) $(CPPFLAGS) $^ -o $@ $(LIB)
+
+$(MODULE): sage/setup.py sage/fp.pyx
+	$(PYTHON) sage/setup.py build_ext --inplace
+
+clean:
+	$(RM) -r build *.so sage/*.cpp $(EXE)

other/trash/float128/src/coefficients.cpp

@@ -0,0 +1,55 @@
+#include "coefficients.hpp"
+
+Coefficients::Coefficients(){
+  mpfr_t* m=new mpfr_t[ncoeffs];
+  for(size_t i=0;i<ncoeffs;++i) mpfr_init2(m[i],prec);
+  mpfr_t a,b,c;
+  mpfr_init2(a,prec);
+  mpfr_init2(b,prec);
+  mpfr_init2(c,prec);
+
+  //C(0,0)=0
+  mpfr_set_zero(m[pos(0,0)],0);
+  //C(1,0)=-1
+  mpfr_set_si(m[pos(1,0)],-1,MPFR_RNDN);
+
+  mpfr_set_si(a,-4,MPFR_RNDN);//a=-4
+  mpfr_const_pi(b,MPFR_RNDN);//b=pi
+  mpfr_div(a,a,b,MPFR_RNDN);//a=-4/pi
+
+  mpfr_set_zero(b,0);//b=0
+  for(size_t i=1;i<=imax;++i){
+    mpfr_set_ui(c,2*i-1,MPFR_RNDN);//c=2i-1
+    mpfr_ui_div(c,1,c,MPFR_RNDN);//c=1/(2i-1)
+    mpfr_add(b,b,c,MPFR_RNDN);// b=sum(1/(2i-1)
+    mpfr_mul(m[pos(i,i)],a,b,MPFR_RNDN);
+  }
+
+  for(size_t j=2;j<=imax;++j){
+    //C(0,j)=4*C(0,j-1)-C(0,j-2)-2*C(1,j-1)
+    mpfr_mul_ui(a,m[pos(0,j-1)],4,MPFR_RNDN);//a=4*C(0,j-1)
+    mpfr_sub(a,a,m[pos(0,j-2)],MPFR_RNDN);//a=4*C(0,j-1)-C(0,j-2)
+    mpfr_mul_ui(b,m[pos(1,j-1)],2,MPFR_RNDN);//b=2*C(1,j-1)
+    mpfr_sub(m[pos(0,j)],a,b,MPFR_RNDN);
+    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)
+      mpfr_mul_ui(a,m[pos(i,j-1)],4,MPFR_RNDN);//a=4*C(i,j-1)
+      mpfr_sub(a,a,m[pos(i,j-2)],MPFR_RNDN);//a=4*C(i,j-1)-C(i,j-2)
+      mpfr_sub(a,a,m[pos(i-1,j-1)],MPFR_RNDN);//a=4*C(i,j-1)-C(i,j-2)-C(i-1,j-1)
+      mpfr_sub(m[pos(i,j)],a,m[pos(i+1,j-1)],MPFR_RNDN);
+    }
+    //C(j-1,j)=2*C(j-1,j-1)-C(j-2,j-1)
+    mpfr_mul_ui(a,m[pos(j-1,j-1)],2,MPFR_RNDN);//a=2*C(j-1,j-1)
+    mpfr_sub(m[pos(j-1,j)],a,m[pos(j-2,j-1)],MPFR_RNDN);
+  }
+
+  for(size_t i=0;i<ncoeffs;++i){
+    coeffs[i]=mpfr_get_d(m[i],MPFR_RNDN);
+    mpfr_clear(m[i]);
+  }
+  mpfr_clear(a);
+  mpfr_clear(b);
+  mpfr_clear(c);
+
+  delete[] m;
+}

other/trash/float128/src/coefficients.hpp

@@ -0,0 +1,30 @@
+#ifndef COEFFICIENTS_HPP
+#define COEFFICIENTS_HPP
+
+#include <mpfr.h>
+
+using Reel = __float128;
+
+class Coefficients{
+private:
+  static constexpr size_t prec=4096;
+  static constexpr size_t imax=1000;
+  static constexpr size_t ncoeffs=((imax+2)*(imax+1))/2;
+  Reel coeffs[ncoeffs];
+  size_t pos(size_t i,size_t j);
+public:
+  Coefficients();
+  Reel get(size_t i,size_t j);
+};
+
+inline size_t
+Coefficients::pos(size_t i,size_t j){
+  return i+j*(j+1)/2;
+}
+
+inline Reel
+Coefficients::get(size_t i,size_t j){
+  return (i<j)?coeffs[pos(i,j)]:coeffs[pos(j,i)];
+}
+
+#endif

other/trash/float128/src/error.cpp

@@ -0,0 +1,20 @@
+#include "error.hpp"
+
+void
+Error::string_error(string msg,string str,size_t pos){
+  Error error;
+  error.msg=msg+" : ";
+  for(size_t i=0;i<min(pos,str.size());++i) error.msg+=str[i];
+  error.msg+="\033[35m";
+  error.msg+=str[pos];
+  error.msg+="\033[0m";
+  for(size_t i=pos+1;i<str.size();++i) error.msg+=str[i];
+  throw(error);
+}
+
+void
+Error::error(string msg){
+  Error error;
+  error.msg=msg;
+  throw(error);
+}

other/trash/float128/src/error.hpp

@@ -0,0 +1,17 @@
+#ifndef ERROR_HPP
+#define ERROR_HPP
+
+#include <string>
+
+using namespace std;
+
+class Error{
+public:
+  string msg;
+  static void string_error(string msg,string str,size_t pos);
+  static void error(string msg);
+
+
+};
+
+#endif

other/trash/float128/src/main.cpp

@@ -0,0 +1,31 @@
+#include "coefficients.hpp"
+#include "polygon.hpp"
+#include <quadmath.h>
+
+int main(){
+  try{
+    Polygon P("u200l200d200r200");
+    printf ("%Qe\n", P.get_fp());
+//    cout<<P.get_fp()<<endl;
+    /*cout<<"Size : "<<P.size()<<endl;
+    cout<<"Graph size : "<<P.graph_size()<<endl;
+    for(size_t i=0;i<P.graph_size();++i){
+      for(size_t j=0;j<P.graph_size();++j){
+        cout<<P.get_coeff_B(i,j)<<' ';
+      }
+      cout<<endl;
+    }
+    cout<<endl;
+    for(size_t i=0;i<P.graph_size();++i){
+      for(size_t j=0;j<P.graph_size();++j){
+        cout<<P.get_coeff_C(i,j)<<' ';
+      }
+      cout<<endl;
+    }*/
+  }
+  catch(const Error& error){
+    cerr<<error.msg<<endl;
+  }
+
+
+}

other/trash/float128/src/matrix.cpp

@@ -0,0 +1,88 @@
+#include "matrix.hpp"
+
+void
+Matrix::init(size_t nrow,size_t ncol){
+  if(data!=nullptr) delete[] data;
+  nr=nrow;
+  nc=ncol;
+  data=new Reel[nr*nc];
+}
+
+//---------------
+// Matrix::clear
+//---------------
+
+void
+Matrix::clear(){
+  for(size_t i=0;i<nr*nc;++i){
+    data[i]=0;
+  }
+}
+
+void
+Matrix::display() const{
+  /*for(size_t i=0;i<nr;++i){
+    for(size_t j=0;j<nc;++j){
+      cout<<get(i,j)<<'\t';
+    }
+    cout<<endl;
+  }*/
+}
+
+void
+Matrix::swap_lines(size_t i,size_t j){
+  for(size_t k=0;k<nc;++k){
+    swap(data[i*nc+k],data[j*nc+k]);
+  }
+}
+
+void
+Matrix::mul_line(size_t i,Reel a){
+  for(size_t k=0;k<nc;++k){
+    data[i*nc+k]*=a;
+  }
+}
+
+void
+Matrix::add_mul_line(size_t i,size_t j,Reel a){
+  for(size_t k=0;k<nc;++k){
+    data[i*nc+k]+=a*data[j*nc+k];
+  }
+}
+
+
+Reel
+Matrix::Gauss(){
+  Reel det=1;
+  size_t np=0; //np=0
+  for(size_t j=0;j<nc;++j){
+    for(size_t p=np;p<nr;++p){
+      Reel c=get(p,j);
+      if(c!=0){
+        det*=c;
+        mul_line(p,1.0/c);
+        for(size_t k=0;k<nr;++k){
+          if(k!=p){
+            add_mul_line(k,p,-get(k,j));
+          }
+        }
+        if(p!=np){
+          swap_lines(np,p);
+          det*=-1;
+        }
+        ++np;
+        break;
+      }
+    }
+  }
+  return det;
+}
+
+Reel
+Matrix::get_diag_square_sym(size_t i) const{
+  Reel res=0;
+  for(size_t k=0;k<nc;++k){
+    res+=get(i,k);
+  }
+  return res;
+}

other/trash/float128/src/matrix.hpp

@@ -0,0 +1,58 @@
+#ifndef MATRIX_HPP
+#define MATRIX_HPP
+
+#include <iostream>
+#include <immintrin.h>
+
+using namespace std;
+
+using Reel = __float128;
+
+class Matrix{
+protected:
+  size_t nr,nc;
+  Reel* data;
+public:
+  Matrix();
+  void init(size_t n);
+  void init(size_t nrow,size_t ncol);
+  Reel get(size_t i,size_t j) const;
+  Reel& get(size_t i,size_t j);
+  void clear();
+  void display() const;
+
+  void swap_lines(size_t i,size_t j);
+  void mul_line(size_t i,Reel a);
+  void add_mul_line(size_t i,size_t j,Reel a);
+  Reel get_diag_square_sym(size_t i) const;
+  Reel Gauss();
+};
+
+//******************
+//* Inline methods *
+//******************
+
+inline
+Matrix::Matrix(){
+  nr=0;
+  nc=0;
+  data=nullptr;
+}
+
+inline void
+Matrix::init(size_t n){
+  return init(n,n);
+}
+
+inline Reel
+Matrix::get(size_t i,size_t j) const{
+  return data[i*nc+j];
+}
+
+inline Reel&
+Matrix::get(size_t i,size_t j){
+  return data[i*nc+j];
+}
+
+
+#endif

other/trash/float128/src/polygon.cpp

@@ -0,0 +1,147 @@
+#include "polygon.hpp"
+
+Coefficients Polygon::coefficients;
+
+Polygon::Polygon(string str){
+  Vertex v(0,0);
+  indices.insert({v,0});
+  vertices.push_back(v);
+  length=0;
+  size_t i=0;
+  size_t l=str.size();
+  if(l==0) return;
+  char c=str[0];
+  while(i<l){
+    char s=c;
+    int64 xs=0;
+    int64 ys=0;
+    switch(s){
+    case 'l':
+      xs=-1;
+      break;
+    case 'r':
+      xs=1;
+      break;
+    case 'u':
+      ys=1;
+      break;
+    case 'd':
+      ys=-1;
+      break;
+    default:
+      Error::string_error("Invalid step",str,i);
+      break;
+    };
+    size_t j=++i;
+    size_t n=0;
+    if(i<l){
+      c=str[i];
+      while(i<l and '0'<=c and c<='9'){
+        n=n*10+c-'0';
+        if(++i==l) break;
+        c=str[i];
+      }
+    }
+    if(i==j) n=1; // No digits after step s
+    for(size_t k=0;k<n;++k){
+      v.x+=xs;
+      v.y+=ys;
+      ++length;
+      auto res=indices.insert({v,length});
+      if(!res.second and (i<l or k<n-1)) Error::string_error("Not self avoiding",str,j-1);
+      if(res.second) vertices.push_back(v);
+    }
+  }
+  if(v.x!=0 or v.y!=0) Error::error("Not a polygon");
+  compute_B();
+  compute_C();
+  compute_M();
+  compute_fp();
+}
+
+void Polygon::add_vertex(const int32& x,const int32& y){
+  Vertex w(x,y);
+  if(indices.find(w)==indices.end()){
+    indices.insert({w,vertices.size()});
+    vertices.push_back(w);
+  }
+}
+
+void
+Polygon::add_neighbours(const int32& x,const int32& y){
+  add_vertex(x,y+1);
+  add_vertex(x,y-1);
+  add_vertex(x-1,y);
+  add_vertex(x+1,y);
+}
+
+void
+Polygon::compute_B(){
+  for(size_t i=0;i<length;++i){
+    int32 x=vertices[i].x;
+    int32 y=vertices[i].y;
+    add_neighbours(x,y);
+  }
+  B.init(vertices.size());
+  B.clear();
+  for(size_t i=0;i<length;++i){
+    Vertex& v=vertices[i];
+    size_t r=indices[Vertex(v.x+1,v.y)];
+    size_t l=indices[Vertex(v.x-1,v.y)];
+    size_t u=indices[Vertex(v.x,v.y+1)];
+    size_t d=indices[Vertex(v.x,v.y-1)];
+    B.get(i,r)=1;
+    B.get(r,i)=1;
+    B.get(i,l)=1;
+    B.get(l,i)=1;
+    B.get(i,u)=1;
+    B.get(u,i)=1;
+    B.get(i,d)=1;
+    B.get(d,i)=1;
+  }
+}
+
+void
+Polygon::compute_C(){
+  C.init(vertices.size());
+  size_t n=vertices.size();
+  for(size_t i=0;i<n;++i){
+    for(size_t j=0;j<n;++j){
+      Vertex& vi=vertices[i];
+      Vertex& vj=vertices[j];
+      double dx=abs(vi.x-vj.x);
+      double dy=abs(vi.y-vj.y);
+      C.get(i,j)=coefficients.get(dx,dy);
+    }
+  }
+}
+
+void
+Polygon::compute_M(){
+  //M=[I+1/4*C*B|1]=[I+1/4*C*tB|1]
+  size_t n=vertices.size();
+  M.init(n,n+1);
+  for(size_t i=0;i<n;++i){
+    for(size_t j=0;j<n;++j){
+      Reel c=0;
+      for(size_t k=0;k<n;++k){
+        c+=C.get(i,k)*B.get(k,j);
+      }
+      M.get(i,j)=c/4;
+    }
+    M.get(i,i)+=1;
+    M.get(i,n)=1;
+  }
+}
+
+void
+Polygon::compute_fp(){
+  size_t n=vertices.size();
+  Reel det=M.Gauss();
+  //cout<<"Det = "<<det<<endl;
+  fp=0;
+  for(size_t i=0;i<n;++i){
+    fp+=B.get_diag_square_sym(i)*M.get(i,n);
+  }
+  fp*=(det*0.25);
+}

other/trash/float128/src/polygon.hpp

@@ -0,0 +1,89 @@
+#ifndef POLYGON_HPP
+#define POLYGON_HPP
+
+
+#include <iostream>
+#include <string>
+#include <vector>
+#include "vertex.hpp"
+#include "matrix.hpp"
+#include "coefficients.hpp"
+#include "error.hpp"
+
+using namespace std;
+
+class Polygon{
+private:
+  static Coefficients coefficients;
+  size_t length;
+  vector<Vertex> vertices;
+  unordered_map<Vertex,uint64> indices;
+  Matrix B;
+  Matrix C;
+  Matrix M;
+  Reel fp;
+  void compute_B();
+  void compute_C();
+  void compute_M();
+  void compute_fp();
+  void add_vertex(const int32& x,const int32& y);
+  void add_neighbours(const int32& x,const int32& y);
+public:
+  static Coefficients& get_coefficients();
+  Polygon();
+  Polygon(string str);
+  size_t size() const;
+  size_t graph_size() const;
+  Vertex vertex(size_t i) const;
+  Reel get_coeff_B(size_t i,size_t j) const;
+  Reel get_coeff_C(size_t i,size_t j) const;
+  Reel get_coeff_M(size_t i,size_t j) const;
+  Reel get_fp() const;
+
+};
+
+inline
+Polygon::Polygon(){
+  length=0;
+}
+
+inline size_t
+Polygon::size() const{
+  return length;
+}
+
+inline size_t
+Polygon::graph_size() const{
+  return vertices.size();
+}
+
+inline Vertex
+Polygon::vertex(size_t i) const{
+  return vertices[i];
+}
+
+inline Coefficients&
+Polygon::get_coefficients(){
+  return coefficients;
+}
+
+inline Reel
+Polygon::get_fp() const{
+  return fp;
+}
+
+inline Reel
+Polygon::get_coeff_B(size_t i,size_t j) const{
+  return B.get(i,j);
+}
+
+inline Reel
+Polygon::get_coeff_C(size_t i,size_t j) const{
+  return C.get(i,j);
+}
+
+inline Reel
+Polygon::get_coeff_M(size_t i,size_t j) const{
+  return M.get(i,j);
+}
+#endif

other/trash/float128/src/vertex.hpp

@@ -0,0 +1,59 @@
+#ifndef VERTEX_HPP
+#define VERTEX_HPP
+
+#include <cstdint>
+#include <unordered_map>
+
+using namespace std;
+
+using int32=int32_t;
+using uint32=uint32_t;
+using uint64=uint64_t;
+using int64=int64_t;
+
+class Vertex{
+public:
+  int32 x,y;
+  Vertex();
+  Vertex(int x,int y);
+  Vertex(const Vertex& v);
+  bool operator==(const Vertex& v);
+ };
+
+inline
+Vertex::Vertex(){
+}
+
+inline
+Vertex::Vertex(int32 _x,int32 _y):x(_x),y(_y){
+}
+
+inline Vertex::Vertex(const Vertex& v){
+  x=v.x;
+  y=v.y;
+}
+
+template<>
+struct std::hash<Vertex>{
+  size_t operator()(const Vertex& v) const;
+};
+
+template<>
+struct std::equal_to<Vertex>{
+  constexpr bool operator()(const Vertex& u,const Vertex& v) const;
+};
+
+
+inline size_t
+std::hash<Vertex>::operator()(const Vertex& v) const{
+  uint64 x=(uint32)v.x;
+  uint64 y=(uint32)v.y;
+  return (x<<32)+y;
+}
+
+inline constexpr bool
+std::equal_to<Vertex>::operator()(const Vertex& u,const Vertex& v) const{
+  return u.x==v.x and u.y==v.y;
+}
+
+#endif

other/trash/mpfr/.ipynb_checkpoints/Untitled1-checkpoint.ipynb

@@ -0,0 +1,6 @@
+{
+ "cells": [],
+ "metadata": {},
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

other/trash/mpfr/.ipynb_checkpoints/Untitled2-checkpoint.ipynb

@@ -0,0 +1,6 @@
+{
+ "cells": [],
+ "metadata": {},
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

other/trash/mpfr/Makefile

@@ -0,0 +1,16 @@
+CPP		= g++
+CPPFLAGS	= -O3 -g
+EXE		= fp
+LIB		= -lmpfr
+OBJS      = obj/adjacency_matrix.o obj/coefficients.o obj/error.o obj/matrix.o obj/polygon.o
+
+all: $(EXE)
+
+obj/%.o: src/%.cpp src/%.hpp
+	$(CPP) $(CPPFLAGS) -c $< -o $@
+
+$(EXE): $(OBJS) src/main.cpp
+	$(CPP) $(CPPFLAGS) $^ $(LIB) -o $@
+
+clean:
+	$(RM) -r obj/*.o $(EXE)

other/trash/mpfr/Makefile~

@@ -0,0 +1,25 @@
+CPP		= g++
+CPPFLAGS	= -O3 -g
+EXE		= fp
+LIB		= -lmpfr
+SAGE		= sage
+PYTHON 		= $(SAGE) -python
+MODULE		= fp.so
+OBJS      = obj/adjacency_matrix.o obj/coefficients.o obj/error.o obj/matrix.o obj/polygon.o
+
+all: $(EXE) #sage
+
+sage: $(MODULE)
+
+
+obj/%.o: src/%.cpp src/%.hpp
+	$(CPP) $(CPPFLAGS) -c $< -o $@
+
+$(EXE): $(OBJS) src/main.cpp
+	$(CPP) $(CPPFLAGS) $^ $(LIB) -o $@
+
+$(MODULE): sage/setup.py sage/fp.pyx
+	$(PYTHON) sage/setup.py build_ext --inplace
+
+clean:
+	$(RM) -r build obj/*.o *.so sage/*.cpp $(EXE)

other/trash/mpfr/src/adjacency_matrix.cpp

@@ -0,0 +1,17 @@
+#include "adjacency_matrix.hpp"
+
+void
+AdjacencyMatrix::init(size_t _n){
+  n=_n;
+  if(data!=nullptr) delete[] data;
+  data=new char[n*n];
+}
+
+unsigned int
+AdjacencyMatrix::get_diag_square_sym(size_t i) const{
+  unsigned int res=0;
+  for(size_t j=0;j<n;++j){
+    res+=get(i,j);
+  }
+  return res;
+}

other/trash/mpfr/src/adjacency_matrix.hpp

@@ -0,0 +1,51 @@
+#ifndef ADJACENCY_MATRIX_HPP
+#define ADJACENCY_MATRIX_HPP
+
+#include <iostream>
+
+using namespace std;
+
+class AdjacencyMatrix{
+protected:
+  size_t n;
+  char* data;
+public:
+  AdjacencyMatrix();
+  ~AdjacencyMatrix();
+  void init(size_t n);
+  void clear();
+  int get(size_t i,size_t j) const;
+  void set(size_t i,size_t j);
+  unsigned int get_diag_square_sym(size_t i) const;
+};
+
+inline
+AdjacencyMatrix::AdjacencyMatrix(){
+  n=0;
+  data=nullptr;
+}
+
+inline
+AdjacencyMatrix::~AdjacencyMatrix(){
+  if(data!=nullptr) delete[] data;
+}
+
+
+inline void
+AdjacencyMatrix::clear(){
+  for(size_t i=0;i<n*n;++i){
+    data[i]=0;
+  }
+}
+
+inline int
+AdjacencyMatrix::get(size_t i,size_t j) const{
+  return data[i*n+j];
+}
+
+inline void
+AdjacencyMatrix::set(size_t i,size_t j){
+  data[i*n+j]=1;
+}
+
+#endif

other/trash/mpfr/src/coefficients.cpp

@@ -0,0 +1,49 @@
+#include "coefficients.hpp"
+
+Coefficients::Coefficients(){
+  coeffs=new mpfr_t[ncoeffs];
+  for(size_t i=0;i<ncoeffs;++i) mpfr_init2(coeffs[i],4096);
+  mpfr_t a,b,c;
+  mpfr_init2(a,prec);
+  mpfr_init2(b,prec);
+  mpfr_init2(c,prec);
+
+  //C(0,0)=0
+  mpfr_set_zero(coeffs[pos(0,0)],0);
+  //C(1,0)=-1
+  mpfr_set_si(coeffs[pos(1,0)],-1,MPFR_RNDN);
+
+  mpfr_set_si(a,-4,MPFR_RNDN);//a=-4
+  mpfr_const_pi(b,MPFR_RNDN);//b=pi
+  mpfr_div(a,a,b,MPFR_RNDN);//a=-4/pi
+
+  mpfr_set_zero(b,0);//b=0
+  for(size_t i=1;i<=imax;++i){
+    mpfr_set_ui(c,2*i-1,MPFR_RNDN);//c=2i-1
+    mpfr_ui_div(c,1,c,MPFR_RNDN);//c=1/(2i-1)
+    mpfr_add(b,b,c,MPFR_RNDN);// b=sum(1/(2i-1)
+    mpfr_mul(coeffs[pos(i,i)],a,b,MPFR_RNDN);
+  }
+
+  for(size_t j=2;j<=imax;++j){
+    //C(0,j)=4*C(0,j-1)-C(0,j-2)-2*C(1,j-1)
+    mpfr_mul_ui(a,coeffs[pos(0,j-1)],4,MPFR_RNDN);//a=4*C(0,j-1)
+    mpfr_sub(a,a,coeffs[pos(0,j-2)],MPFR_RNDN);//a=4*C(0,j-1)-C(0,j-2)
+    mpfr_mul_ui(b,coeffs[pos(1,j-1)],2,MPFR_RNDN);//b=2*C(1,j-1)
+    mpfr_sub(coeffs[pos(0,j)],a,b,MPFR_RNDN);
+    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)
+      mpfr_mul_ui(a,coeffs[pos(i,j-1)],4,MPFR_RNDN);//a=4*C(i,j-1)
+      mpfr_sub(a,a,coeffs[pos(i,j-2)],MPFR_RNDN);//a=4*C(i,j-1)-C(i,j-2)
+      mpfr_sub(a,a,coeffs[pos(i-1,j-1)],MPFR_RNDN);//a=4*C(i,j-1)-C(i,j-2)-C(i-1,j-1)
+      mpfr_sub(coeffs[pos(i,j)],a,coeffs[pos(i+1,j-1)],MPFR_RNDN);
+    }
+    //C(j-1,j)=2*C(j-1,j-1)-C(j-2,j-1)
+    mpfr_mul_ui(a,coeffs[pos(j-1,j-1)],2,MPFR_RNDN);//a=2*C(j-1,j-1)
+    mpfr_sub(coeffs[pos(j-1,j)],a,coeffs[pos(j-2,j-1)],MPFR_RNDN);
+  }
+
+  mpfr_clear(a);
+  mpfr_clear(b);
+  mpfr_clear(c);
+}

other/trash/mpfr/src/coefficients.hpp

@@ -0,0 +1,36 @@
+#ifndef COEFFICIENTS_HPP
+#define COEFFICIENTS_HPP
+
+#include <mpfr.h>
+#include "config.hpp"
+
+
+class Coefficients{
+private:
+
+  static constexpr size_t imax=1000;
+  static constexpr size_t ncoeffs=((imax+2)*(imax+1))/2;
+  mpfr_t* coeffs;
+  size_t pos(size_t i,size_t j);
+public:
+  Coefficients();
+  void set(mpfr_t c,size_t i,size_t j);
+  mpfr_ptr get(size_t i,size_t j);
+
+};
+
+inline size_t
+Coefficients::pos(size_t i,size_t j){
+  return i+j*(j+1)/2;
+}
+
+inline void
+Coefficients::set(mpfr_t c,size_t i,size_t j){
+  mpfr_set(c,(i<j)?coeffs[pos(i,j)]:coeffs[pos(j,i)],MPFR_RNDN);
+}
+
+inline mpfr_ptr
+Coefficients::get(size_t i,size_t j){
+  return(i<j)?coeffs[pos(i,j)]:coeffs[pos(j,i)];
+}
+#endif

other/trash/mpfr/src/config.hpp

@@ -0,0 +1,10 @@
+#ifndef CONFIG_HPP
+#define CONFIG_HPP
+
+#include <iostream>
+
+using namespace std;
+
+constexpr size_t prec=1024;
+
+#endif

other/trash/mpfr/src/error.cpp

@@ -0,0 +1,20 @@
+#include "error.hpp"
+
+void
+Error::string_error(string msg,string str,size_t pos){
+  Error error;
+  error.msg=msg+" : ";
+  for(size_t i=0;i<min(pos,str.size());++i) error.msg+=str[i];
+  error.msg+="\033[35m";
+  error.msg+=str[pos];
+  error.msg+="\033[0m";
+  for(size_t i=pos+1;i<str.size();++i) error.msg+=str[i];
+  throw(error);
+}
+
+void
+Error::error(string msg){
+  Error error;
+  error.msg=msg;
+  throw(error);
+}

other/trash/mpfr/src/error.hpp

@@ -0,0 +1,17 @@
+#ifndef ERROR_HPP
+#define ERROR_HPP
+
+#include <string>
+
+using namespace std;
+
+class Error{
+public:
+  string msg;
+  static void string_error(string msg,string str,size_t pos);
+  static void error(string msg);
+
+
+};
+
+#endif

other/trash/mpfr/src/main.cpp

@@ -0,0 +1,14 @@
+#include "coefficients.hpp"
+#include "polygon.hpp"
+
+int main(){
+  try{
+    Polygon P("u100l100d100r100");
+    mpfr_printf("%Re\n",P.get_fp());
+  }
+  catch(const Error& error){
+    cerr<<error.msg<<endl;
+  }
+
+
+}

other/trash/mpfr/src/matrix.cpp

@@ -0,0 +1,79 @@
+#include "matrix.hpp"
+
+void
+Matrix::init(size_t nrow,size_t ncol){
+  if(data!=nullptr){
+    for(size_t i=0;i<nr*nc;++i) mpfr_clear(data[i]);
+    delete[] data;
+  }
+  nr=nrow;
+  nc=ncol;
+  data=new mpfr_t[nr*nc];
+  for(size_t i=0;i<nr*nc;++i) mpfr_init2(data[i],prec);
+}
+
+Matrix::~Matrix(){
+  if(data!=nullptr){
+    for(size_t i=0;i<nr*nc;++i) mpfr_clear(data[i]);
+    delete[] data;
+  }
+}
+//---------------
+// Matrix::clear
+//---------------
+
+void
+Matrix::clear(){
+  for(size_t i=0;i<nr*nc;++i) mpfr_set_zero(data[i],0);
+}
+
+void
+Matrix::swap_lines(size_t i,size_t j){
+  for(size_t k=0;k<nc;++k){
+    mpfr_swap(data[i*nc+k],data[j*nc+k]);
+  }
+}
+
+void
+Matrix::mul_line(size_t i,mpfr_t x){
+  for(size_t k=0;k<nc;++k){
+    mpfr_mul(data[i*nc+k],data[i*nc+k],x,MPFR_RNDN);
+  }
+}
+
+void
+Matrix::add_mul_line(size_t i,size_t j,mpfr_t x){
+  for(size_t k=0;k<nc;++k){
+    mpfr_fma(data[i*nc+k],data[j*nc+k],x,data[i*nc+k],MPFR_RNDN);
+  }
+}
+
+void
+Matrix::Gauss(mpfr_t det){
+  mpfr_set_ui(det,1,MPFR_RNDN);
+  mpfr_t c;
+  mpfr_init2(c,prec);
+  size_t np=0; //np=0
+  for(size_t j=0;j<nc;++j){
+    for(size_t p=np;p<nr;++p){
+      mpfr_set(c,get(p,j),MPFR_RNDN);
+      if(mpfr_zero_p(c)==0){
+        mpfr_mul(det,det,c,MPFR_RNDN);
+        mpfr_ui_div(c,1,c,MPFR_RNDN);
+        mul_line(p,c);
+        for(size_t k=0;k<nr;++k){
+          if(k!=p){
+            mpfr_neg(c,get(k,j),MPFR_RNDN);
+            add_mul_line(k,p,c);
+          }
+        }
+        if(p!=np){
+          swap_lines(np,p);
+          mpfr_neg(det,det,MPFR_RNDN);
+        }
+        ++np;
+        break;
+      }
+    }
+  }
+}

other/trash/mpfr/src/matrix.hpp

@@ -0,0 +1,64 @@
+#ifndef MATRIX_HPP
+#define MATRIX_HPP
+
+#include <iostream>
+#include <immintrin.h>
+#include <mpfr.h>
+#include "config.hpp"
+
+using namespace std;
+
+class Matrix{
+protected:
+  size_t nr,nc;
+  mpfr_t* data;
+public:
+  Matrix();
+  ~Matrix();
+  void init(size_t n);
+  void init(size_t nrow,size_t ncol);
+  mpfr_srcptr get(size_t i,size_t j) const;
+  mpfr_ptr get(size_t i,size_t j);
+  void set(size_t i,size_t j,mpfr_t x);
+  void clear();
+  void display();
+
+  void swap_lines(size_t i,size_t j);
+  void mul_line(size_t i,mpfr_t a);
+  void add_mul_line(size_t i,size_t j,mpfr_t a);
+  void Gauss(mpfr_t det);
+};
+
+
+//******************
+//* Inline methods *
+//******************
+
+inline
+Matrix::Matrix(){
+  nr=0;
+  nc=0;
+  data=nullptr;
+}
+
+inline void
+Matrix::init(size_t n){
+  return init(n,n);
+}
+
+inline mpfr_srcptr
+Matrix::get(size_t i,size_t j) const{
+  return data[i*nc+j];
+}
+
+inline mpfr_ptr
+Matrix::get(size_t i,size_t j){
+  return data[i*nc+j];
+}
+
+inline void
+Matrix::set(size_t i,size_t j,mpfr_t x){
+  mpfr_set(data[i*nc+j],x,MPFR_RNDN);
+}
+
+#endif

other/trash/mpfr/src/polygon.cpp

@@ -0,0 +1,162 @@
+#include "polygon.hpp"
+
+Coefficients Polygon::coefficients;
+
+Polygon::Polygon(string str){
+  Vertex v(0,0);
+  indices.insert({v,0});
+  vertices.push_back(v);
+  length=0;
+  size_t i=0;
+  size_t l=str.size();
+  if(l==0) return;
+  char c=str[0];
+  while(i<l){
+    char s=c;
+    int64 xs=0;
+    int64 ys=0;
+    switch(s){
+    case 'l':
+      xs=-1;
+      break;
+    case 'r':
+      xs=1;
+      break;
+    case 'u':
+      ys=1;
+      break;
+    case 'd':
+      ys=-1;
+      break;
+    default:
+      Error::string_error("Invalid step",str,i);
+      break;
+    };
+    size_t j=++i;
+    size_t n=0;
+    if(i<l){
+      c=str[i];
+      while(i<l and '0'<=c and c<='9'){
+        n=n*10+c-'0';
+        if(++i==l) break;
+        c=str[i];
+      }
+    }
+    if(i==j) n=1; // No digits after step s
+    for(size_t k=0;k<n;++k){
+      v.x+=xs;
+      v.y+=ys;
+      ++length;
+      auto res=indices.insert({v,length});
+      if(!res.second and (i<l or k<n-1)) Error::string_error("Not self avoiding",str,j-1);
+      if(res.second) vertices.push_back(v);
+    }
+  }
+  if(v.x!=0 or v.y!=0) Error::error("Not a polygon");
+  compute_B();
+  compute_fp();
+}
+
+void Polygon::add_vertex(const int32& x,const int32& y){
+  Vertex w(x,y);
+  if(indices.find(w)==indices.end()){
+    indices.insert({w,vertices.size()});
+    vertices.push_back(w);
+  }
+}
+
+void
+Polygon::add_neighbours(const int32& x,const int32& y){
+  add_vertex(x,y+1);
+  add_vertex(x,y-1);
+  add_vertex(x-1,y);
+  add_vertex(x+1,y);
+}
+
+void
+Polygon::compute_B(){
+  for(size_t i=0;i<length;++i){
+    int32 x=vertices[i].x;
+    int32 y=vertices[i].y;
+    add_neighbours(x,y);
+  }
+  B.init(vertices.size());
+  B.clear();
+  for(size_t i=0;i<length;++i){
+    Vertex& v=vertices[i];
+    size_t r=indices[Vertex(v.x+1,v.y)];
+    size_t l=indices[Vertex(v.x-1,v.y)];
+    size_t u=indices[Vertex(v.x,v.y+1)];
+    size_t d=indices[Vertex(v.x,v.y-1)];
+    B.set(i,r);
+    B.set(r,i);
+    B.set(i,l);
+    B.set(l,i);
+    B.set(i,u);
+    B.set(u,i);
+    B.set(i,d);
+    B.set(d,i);
+  }
+}
+
+void
+Polygon::compute_C(Matrix& C){
+  C.init(vertices.size());
+  size_t n=vertices.size();
+  for(size_t i=0;i<n;++i){
+    for(size_t j=0;j<n;++j){
+      Vertex& vi=vertices[i];
+      Vertex& vj=vertices[j];
+      double dx=abs(vi.x-vj.x);
+      double dy=abs(vi.y-vj.y);
+      C.set(i,j,coefficients.get(dx,dy));
+    }
+  }
+}
+
+void
+Polygon::compute_M(Matrix& M,const Matrix& C){
+  //M=[I+1/4*C*B|1]=[I+1/4*C*tB|1]
+  size_t n=vertices.size();
+  mpfr_t temp;
+  mpfr_init2(temp,prec);
+  M.init(n,n+1);
+  M.clear();
+  for(size_t i=0;i<n;++i){
+    for(size_t j=0;j<n;++j){
+      mpfr_set_zero(temp,0);
+      for(size_t k=0;k<n;++k){
+        if(B.get(k,j)!=0){
+          mpfr_add(temp,temp,C.get(i,k),MPFR_RNDN);
+        }
+      }
+      mpfr_div_ui(temp,temp,4,MPFR_RNDN);
+      M.set(i,j,temp);
+    }
+    mpfr_add_ui(M.get(i,i),M.get(i,i),1,MPFR_RNDN);
+    mpfr_set_ui(M.get(i,n),1,MPFR_RNDN);
+  }
+}
+
+void
+Polygon::compute_fp(){
+  Matrix C,M;
+  compute_C(C);
+  compute_M(M,C);
+  size_t n=vertices.size();
+  mpfr_t det,temp;
+  mpfr_init2(det,prec);
+  mpfr_init2(fp,prec);
+  mpfr_init2(temp,prec);
+  mpfr_set_zero(fp,0);
+  M.Gauss(det);
+  mpfr_printf("Det = %Re\n",det);
+  for(size_t i=0;i<n;++i){
+    mpfr_set_ui(temp,B.get_diag_square_sym(i),MPFR_RNDN);
+    mpfr_fma(fp,temp,M.get(i,n),fp,MPFR_RNDN);
+  }
+  mpfr_mul(fp,fp,det,MPFR_RNDN);
+  mpfr_div_ui(fp,fp,4,MPFR_RNDN);
+  mpfr_clear(det);
+  mpfr_clear(temp);
+}

other/trash/mpfr/src/polygon.hpp

@@ -0,0 +1,77 @@
+#ifndef POLYGON_HPP
+#define POLYGON_HPP
+
+
+#include <iostream>
+#include <string>
+#include <vector>
+#include "vertex.hpp"
+#include "matrix.hpp"
+#include "adjacency_matrix.hpp"
+#include "coefficients.hpp"
+#include "error.hpp"
+
+using namespace std;
+
+class Polygon{
+private:
+  static Coefficients coefficients;
+  size_t length;
+  vector<Vertex> vertices;
+  unordered_map<Vertex,uint64> indices;
+  AdjacencyMatrix B;
+  mpfr_t fp;
+  void compute_B();
+  void compute_C(Matrix& C);
+  void compute_M(Matrix& M,const Matrix& C);
+  void compute_fp();
+  void add_vertex(const int32& x,const int32& y);
+  void add_neighbours(const int32& x,const int32& y);
+public:
+  static Coefficients& get_coefficients();
+  Polygon();
+  Polygon(string str);
+  size_t size() const;
+  size_t graph_size() const;
+  Vertex vertex(size_t i) const;
+  int get_coeff_B(size_t i,size_t j) const;
+  mpfr_srcptr get_fp() const;
+
+};
+
+inline
+Polygon::Polygon(){
+  length=0;
+}
+
+inline size_t
+Polygon::size() const{
+  return length;
+}
+
+inline size_t
+Polygon::graph_size() const{
+  return vertices.size();
+}
+
+inline Vertex
+Polygon::vertex(size_t i) const{
+  return vertices[i];
+}
+
+inline Coefficients&
+Polygon::get_coefficients(){
+  return coefficients;
+}
+
+inline  mpfr_srcptr
+Polygon::get_fp() const{
+  return fp;
+}
+
+inline int
+Polygon::get_coeff_B(size_t i,size_t j) const{
+  return B.get(i,j);
+}
+
+#endif

other/trash/mpfr/src/vertex.hpp

@@ -0,0 +1,59 @@
+#ifndef VERTEX_HPP
+#define VERTEX_HPP
+
+#include <cstdint>
+#include <unordered_map>
+
+using namespace std;
+
+using int32=int32_t;
+using uint32=uint32_t;
+using uint64=uint64_t;
+using int64=int64_t;
+
+class Vertex{
+public:
+  int32 x,y;
+  Vertex();
+  Vertex(int x,int y);
+  Vertex(const Vertex& v);
+  bool operator==(const Vertex& v);
+ };
+
+inline
+Vertex::Vertex(){
+}
+
+inline
+Vertex::Vertex(int32 _x,int32 _y):x(_x),y(_y){
+}
+
+inline Vertex::Vertex(const Vertex& v){
+  x=v.x;
+  y=v.y;
+}
+
+template<>
+struct std::hash<Vertex>{
+  size_t operator()(const Vertex& v) const;
+};
+
+template<>
+struct std::equal_to<Vertex>{
+  constexpr bool operator()(const Vertex& u,const Vertex& v) const;
+};
+
+
+inline size_t
+std::hash<Vertex>::operator()(const Vertex& v) const{
+  uint64 x=(uint32)v.x;
+  uint64 y=(uint32)v.y;
+  return (x<<32)+y;
+}
+
+inline constexpr bool
+std::equal_to<Vertex>::operator()(const Vertex& u,const Vertex& v) const{
+  return u.x==v.x and u.y==v.y;
+}
+
+#endif

sage/coefficients.pxi

@@ -0,0 +1,11 @@
+from cpp_coefficients cimport cpp_Coefficients
+
+cdef class Coefficients:
+  cdef cpp_Coefficients cpp
+
+  def __cinit__(self,init=True):
+    if init:
+      self.cpp=cpp_Coefficients()
+
+  def get(self,i,j):
+    return self.cpp.get(i,j)

sage/cpp_coefficients.pxd

@@ -0,0 +1,7 @@
+cdef extern from "coefficients.cpp":
+  pass
+
+cdef extern from "coefficients.hpp":
+  cdef cppclass cpp_Coefficients "Coefficients":
+    cpp_Coefficients()
+    double get(size_t i,size_t j)

sage/cpp_polygon.pxd

@@ -0,0 +1,37 @@
+from libcpp.string cimport string
+from cpp_vertex cimport cpp_Vertex
+from cpp_coefficients cimport cpp_Coefficients
+cdef extern from "polygon.cpp":
+  pass
+
+cdef extern from "matrix.hpp":
+  pass
+
+cdef extern from "matrix.cpp":
+  pass
+
+cdef extern from "error.hpp":
+  pass
+
+cdef extern from "error.cpp":
+  pass
+
+cdef extern from "py_error.hpp":
+  cdef int raise_py_error()
+
+cdef extern from "polygon.hpp":
+  cdef cppclass cpp_Polygon "Polygon":
+    cpp_Polygon()
+    cpp_Polygon(string) except +raise_py_error
+    size_t size()
+    size_t graph_size()
+    cpp_Vertex vertex(size_t i)
+    double get_coeff_B(size_t i,size_t j)
+    double get_coeff_C(size_t i,size_t j)
+    double get_coeff_M(size_t i,size_t j)
+    double get_fp()
+    @staticmethod
+    cpp_Coefficients get_coefficients()
+
+#cdef extern from "polygon.cpp" namespace "Polygon":
+#  cpp_Coefficients get_coefficients()

sage/cpp_vertex.pxd

@@ -0,0 +1,8 @@
+from libc.stdint cimport int32_t
+
+ctypedef int32_t int32
+
+cdef extern from "vertex.hpp":
+  cdef cppclass cpp_Vertex "Vertex":
+    int32 x
+    int32 y

sage/fp.pyx

@@ -0,0 +1,5 @@
+#distutils: language = c++
+
+include "coefficients.pxi"
+coefficients=Coefficients()
+include "polygon.pxi"

sage/polygon.pxi

@@ -0,0 +1,107 @@
+from cpp_polygon cimport *
+from cpp_vertex cimport cpp_Vertex
+from cpp_coefficients cimport cpp_Coefficients
+from libcpp.string cimport string
+from sage.rings.integer_ring import ZZ
+from sage.rings.real_mpfr import RR
+from sage.matrix.constructor import matrix
+from sage.plot.graphics import Graphics
+from sage.plot.line import line2d
+
+cdef class Polygon:
+  cdef cpp_Polygon cpp
+
+  def __cinit__(self,str=None):
+    if str==None:
+      self.cpp=cpp_Polygon()
+    else:
+      self.cpp=cpp_Polygon(str.encode('utf-8'))
+
+  cpdef size(self):
+    return self.cpp.size()
+
+  cpdef graph_size(self):
+    return self.cpp.graph_size()
+
+  def vertex(self,i):
+    assert(i<self.graph_size())
+    cdef cpp_Vertex v=self.cpp.vertex(i)
+    return (v.x,v.y)
+
+  cpdef vertices(self):
+    res=[]
+    cdef ng=self.cpp.graph_size()
+    cdef cpp_Vertex v
+    for i in range(ng):
+      v=self.cpp.vertex(i)
+      res.append((v.x,v.y))
+    return res
+
+  @staticmethod
+  def get_coefficients():
+    C=Coefficients(init=False)
+    C.cpp=cpp_Polygon.get_coefficients()
+    return C
+
+  cpdef matrix_B(self):
+    cdef ng=self.graph_size()
+    cdef double coeff
+    B=matrix(ZZ,ng)
+    for i in range(ng):
+      for j in range(i,ng):
+        coeff=self.cpp.get_coeff_B(i,j)
+        B[i,j]=int(coeff)
+        B[j,i]=int(coeff)
+    return B
+
+  cpdef matrix_C(self):
+    cdef ng=self.graph_size()
+    cdef double coeff
+    C=matrix(RR,ng)
+    for i in range(ng):
+      for j in range(i,ng):
+        coeff=self.cpp.get_coeff_C(i,j)
+        C[i,j]=coeff
+        C[j,i]=coeff
+    return C
+
+  cpdef matrix_M(self):
+      cdef ng=self.graph_size()
+      cdef double coeff
+      M=matrix(RR,ng,ng+1)
+      for i in range(ng):
+        for j in range(ng+1):
+          M[i,j]=self.cpp.get_coeff_M(i,j)
+      return M
+
+  def fp(self):
+    return self.cpp.get_fp()
+
+  def show(self):
+    B=self.matrix_B()
+    V=self.vertices()
+    X=[v[0] for v in V]
+    Y=[v[1] for v in V]
+    xmin=min(X)
+    xmax=max(X)
+    ymin=min(Y)
+    ymax=max(Y)
+    G=Graphics()
+    for x in range(xmin,xmax+1):
+      G+=line2d([(x,ymin),(x,ymax)],thickness=0.5,rgbcolor=(0.8,0.8,0.8),linestyle="dotted")
+    for y in range(ymin,ymax+1):
+      G+=line2d([(xmin,y),(xmax,y)],thickness=0.5,rgbcolor=(0.8,0.8,0.8),linestyle="dotted")
+    np=self.size()
+    ng=self.graph_size()
+    for i in range(ng):
+      for j in range(i,ng):
+        if B[i,j]==1:
+          if j<np and j==i+1 or (i==0 and j==np-1):
+            G+=line2d([V[i],V[j]],thickness=1.2,rgbcolor=(0.8,0,0))
+          else:
+            G+=line2d([V[i],V[j]],thickness=0.7,rgbcolor=(0.4,0.4,1))
+    G.show(axes=False,aspect_ratio=1)
+
+def Square(l):
+    sl=repr(l)
+    return Polygon('u'+sl+'r'+sl+'d'+sl+'l'+sl)

sage/py_error.hpp

@@ -0,0 +1,13 @@
+#include <Python.h>
+
+#include "error.hpp"
+
+
+void raise_py_error(){
+  try{
+    throw;
+  } catch(Error& error){
+    PyErr_SetString(PyExc_RuntimeError,error.msg.c_str());
+  }
+
+}

sage/setup.py

@@ -0,0 +1,14 @@
+from setuptools import Extension,setup
+from Cython.Build import cythonize
+
+extensions = [
+  Extension("fp",
+    sources=["sage/fp.pyx"],
+    include_dirs=["src"],
+    libraries=[],
+    library_dirs=[],
+    extra_compile_args = ['-O3','-mavx2','-mfma'],
+  )
+]
+
+setup(ext_modules=cythonize(extensions))

sage/walk.py

@@ -0,0 +1,93 @@
+class Walk:
+    def __init__(self,str=None):
+        if str!=None:
+            self.__decompress__(str)
+
+    def __decompress__(self,str):
+        self.xmin=0
+        self.xmax=0
+        self.ymin=0
+        self.ymax=0
+        self.length=0
+        x=0
+        y=0
+        self.steps=[]
+        valid_steps=['l','r','u','b']
+        i=0
+        while i<len(str):
+            s=str[i]
+            if not s in valid_steps:
+                raise AttributeError("Step '"+c+"' unkown")
+            i+=1
+            j=i
+            if i<len(str):
+                j=i
+                v=0
+                while i<len(str) and '0'<=str[i] and str[i]<='9':
+                    v=v*10+ord(str[i])-ord('0')
+                    i+=1
+            if j==i:
+                #No digits
+                v=1
+            self.length+=v
+            for k in range(v):
+                self.steps.append(s)
+            if s=='l':
+                x+=v
+                self.xmax=max(self.xmax,x)
+            elif s=='r':
+                x-=v
+                self.xmin=min(self.xmin,x)
+            elif s=='u':
+                y+=v
+                self.ymax=max(self.ymax,y)
+            else: # s=='b'
+                y-=v
+                self.ymin=min(self.ymin,y)
+
+    def show(self):
+        G=Graphics()
+        for x in range(self.xmin-1,self.xmax+2):
+            G+=line2d([(x,self.ymin-1),(x,self.ymax+1)],thickness=1,rgbcolor=(0.5,0.5,0.5))
+        for y in range(self.ymin-1,self.ymax+2):
+            G+=line2d([(self.xmin-1,y),(self.xmax+1,y)],thickness=1,rgbcolor=(0.5,0.5,0.5))
+        P=[(0,0)]
+        x=0
+        y=0
+        for s in self.steps:
+            if s=='l':
+                x+=1
+            elif s=='r':
+                x-=1
+            elif s=='u':
+                y+=1
+            else:
+                y-=1
+            P.append((x,y))
+        G+=line2d(P,thickness=3)
+        G.show(axes=False,aspect_ratio=1)
+
+    def is_self_avoiding_polygon(self):
+        x=0
+        y=0
+        v=set([(0,0)])
+        l=0
+        for s in self.steps:
+            l+=1
+            if s=='r':
+                x+=1
+            elif s=='l':
+                x-=1
+            elif s=='u':
+                y+=1
+            else:
+                y-=1
+            if l<self.length and (x,y) in v:
+                return False
+            if l==self.length:
+                return (x,y)==(0,0)
+
+    def polygon(self):
+        if not self.is_self_avoiding_polygon():
+            raise RuntimeError("This walk is not a self avoiding polygon")
+        print("Ok lets go")

src/coefficients.cpp

@@ -0,0 +1,55 @@
+#include "coefficients.hpp"
+
+Coefficients::Coefficients(){
+  mpfr_t* m=new mpfr_t[ncoeffs];
+  for(size_t i=0;i<ncoeffs;++i) mpfr_init2(m[i],prec);
+  mpfr_t a,b,c;
+  mpfr_init2(a,prec);
+  mpfr_init2(b,prec);
+  mpfr_init2(c,prec);
+
+  //C(0,0)=0
+  mpfr_set_zero(m[pos(0,0)],0);
+  //C(1,0)=-1
+  mpfr_set_si(m[pos(1,0)],-1,MPFR_RNDN);
+
+  mpfr_set_si(a,-4,MPFR_RNDN);//a=-4
+  mpfr_const_pi(b,MPFR_RNDN);//b=pi
+  mpfr_div(a,a,b,MPFR_RNDN);//a=-4/pi
+
+  mpfr_set_zero(b,0);//b=0
+  for(size_t i=1;i<=imax;++i){
+    mpfr_set_ui(c,2*i-1,MPFR_RNDN);//c=2i-1
+    mpfr_ui_div(c,1,c,MPFR_RNDN);//c=1/(2i-1)
+    mpfr_add(b,b,c,MPFR_RNDN);// b=sum(1/(2i-1)
+    mpfr_mul(m[pos(i,i)],a,b,MPFR_RNDN);
+  }
+
+  for(size_t j=2;j<=imax;++j){
+    //C(0,j)=4*C(0,j-1)-C(0,j-2)-2*C(1,j-1)
+    mpfr_mul_ui(a,m[pos(0,j-1)],4,MPFR_RNDN);//a=4*C(0,j-1)
+    mpfr_sub(a,a,m[pos(0,j-2)],MPFR_RNDN);//a=4*C(0,j-1)-C(0,j-2)
+    mpfr_mul_ui(b,m[pos(1,j-1)],2,MPFR_RNDN);//b=2*C(1,j-1)
+    mpfr_sub(m[pos(0,j)],a,b,MPFR_RNDN);
+    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)
+      mpfr_mul_ui(a,m[pos(i,j-1)],4,MPFR_RNDN);//a=4*C(i,j-1)
+      mpfr_sub(a,a,m[pos(i,j-2)],MPFR_RNDN);//a=4*C(i,j-1)-C(i,j-2)
+      mpfr_sub(a,a,m[pos(i-1,j-1)],MPFR_RNDN);//a=4*C(i,j-1)-C(i,j-2)-C(i-1,j-1)
+      mpfr_sub(m[pos(i,j)],a,m[pos(i+1,j-1)],MPFR_RNDN);
+    }
+    //C(j-1,j)=2*C(j-1,j-1)-C(j-2,j-1)
+    mpfr_mul_ui(a,m[pos(j-1,j-1)],2,MPFR_RNDN);//a=2*C(j-1,j-1)
+    mpfr_sub(m[pos(j-1,j)],a,m[pos(j-2,j-1)],MPFR_RNDN);
+  }
+
+  for(size_t i=0;i<ncoeffs;++i){
+    coeffs[i]=mpfr_get_d(m[i],MPFR_RNDN);
+    mpfr_clear(m[i]);
+  }
+  mpfr_clear(a);
+  mpfr_clear(b);
+  mpfr_clear(c);
+
+  delete[] m;
+}

src/coefficients.hpp

@@ -0,0 +1,28 @@
+#ifndef COEFFICIENTS_HPP
+#define COEFFICIENTS_HPP
+
+#include <mpfr.h>
+
+class Coefficients{
+private:
+  static constexpr size_t prec=4096;
+  static constexpr size_t imax=1000;
+  static constexpr size_t ncoeffs=((imax+2)*(imax+1))/2;
+  double coeffs[ncoeffs];
+  size_t pos(size_t i,size_t j);
+public:
+  Coefficients();
+  double get(size_t i,size_t j);
+};
+
+inline size_t
+Coefficients::pos(size_t i,size_t j){
+  return i+j*(j+1)/2;
+}
+
+inline double
+Coefficients::get(size_t i,size_t j){
+  return (i<j)?coeffs[pos(i,j)]:coeffs[pos(j,i)];
+}
+
+#endif

src/error.cpp

@@ -0,0 +1,20 @@
+#include "error.hpp"
+
+void
+Error::string_error(string msg,string str,size_t pos){
+  Error error;
+  error.msg=msg+" : ";
+  for(size_t i=0;i<min(pos,str.size());++i) error.msg+=str[i];
+  error.msg+="\033[35m";
+  error.msg+=str[pos];
+  error.msg+="\033[0m";
+  for(size_t i=pos+1;i<str.size();++i) error.msg+=str[i];
+  throw(error);
+}
+
+void
+Error::error(string msg){
+  Error error;
+  error.msg=msg;
+  throw(error);
+}

src/error.hpp

@@ -0,0 +1,17 @@
+#ifndef ERROR_HPP
+#define ERROR_HPP
+
+#include <string>
+
+using namespace std;
+
+class Error{
+public:
+  string msg;
+  static void string_error(string msg,string str,size_t pos);
+  static void error(string msg);
+
+
+};
+
+#endif

src/main.cpp

@@ -0,0 +1,29 @@
+#include "coefficients.hpp"
+#include "polygon.hpp"
+
+int main(){
+  try{
+    Polygon P("u200l200d200r200");
+    cout<<P.get_fp()<<endl;
+    /*cout<<"Size : "<<P.size()<<endl;
+    cout<<"Graph size : "<<P.graph_size()<<endl;
+    for(size_t i=0;i<P.graph_size();++i){
+      for(size_t j=0;j<P.graph_size();++j){
+        cout<<P.get_coeff_B(i,j)<<' ';
+      }
+      cout<<endl;
+    }
+    cout<<endl;
+    for(size_t i=0;i<P.graph_size();++i){
+      for(size_t j=0;j<P.graph_size();++j){
+        cout<<P.get_coeff_C(i,j)<<' ';
+      }
+      cout<<endl;
+    }*/
+  }
+  catch(const Error& error){
+    cerr<<error.msg<<endl;
+  }
+
+
+}

src/main.cpp~

@@ -0,0 +1,29 @@
+#include "coefficients.hpp"
+#include "polygon.hpp"
+
+int main(){
+  try{
+    Polygon P("u100l100d100r100");
+    cout<<P.get_fp()<<endl;
+    /*cout<<"Size : "<<P.size()<<endl;
+    cout<<"Graph size : "<<P.graph_size()<<endl;
+    for(size_t i=0;i<P.graph_size();++i){
+      for(size_t j=0;j<P.graph_size();++j){
+        cout<<P.get_coeff_B(i,j)<<' ';
+      }
+      cout<<endl;
+    }
+    cout<<endl;
+    for(size_t i=0;i<P.graph_size();++i){
+      for(size_t j=0;j<P.graph_size();++j){
+        cout<<P.get_coeff_C(i,j)<<' ';
+      }
+      cout<<endl;
+    }*/
+  }
+  catch(const Error& error){
+    cerr<<error.msg<<endl;
+  }
+
+
+}

src/matrix.cpp

@@ -0,0 +1,108 @@
+#include "matrix.hpp"
+
+void
+Matrix::init(size_t nrow,size_t ncol){
+  if(data!=nullptr) delete[] data;
+  nr=nrow;
+  nc=ncol;
+  nc_avx=(nc-1)/4+1;
+  nc_full=4*nc_avx;
+  data=(double*)new __m256d[nr*nc_full];
+}
+
+//---------------
+// Matrix::clear
+//---------------
+
+void
+Matrix::clear(){
+  __m256d* avx=(__m256d*)data;
+  for(size_t i=0;i<nr*nc_avx;++i){
+    avx[i]=zeros;
+  }
+}
+
+void
+Matrix::display() const{
+  for(size_t i=0;i<nr;++i){
+    for(size_t j=0;j<nc;++j){
+      cout<<get(i,j)<<'\t';
+    }
+    cout<<endl;
+  }
+}
+
+void
+Matrix::swap_lines(size_t i,size_t j){
+  __m256d* avx_i=get_avx_row(i);
+  __m256d* avx_j=get_avx_row(j);
+  for(size_t k=0;k<nc_avx;++k){
+    __m256d a=*avx_i;
+    *avx_i=*avx_j;
+    *avx_j=a;
+    ++avx_i;
+    ++avx_j;
+  }
+}
+
+void
+Matrix::mul_line(size_t i,double a){
+  __m256d b=_mm256_set1_pd(a);
+  __m256d* avx=get_avx_row(i);
+  for(size_t k=0;k<nc_avx;++k){
+    *avx=_mm256_mul_pd(*avx,b);
+    ++avx;
+  }
+}
+
+void
+Matrix::add_mul_line(size_t i,size_t j,double a){
+  __m256d b=_mm256_set1_pd(a);
+  __m256d* avx_i=get_avx_row(i);
+  __m256d* avx_j=get_avx_row(j);
+  for(size_t k=0;k<nc_avx;++k){
+    *avx_i=_mm256_fmadd_pd(*avx_j,b,*avx_i);
+    ++avx_i;
+    ++avx_j;
+  }
+}
+
+
+double
+Matrix::Gauss(){
+  double det=1;
+  size_t np=0; //np=0
+  for(size_t j=0;j<nc;++j){
+    for(size_t p=np;p<nr;++p){
+      double c=get(p,j);
+      if(c!=0){
+        det*=c;
+        mul_line(p,1.0/c);
+        for(size_t k=0;k<nr;++k){
+          if(k!=p){
+            add_mul_line(k,p,-get(k,j));
+          }
+        }
+        if(p!=np){
+          swap_lines(np,p);
+          det*=-1;
+        }
+        ++np;
+        break;
+      }
+    }
+  }
+  return det;
+}
+
+double
+Matrix::get_diag_square_sym(size_t i) const{
+  AvxBlock b;
+  b.avx=zeros;
+  const __m256d* avx=get_avx_row(i);
+  for(size_t k=0;k<nc_avx;++k){
+    b.avx=_mm256_fmadd_pd(*avx,*avx,b.avx);
+    ++avx;
+  }
+  return b.data[0]+b.data[1]+b.data[2]+b.data[3];
+}

src/matrix.hpp

@@ -0,0 +1,87 @@
+#ifndef MATRIX_HPP
+#define MATRIX_HPP
+
+#include <iostream>
+#include <immintrin.h>
+
+using namespace std;
+
+class Matrix{
+protected:
+  size_t nr,nc,nc_avx,nc_full;
+  double* data;
+public:
+  Matrix();
+  void init(size_t n);
+  void init(size_t nrow,size_t ncol);
+  double get(size_t i,size_t j) const;
+  double& get(size_t i,size_t j);
+  void clear();
+  void display() const;
+  __m256d* get_avx_row(size_t i);
+  const __m256d* get_avx_row(size_t i) const;
+
+  void swap_lines(size_t i,size_t j);
+  void mul_line(size_t i,double a);
+  void add_mul_line(size_t i,size_t j,double a);
+  double get_diag_square_sym(size_t i) const;
+  double Gauss();
+};
+
+//*****************
+//* 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  AvxBlock{
+  __m256d avx;
+ double data[4];
+};
+
+//******************
+//* Inline methods *
+//******************
+
+inline
+Matrix::Matrix(){
+  nr=0;
+  nc=0;
+  nc_avx=0;
+  nc_full=0;
+  data=nullptr;
+}
+
+inline void
+Matrix::init(size_t n){
+  return init(n,n);
+}
+
+inline double
+Matrix::get(size_t i,size_t j) const{
+  return data[i*nc_full+j];
+}
+
+inline double&
+Matrix::get(size_t i,size_t j){
+  return data[i*nc_full+j];
+}
+
+inline __m256d*
+Matrix::get_avx_row(size_t r){
+  return (__m256d*)(&data[r*nc_full]);
+}
+
+inline const __m256d*
+Matrix::get_avx_row(size_t r) const{
+  return (__m256d*)(&data[r*nc_full]);
+}
+
+#endif

src/polygon.cpp

@@ -0,0 +1,153 @@
+#include "polygon.hpp"
+
+Coefficients Polygon::coefficients;
+
+Polygon::Polygon(string str){
+  Vertex v(0,0);
+  indices.insert({v,0});
+  vertices.push_back(v);
+  length=0;
+  size_t i=0;
+  size_t l=str.size();
+  if(l==0) return;
+  char c=str[0];
+  while(i<l){
+    char s=c;
+    int64 xs=0;
+    int64 ys=0;
+    switch(s){
+    case 'l':
+      xs=-1;
+      break;
+    case 'r':
+      xs=1;
+      break;
+    case 'u':
+      ys=1;
+      break;
+    case 'd':
+      ys=-1;
+      break;
+    default:
+      Error::string_error("Invalid step",str,i);
+      break;
+    };
+    size_t j=++i;
+    size_t n=0;
+    if(i<l){
+      c=str[i];
+      while(i<l and '0'<=c and c<='9'){
+        n=n*10+c-'0';
+        if(++i==l) break;
+        c=str[i];
+      }
+    }
+    if(i==j) n=1; // No digits after step s
+    for(size_t k=0;k<n;++k){
+      v.x+=xs;
+      v.y+=ys;
+      ++length;
+      auto res=indices.insert({v,length});
+      if(!res.second and (i<l or k<n-1)) Error::string_error("Not self avoiding",str,j-1);
+      if(res.second) vertices.push_back(v);
+    }
+  }
+  if(v.x!=0 or v.y!=0) Error::error("Not a polygon");
+  compute_B();
+  compute_C();
+  compute_M();
+  compute_fp();
+}
+
+void Polygon::add_vertex(const int32& x,const int32& y){
+  Vertex w(x,y);
+  if(indices.find(w)==indices.end()){
+    indices.insert({w,vertices.size()});
+    vertices.push_back(w);
+  }
+}
+
+void
+Polygon::add_neighbours(const int32& x,const int32& y){
+  add_vertex(x,y+1);
+  add_vertex(x,y-1);
+  add_vertex(x-1,y);
+  add_vertex(x+1,y);
+}
+
+void
+Polygon::compute_B(){
+  for(size_t i=0;i<length;++i){
+    int32 x=vertices[i].x;
+    int32 y=vertices[i].y;
+    add_neighbours(x,y);
+  }
+  B.init(vertices.size());
+  B.clear();
+  for(size_t i=0;i<length;++i){
+    Vertex& v=vertices[i];
+    size_t r=indices[Vertex(v.x+1,v.y)];
+    size_t l=indices[Vertex(v.x-1,v.y)];
+    size_t u=indices[Vertex(v.x,v.y+1)];
+    size_t d=indices[Vertex(v.x,v.y-1)];
+    B.get(i,r)=1;
+    B.get(r,i)=1;
+    B.get(i,l)=1;
+    B.get(l,i)=1;
+    B.get(i,u)=1;
+    B.get(u,i)=1;
+    B.get(i,d)=1;
+    B.get(d,i)=1;
+  }
+}
+
+void
+Polygon::compute_C(){
+  C.init(vertices.size());
+  size_t n=vertices.size();
+  for(size_t i=0;i<n;++i){
+    for(size_t j=0;j<n;++j){
+      Vertex& vi=vertices[i];
+      Vertex& vj=vertices[j];
+      double dx=abs(vi.x-vj.x);
+      double dy=abs(vi.y-vj.y);
+      C.get(i,j)=coefficients.get(dx,dy);
+    }
+  }
+}
+
+void
+Polygon::compute_M(){
+  //M=[I+1/4*C*B|1]=[I+1/4*C*tB|1]
+  size_t n=vertices.size();
+  M.init(n,n+1);
+  size_t n_avx=(n-1)/4+1;
+  AvxBlock b;
+  for(size_t i=0;i<n;++i){
+    for(size_t j=0;j<n;++j){
+      __m256d* avx_C=C.get_avx_row(i);
+      __m256d* avx_B=B.get_avx_row(j);
+      b.avx=zeros;
+      for(size_t k=0;k<n_avx;++k){
+        b.avx=_mm256_fmadd_pd(*avx_C,*avx_B,b.avx);
+        ++avx_C;
+        ++avx_B;
+      }
+      M.get(i,j)=(0.25)*(b.data[0]+b.data[1]+b.data[2]+b.data[3]);
+    }
+    M.get(i,i)+=1;
+    M.get(i,n)=1;
+  }
+}
+
+void
+Polygon::compute_fp(){
+  size_t n=vertices.size();
+  double det=M.Gauss();
+  cout<<"Det = "<<det<<endl;
+  fp=0;
+  for(size_t i=0;i<n;++i){
+    fp+=B.get_diag_square_sym(i)*M.get(i,n);
+  }
+  fp*=(det*0.25);
+}