/*-*- compile-command: "/usr/bin/gcc -c -o DIFprecS.gp.o -O3 -Wall -fno-strict-aliasing -fomit-frame-pointer -fPIC -I\"/usr/local/include\" DIFprecS.gp.c && /usr/bin/gcc -o DIFprecS.gp.so -bundle -undefined dynamic_lookup -O3 -Wall -fno-strict-aliasing -fomit-frame-pointer -fPIC DIFprecS.gp.o "; -*-*/
#include <pari/pari.h>
/*
GP;install("my_init_DIFprecS","v","init_DIFprecS","./DIFprecS.gp.so");
GP;install("my_DIFprecS","vD0,G,D0,G,p","DIFprecS","./DIFprecS.gp.so");
*/
void my_init_DIFprecS(void);
void my_DIFprecS(GEN my_x, GEN my_y, long prec);
/*End of prototype*/

void
my_init_DIFprecS(void)	  /* void */
{
  pari_sp ltop = avma;
  avma = ltop;
  return;
}

static GEN
my_anon_0(GEN my_n, GEN my_aoverq, GEN my_boverq, long prec)
{
  pari_sp ltop = avma;
  GEN p1 = gen_0;
  p1 = gsub(gsub(gsub(gneg(gsqr(glog(gaddsg(1, gdiv(my_aoverq, my_n)), prec))), gmul(gmulsg(2, glog(my_n, prec)), gsub(glog(gaddsg(1, gdiv(my_aoverq, my_n)), prec), gdiv(my_aoverq, my_n)))), gsqr(glog(gaddsg(1, gdiv(my_boverq, my_n)), prec))), gmul(gmulsg(2, glog(my_n, prec)), gsub(glog(gaddsg(1, gdiv(my_boverq, my_n)), prec), gdiv(my_boverq, my_n))));
  p1 = gerepilecopy(ltop, p1);
  return p1;
}

static GEN
wrap_my_anon_0(void * _cargs, GEN my_n)
{
  GEN _args = (GEN) _cargs;
  GEN _res = my_anon_0(my_n, gel(_args,1), gel(_args,2), gtos(gel(_args,3)));
  return _res;
}

/*
Copyright (C) 2019 Alessandro Languasco
*/

/****************   	A. LANGUASCO  		******************** 
************* COMPUTATION OF THE EULER KRONECKER CONSTANTS MOD q(PRIME) *******/


/* Here we compute the values of -S(a/q) */

/*****************   Precomputations for the external C program   ****************/
/* \p needed to decide the precision  default (format, f)  for the floating point  */

void
my_DIFprecS(GEN my_x, GEN my_y, long prec)	  /* void */
{
  pari_sp ltop = avma;
  GEN my_vec = gen_0, my_aoverq = gen_0, my_boverq = gen_0, my_S = gen_0, my_q = gen_0, my_nhalf = gen_0, my_bound = gen_0, my_g = gen_0, my_a = gen_0, my_b = gen_0, my_C1 = gen_0, my_minutes = gen_0, my_millisec = gen_0, my_seconds = gen_0, my_elaptimeprecomp = gen_0, my_Sfile = gen_0, my_tab = gen_0, my_name = gen_0, my_namefile = gen_0, my_f = pol_x(fetch_user_var("f"));
  default0("format", GENtostr_unquoted(my_f));
  /* print numbers with floating point notation (not E notation) */
  
  my_vec = gp_readvec_file("./primroot.res");
  /*getting q and g from primroot.res */
  
  my_q = gcopy(gel(my_vec, 1));
  my_g = gcopy(gel(my_vec, 2));
  /* print(q);print(g); */
  
  /* initialization */
  my_C1 = negr(strtor("0.145631690967353449721172751749802638275472676668675905198013119482802867143022969756173856489688029208", prec));
  my_tab = sumnuminit(NULL, prec);
  /* needed for the sumnum function */
  
  pari_printf("Precomputation (decimated in frequency) S-values for q = %Ps and saving on file\n", my_q);
  my_nhalf = gdivgs(gsubgs(my_q, 1), 2);
  my_bound = gsubgs(my_nhalf, 1);
  if (gcmpgs(my_x, 0) < 0)
    my_x = gen_0;
  /* needed to avoid problem with the first interval */
  if (gcmp(my_x, my_y) >= 0)
    pari_err(e_MISC, "first exponent should be less than the second");
  if (gcmp(my_y, my_bound) > 0)
    my_y = gcopy(my_bound);
  /* needed to avoid problem with the last interval */
  
  
  my_name = strtoGENstr("./DIFprecomp_S-");
  my_namefile = Strprintf("%s-%012d-%012d.res", mkvec3(my_name, my_x, my_y));
  my_Sfile = stoi(gp_fileopen(GENtostr_unquoted(my_namefile), "w"));
  if (gequal0(my_x))
    gp_filewrite(gtos(my_Sfile), GENtostr_unquoted(my_q));
  gettime();
  my_a = gmod(gpow(my_g, my_x, prec), my_q);
  my_b = gmod(gpow(my_g, gadd(my_x, my_nhalf), prec), my_q);
  /* For decimation in frequency of S; we'll need just the even positions of fft(S) */
  {
    pari_sp btop = avma;
    GEN my_k = gen_0;
    for (my_k = gcopy(my_x); gcmp(my_k, my_y) <= 0; my_k = gaddgs(my_k, 1))
    {
      my_aoverq = gdiv(my_a, my_q);
      /* already sorted for the final summation  in external program */
      my_boverq = gdiv(my_b, my_q);
      /* DECIMATION IN FREQUENCY  (with two sums of series (sum a_n) + (sum b_n) )  */
      /* Saq = sumnum(n=2, -(log(1+aoverq/n))^2-2*log(n)*(log(1+aoverq/n)-aoverq/n),tab)-(log(1+aoverq))^2 -C1*aoverq-(log(aoverq))^2; */
      /* Sbq = sumnum(n=2, -(log(1+boverq/n))^2-2*log(n)*(log(1+boverq/n)-boverq/n),tab)-(log(1+boverq))^2 -C1*boverq-(log(boverq))^2; */
      /* filewrite(Sfile, Saq + Sbq);  \\ S(g^k/q)+S(g^(k+nhalf)/q)    */
      
      /* DECIMATION IN FREQUENCY (with sum (a_n +b_n) ) */
      my_S = gsub(gsub(gsub(sumnum(mkvec3(my_aoverq, my_boverq, stoi(prec)), wrap_my_anon_0, gen_1, my_tab, prec), gmul(my_C1, gadd(my_aoverq, my_boverq))), gsqr(glog(my_boverq, prec))), gsqr(glog(my_aoverq, prec)));
      gp_filewrite(gtos(my_Sfile), GENtostr_unquoted(my_S));
      my_a = gmod(gmul(my_a, my_g), my_q);
      /*a=g^k%q;  */
      my_b = gmod(gmul(my_b, my_g), my_q);
      if (gc_needed(btop, 1))
        gerepileall(btop, 6, &my_aoverq, &my_boverq, &my_S, &my_a, &my_b, &my_k);
    }
  }
  /*b=(g^(k+nhalf))%q;  */
  my_elaptimeprecomp = stoi(gettime());
  gp_fileclose(gtos(my_Sfile));
  /* print computation time */
  my_seconds = gmodgs(gfloor(gdivgs(my_elaptimeprecomp, 1000)), 60);
  my_minutes = gfloor(gdivgs(my_elaptimeprecomp, 60000));
  my_millisec = gsub(gsub(my_elaptimeprecomp, gmulgs(my_minutes, 60000)), gmulgs(my_seconds, 1000));
  /*print(elaptimeprecomp);  */
  pari_printf("Precomputation time (I/O included): %Ps min, %Ps sec, %Ps millisec\n", my_minutes, my_seconds, my_millisec);
  avma = ltop;
  return;
}