self_avoiding_polygon/other/trash/mpfr/src/polygon.cpp

#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);
}