self_avoiding_polygon/other/old/exact/coeffs.cpp

#include "coeffs.hpp"

//*****************************
//* Defintion of externs data *
//****************************

Reel coeffs[num_coeffs];
fmpz_poly_t coeffs_poly[num_coeffs];
fmpz_t coeffs_den;

//******************
//* compute_coeffs *
//******************

void compute_coeffs(){
  CoeffAnalytic ca[num_coeffs];
  
  for(size_t i=0;i<num_coeffs;++i){
    fmpz_poly_init(coeffs_poly[i]);
    fmpq_init(ca[i].a);
    fmpq_init(ca[i].b);
  }
  //C_(0,0)=0
  fmpq_zero(ca[pos(0,0)].a); 
  fmpq_zero(ca[pos(0,0)].b);

  //C_(0,1)=-1
  fmpq_set_si(ca[pos(0,1)].a,-1,1);
  fmpq_zero(ca[pos(0,1)].b);
  
  //Compute C_(i,i)
  fmpq_t r,t,q;
  fmpq_init(r);
  fmpq_init(t);
  fmpq_init(q);
  fmpq_zero(r);
  fmpq_set_si(q,-4,1);
  for(size_t i=1;i<=max_ind_coeffs;++i){
    fmpq_set_si(t,1,2*i-1);
    fmpq_add(r,r,t);
    fmpq_zero(ca[pos(i,i)].a);
    fmpq_mul(t,r,q);
    fmpq_set(ca[pos(i,i)].b,t);
  }
  fmpz_init(coeffs_den);
  fmpz_set(coeffs_den,fmpq_denref(ca[pos(max_ind_coeffs,max_ind_coeffs)].b));
  
  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)
    fmpq* dest=ca[pos(0,j)].a; 
    fmpq_zero(dest);
    fmpq_add(dest,ca[pos(0,j-1)].a,ca[pos(0,j-1)].a);
    fmpq_sub(dest,dest,ca[pos(1,j-1)].a);
    fmpq_add(dest,dest,dest);
    fmpq_sub(dest,dest,ca[pos(0,j-2)].a);
    dest=ca[pos(0,j)].b;
    fmpq_zero(dest);
    fmpq_add(dest,ca[pos(0,j-1)].b,ca[pos(0,j-1)].b);
    fmpq_sub(dest,dest,ca[pos(1,j-1)].b);
    fmpq_add(dest,dest,dest);
    fmpq_sub(dest,dest,ca[pos(0,j-2)].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)
      fmpq_set_si(r,4,1);
      dest=ca[pos(i,j)].a;
      fmpq_zero(dest);
      fmpq_mul(dest,r,ca[pos(i,j-1)].a);
      fmpq_sub(dest,dest,ca[pos(i,j-2)].a);
      fmpq_sub(dest,dest,ca[pos(i-1,j-1)].a);
      fmpq_sub(dest,dest,ca[pos(i+1,j-1)].a);
      dest=ca[pos(i,j)].b;
      fmpq_zero(dest);
      fmpq_mul(dest,r,ca[pos(i,j-1)].b);
      fmpq_sub(dest,dest,ca[pos(i,j-2)].b);
      fmpq_sub(dest,dest,ca[pos(i-1,j-1)].b);
      fmpq_sub(dest,dest,ca[pos(i+1,j-1)].b);
    }
    //C(j-1,j)=2*C(j-1,j-1)-C(j-2,j-1)
    dest=ca[pos(j-1,j)].a;
    fmpq_zero(dest);
    fmpq_add(dest,ca[pos(j-1,j-1)].a,ca[pos(j-1,j-1)].a);
    fmpq_sub(dest,dest,ca[pos(j-2,j-1)].a);
    dest=ca[pos(j-1,j)].b;
    fmpq_zero(dest);
    fmpq_add(dest,ca[pos(j-1,j-1)].b,ca[pos(j-1,j-1)].b);
    fmpq_sub(dest,dest,ca[pos(j-2,j-1)].b);
  }

  fmpz_t temp;
  fmpz_init(temp);
  for(size_t i=0;i<num_coeffs;++i){
 
    fmpz_mul(temp,fmpq_numref(ca[i].a),coeffs_den);
    fmpz_divexact(temp,temp,fmpq_denref(ca[i].a));
    fmpz_poly_set_coeff_fmpz(coeffs_poly[i],0,temp);
    
    fmpz_mul(temp,fmpq_numref(ca[i].b),coeffs_den);
    fmpz_divexact(temp,temp,fmpq_denref(ca[i].b));
    fmpz_poly_set_coeff_fmpz(coeffs_poly[i],1,temp);

  }

}