/*
    Copyright (C) 2019 Alessandro Languasco
*/

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

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

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


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

void
my_precpsi(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_minutes = gen_0, my_millisec = gen_0, my_seconds = gen_0, my_elaptimeprecomp = gen_0, my_psifile = 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 psi-values for q = %Ps and saving on file\n", my_q);
  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_psi-");
  my_namefile = Str(mkvecn(5, my_name, my_x, strtoGENstr("-"), my_y, strtoGENstr(".res")));
  my_psifile = stoi(gp_fileopen(GENtostr_unquoted(my_namefile), "w"));
  if (gequal0(my_x))
    gp_filewrite(gtos(my_psifile), 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_psifile), GENtostr_unquoted(gpsi(my_aoverq, 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_psifile));
  /* 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_precpsi();
   printf("x = ");
  scanf("%ld",&a);
  printf("y = ");
  scanf("%ld",&b);
  x = stoi(a);
  y = stoi(b);

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

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