self_avoiding_polygon/other/old/exact/polygon.cpp

#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_mul_si(z,coeffs_den,4);
  fmpz_poly_set_fmpz(cc,z);
  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_mul_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;
}