12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849 |
- #include "coefficients.hpp"
- Coefficients::Coefficients(){
- coeffs=new mpfr_t[ncoeffs];
- for(size_t i=0;i<ncoeffs;++i) mpfr_init2(coeffs[i],4096);
- mpfr_t a,b,c;
- mpfr_init2(a,prec);
- mpfr_init2(b,prec);
- mpfr_init2(c,prec);
- //C(0,0)=0
- mpfr_set_zero(coeffs[pos(0,0)],0);
- //C(1,0)=-1
- mpfr_set_si(coeffs[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(coeffs[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,coeffs[pos(0,j-1)],4,MPFR_RNDN);//a=4*C(0,j-1)
- mpfr_sub(a,a,coeffs[pos(0,j-2)],MPFR_RNDN);//a=4*C(0,j-1)-C(0,j-2)
- mpfr_mul_ui(b,coeffs[pos(1,j-1)],2,MPFR_RNDN);//b=2*C(1,j-1)
- mpfr_sub(coeffs[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,coeffs[pos(i,j-1)],4,MPFR_RNDN);//a=4*C(i,j-1)
- mpfr_sub(a,a,coeffs[pos(i,j-2)],MPFR_RNDN);//a=4*C(i,j-1)-C(i,j-2)
- mpfr_sub(a,a,coeffs[pos(i-1,j-1)],MPFR_RNDN);//a=4*C(i,j-1)-C(i,j-2)-C(i-1,j-1)
- mpfr_sub(coeffs[pos(i,j)],a,coeffs[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,coeffs[pos(j-1,j-1)],2,MPFR_RNDN);//a=2*C(j-1,j-1)
- mpfr_sub(coeffs[pos(j-1,j)],a,coeffs[pos(j-2,j-1)],MPFR_RNDN);
- }
- mpfr_clear(a);
- mpfr_clear(b);
- mpfr_clear(c);
- }
|