self_avoiding_polygon/exact/main.cpp

#include <iomanip>
#include <cilk/cilk.h>
#include <cilk/cilk_api.h>
#include "config.hpp"
#include "coeffs.hpp"
#include "polygon_generator.hpp"
#include "results.hpp"
#include "flint/fmpq_poly.h"

ResultsReducer cilk_result;

void box(string str){
  size_t w=str.size();
  cout<<"\u250c";
  for(size_t i=0;i<w+2;++i) cout<<"\u2500";
  cout<<"\u2510"<<endl;
  cout<<"\u2502 "<<str<<" \u2502"<<endl;
  cout<<"\u2514";
  for(size_t i=0;i<w+2;++i) cout<<"\u2500";
  cout<<"\u2518"<<endl;
}

void disp_info(){
  box("Compute Fp's series of Self Avoiding Polygon");
  cout<<endl;
  cout<<" Maximal length is "<<max_len<<endl;
  cout<<" Workers number is "<<__cilkrts_get_nworkers()<<endl;
}

void treat(Polygon* P){
  size_t l=P->length;
  size_t i=l/2-1;
  ++cilk_result.numbers(i);
  fmpq_poly_t fp;
  fmpq_poly_init(fp);
  P->fp(fp);
  cilk_result.add_fp(i,fp);
  delete P;
}

int main(){
  disp_info();
  compute_coeffs();
  size_t nb[max_len/2];
  PolygonGenerator gen;
  size_t total=0;
  Polygon* P;
  while(gen.next()){
    P=new Polygon;
    gen.set(*P);
    cilk_spawn treat(P);
    ++total;
  }
  cilk_sync;

  //! Mutiplicative coeffcient for the sum of Fp 
  fmpq_t c;
  fmpq_init(c);
  
  //! Numerator and denominator of the coefficient c
  fmpz_t d,n;
  fmpz_init(n);
  fmpz_init(d);

  //! Init denominator to 256
  fmpz_set_si(d,256);
  for(size_t i=1;2*i<max_len;++i){

    //! Length of the SAP
    Int l=2*i+2;
    cout<<endl<<"=== Length "<<l<<" ==="<<endl;

    //! Set numerator of c
    fmpz_set_si(n,2*l);
    //! Set c
    fmpq_set_fmpz_frac(c,n,d);

    //! Retrieve Fp sum in res
    fmpq_poly_t res;
    fmpq_poly_init(res);
    fmpq_poly_set(res,cilk_result.get_fp(i));

    //! Multiply res by c
    fmpq_poly_scalar_mul_fmpq(res,res,c);

    //! Display res
    char* str=fmpq_poly_get_str_pretty(res,"x");
    
    cout<<" > number : "<<cilk_result.numbers(i)<<endl;
    cout<<" > coeff  : "<<str<<endl; 

    delete str;
    
    //! Update denominator
    fmpz_mul_si(d,d,16);
  }
  cout<<endl<<">>> Total of SAP computed : "<<total<<endl;
  cout<<endl;
}