/*
    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) + S(1-a/q)) */
 
 /*****************   Precomputations for the external C program   ****************/
 \\ \p needed to decide the precision  default (format, f)  for the floating point 
 
{DIFprecS(x,y)= local(vec, aoverq, S, eps, ord, oneminusaoverq, q, nhalf, bound, g, a, minutes, millisec, seconds, elaptimeprecomp, Sfile, tab, name, namefile); 

default(format , f);  \\ print numbers with floating point notation (not E notation)

vec = readvec("./primroot.res"); \\getting q and g from primroot.res

q=vec[1];
g=vec[2];

\\ print(q);print(g);

/* initialization */
eps=1/100.; \\ parameter needed to switch between the two definitions of S(x)+S(1-x)
tab=sumnuminit();  \\ needed for the sumnum function
 
print("Precomputation (decimated in frequency) S-values for q = ",q, " and saving on file"); 

nhalf=(q-1)/2;
bound= nhalf-1;
if (x<0, x=0); \\ needed to avoid problem with the first interval
if (x>=y, error("first exponent should be less than the second"));
if (y>bound, y=bound); \\ needed to avoid problem with the last interval


name="./precomp_S-DIF-";
namefile = Strprintf("%s-%012d-%012d.res", name, x, y); 

Sfile = fileopen(namefile, "w");  
if (x==0,filewrite(Sfile,q)); 

gettime();

a= g^x%q;	
for(k=x, y,  \\ S = - (S(a/q) + S(1-a/q)) 
		aoverq = a/q; \\ already sorted for the final summation  in external program
		oneminusaoverq = 1 - aoverq;
				
		\\ DECIMATION IN FREQUENCY (with Dilcher sum formula about log(Gamma_1) or Deninger integral formula ) 
                
        ord=min(aoverq, oneminusaoverq); \\rate of decay for the integral formula of S-DIF
		if (ord < eps, 
				\\ if ord is small, it uses the series def. of S-DIF; otherwise it uses the integral def. of S-DIF
		        \\ using log1p in the second log gives large differences; bug ?
				\\ sqr(x) is a bit faster than x^2
		    S = sumnum(n=1, -sqr(log1p(aoverq/n)) - sqr(log(1-aoverq/n)) - 2*log(n)*log1p(-sqr(aoverq/n)) , tab) - sqr(log(aoverq)),
		    
		    S = (-2)*intnum(t=0,[+oo,ord],(-3 + exp(-t) + exp(aoverq*t) + exp(oneminusaoverq*t))*(Euler+log(t))/(t*(exp(t)-1)));
			);
		filewrite(Sfile, S);
		a=(a*g)%q;  \\a=g^k%q;    
	);

elaptimeprecomp=gettime();
fileclose(Sfile);  

/* print computation time */
seconds=floor(elaptimeprecomp/1000)%60;
minutes=floor(elaptimeprecomp/60000);
millisec=elaptimeprecomp- minutes*60000 - seconds*1000;
\\print(elaptimeprecomp); 
print("Precomputation time (I/O included): ", minutes, " min, ", seconds, " sec, ", millisec, " millisec");	

} 
 
/************************************  
 
DA FARE SU OGNI MACCHINA (le librerie create da gp2c dipendono dalla macchina che si usa)
 air:Desktop languasc$ gp2c-run -pmy_ -g -W DIFprecS.gp 
 
per creare le librerie .so e il file .gp2c
 
*******/