self_avoiding_polygon/coeffs.cpp

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