#include "coefficients.hpp"

Coefficients::Coefficients(){
  mpfr_t* m=new mpfr_t[ncoeffs];
  for(size_t i=0;i<ncoeffs;++i) mpfr_init2(m[i],prec);
  mpfr_t a,b,c;
  mpfr_init2(a,prec);
  mpfr_init2(b,prec);
  mpfr_init2(c,prec);

  //C(0,0)=0
  mpfr_set_zero(m[pos(0,0)],0);
  //C(1,0)=-1
  mpfr_set_si(m[pos(1,0)],-1,MPFR_RNDN);

  mpfr_set_si(a,-4,MPFR_RNDN);//a=-4
  mpfr_const_pi(b,MPFR_RNDN);//b=pi
  mpfr_div(a,a,b,MPFR_RNDN);//a=-4/pi

  mpfr_set_zero(b,0);//b=0
  for(size_t i=1;i<=imax;++i){
    mpfr_set_ui(c,2*i-1,MPFR_RNDN);//c=2i-1
    mpfr_ui_div(c,1,c,MPFR_RNDN);//c=1/(2i-1)
    mpfr_add(b,b,c,MPFR_RNDN);// b=sum(1/(2i-1)
    mpfr_mul(m[pos(i,i)],a,b,MPFR_RNDN);
  }

  for(size_t j=2;j<=imax;++j){
    //C(0,j)=4*C(0,j-1)-C(0,j-2)-2*C(1,j-1)
    mpfr_mul_ui(a,m[pos(0,j-1)],4,MPFR_RNDN);//a=4*C(0,j-1)
    mpfr_sub(a,a,m[pos(0,j-2)],MPFR_RNDN);//a=4*C(0,j-1)-C(0,j-2)
    mpfr_mul_ui(b,m[pos(1,j-1)],2,MPFR_RNDN);//b=2*C(1,j-1)
    mpfr_sub(m[pos(0,j)],a,b,MPFR_RNDN);
    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)
      mpfr_mul_ui(a,m[pos(i,j-1)],4,MPFR_RNDN);//a=4*C(i,j-1)
      mpfr_sub(a,a,m[pos(i,j-2)],MPFR_RNDN);//a=4*C(i,j-1)-C(i,j-2)
      mpfr_sub(a,a,m[pos(i-1,j-1)],MPFR_RNDN);//a=4*C(i,j-1)-C(i,j-2)-C(i-1,j-1)
      mpfr_sub(m[pos(i,j)],a,m[pos(i+1,j-1)],MPFR_RNDN);
    }
    //C(j-1,j)=2*C(j-1,j-1)-C(j-2,j-1)
    mpfr_mul_ui(a,m[pos(j-1,j-1)],2,MPFR_RNDN);//a=2*C(j-1,j-1)
    mpfr_sub(m[pos(j-1,j)],a,m[pos(j-2,j-1)],MPFR_RNDN);
  }

  for(size_t i=0;i<ncoeffs;++i){
    coeffs[i]=mpfr_get_d(m[i],MPFR_RNDN);
    mpfr_clear(m[i]);
  }
  mpfr_clear(a);
  mpfr_clear(b);
  mpfr_clear(c);

  delete[] m;
}