Replat Int (Int64) as fmpz in computation of coeff

exact/coeffs.cpp

@@ -6,55 +6,100 @@
 
 Reel coeffs[num_coeffs];
 fmpz_poly_t coeffs_poly[num_coeffs];
-Int coeffs_den;
+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);
   }
-  CoeffAnalytic ca[num_coeffs];
-  ca[pos(0,0)]={0,0};  //C_(0,0)=0
-  ca[pos(0,1)]={-1,0}; //C_(1,0)=-1
+  //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)
-  Rationnal r=0;
+  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){
-    Rationnal t(1,2*i-1);
-    r=r+t;
-    ca[pos(i,i)]={0,-4*r};
+    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);
   }
-  coeffs_den=ca[pos(max_ind_coeffs,max_ind_coeffs)].b.denominator();
+  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)
-    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;
+    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)
-      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;
+      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)
-    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;
-    }
+    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_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());
+ 
+    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);
+
   }
+
 }

exact/coeffs.hpp

@@ -5,6 +5,7 @@
 #include "config.hpp"
 #include "rationnal.hpp"
 #include "flint/fmpz.h"
+#include "flint/fmpq.h"
 #include "flint/fmpz_poly.h"
 
 
@@ -21,7 +22,7 @@ static const size_t num_coeffs=((max_ind_coeffs+1)*(max_ind_coeffs+2))/2;
 //! Array of coefficients
 extern Reel coeffs[num_coeffs];
 extern fmpz_poly_t coeffs_poly[num_coeffs];
-extern Int coeffs_den;
+extern fmpz_t coeffs_den;
 
 //! Return the indice of coefficient C_{i,j}
 size_t pos(size_t i,size_t j);
@@ -38,7 +39,7 @@ void compute_coeffs();
 //! Represent an analytic coefficient
 struct CoeffAnalytic{
   //! The coefficient is a-4b/pi
-  Rationnal a,b;
+  fmpq_t a,b;
 };
 
 //********************

exact/polygon.cpp

@@ -153,7 +153,8 @@ Polygon::fp(fmpq_poly_t res) const{
 
   //! M=4*coeffs_denI+CB
   //! U[i][0]=1;
-  fmpz_poly_set_si(cc,4*coeffs_den);
+  fmpz_mul_si(z,coeffs_den,4);
+  fmpz_poly_set_fmpz(cc,z);
   for(size_t i=0;i<ne;++i){
     fmpz_poly_add(fmpz_poly_mat_entry(M,i,i),fmpz_poly_mat_entry(M,i,i),cc);
     fmpz_poly_set_si(fmpz_poly_mat_entry(U,i,0),1);
@@ -179,7 +180,7 @@ Polygon::fp(fmpq_poly_t res) const{
   fmpq_poly_set_fmpz_poly(res,fmpz_poly_mat_entry(R,0,0));
 
   //! Compute the multiplicative scalar factor
-  fmpz_set_si(z,coeffs_den*4);
+  fmpz_mul_si(z,coeffs_den,4);
   fmpz_pow_ui(z,z,ne-1);
   fmpz_mul_si(z,z,4*s);