/*
    Copyright (C) 2019 Alessandro Languasco
*/

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

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

static GEN
my_anon_0(GEN my_m, GEN my_x, long prec)
{
  pari_sp ltop = avma;
  GEN p1 = gen_0;
  p1 = gsub(gdiv(glog(gadd(my_x, my_m), prec), gadd(my_x, my_m)), gdiv(glog(gaddsg(1, my_m), prec), gaddsg(1, my_m)));
  p1 = gerepilecopy(ltop, p1);
  return p1;
}

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

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


GEN
my_DilcherTtab(GEN my_x, GEN my_tabstandalone, long prec)
{
  pari_sp ltop = avma;
  GEN my_T = gen_0, my_m = gen_0;
  /* for the standalone precomputations */
  my_T = sumnum(mkvec2(my_x, stoi(prec)), wrap_my_anon_0, gen_0, my_tabstandalone, prec);
  /*tabstandalone is passed  */
  my_T = gerepilecopy(ltop, my_T);
  return my_T;
}

static GEN
my_anon_1(GEN my_m, GEN my_aoverq, long prec)
{
  pari_sp ltop = avma;
  GEN p1 = gen_0;
  p1 = gsub(gdiv(glog(gadd(my_aoverq, my_m), prec), gadd(my_aoverq, my_m)), gdiv(glog(gaddsg(1, my_m), prec), gaddsg(1, my_m)));
  p1 = gerepilecopy(ltop, p1);
  return p1;
}

static GEN
wrap_my_anon_1(void * _cargs, GEN my_m)
{
  GEN _args = (GEN) _cargs;
  GEN _res = my_anon_1(my_m, gel(_args,1), gtos(gel(_args,2)));
  return _res;
}

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

void
my_precT(GEN my_x, GEN my_y, long prec)	  /* void */
{
  pari_sp ltop = avma;
  GEN my_vec = gen_0, my_aoverq = gen_0, my_q = gen_0, my_g = gen_0, my_a = gen_0, my_C1 = gen_0, my_minutes = gen_0, my_millisec = gen_0, my_seconds = gen_0, my_elaptimeprecomp = gen_0, my_Tfile = 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); */
  
  pari_printf("Precomputation T-values for q = %Ps and saving on file\n", my_q);
  /* initialization */
  my_C1 = negr(strtor("0.145631690967353449721172751749802638275472676668675905198013119482802867143022969756173856489688029208", prec));
  my_tab = sumnuminit(NULL, prec);
  /* needed for the sumnum function */
  
  
  if ((gcmpgs(my_x, 0) < 0) || (gcmp(my_y, gsubgs(my_q, 2)) > 0))
    pari_err(e_MISC, "the exponents are from 0 to q-2");
  if (gcmp(my_x, my_y) > 0)
    pari_err(e_MISC, "first exponent should be less than the second");
  my_name = strtoGENstr("./precomp_T-");
  my_namefile = Str(mkvecn(5, my_name, my_x, strtoGENstr("-"), my_y, strtoGENstr(".res")));
  my_Tfile = stoi(gp_fileopen(GENtostr_unquoted(my_namefile), "w"));
  if (gequal0(my_x))
    gp_filewrite(gtos(my_Tfile), GENtostr_unquoted(my_q));
  gettime();
  my_a = gmod(gpow(my_g, my_x, prec), my_q);
  {
    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  */
      gp_filewrite(gtos(my_Tfile), GENtostr_unquoted(sumnum(mkvec2(my_aoverq, stoi(prec)), wrap_my_anon_1, gen_0, my_tab, prec)));
      my_a = gmod(gmul(my_a, my_g), my_q);
      if (gc_needed(btop, 1))
        gerepileall(btop, 3, &my_aoverq, &my_a, &my_k);
    }
  }
  /*a=g^k%q;  */
  my_elaptimeprecomp = stoi(gettime());
  gp_fileclose(gtos(my_Tfile));
  /* 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;
}






int main() {
  long a, b;
  GEN x, y;
  pari_init(40000000,0);  /* memory size, primes precomputed*/
  my_init_precT();
   printf("x = ");
  scanf("%ld",&a);
  printf("y = ");
  scanf("%ld",&b);
  x = stoi(a);
  y = stoi(b);

  my_precT(x,y,38);
 
  pari_close();
  return 0;
}

/*
clang -o precT.exe -g -O3 -Wall -fomit-frame-pointer -fno-strict-aliasing -fPIC -I"/usr/include" precT.c -lc -lm -L/usr/lib -lpari
*/