self_avoiding_polygon/exact/coeffs.cpp

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