/*
    Copyright (C) 2020 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_new(x,y)= local(vec, aoverq, S, oneminusaoverq, q, nhalf, harmonicvector,
bound, g, a, u , coeff, r1, m, maxindex, datavector, minutes, millisec, logaoverq, logoneminusaoverq,
seconds, elaptimeprecomp, Sfile, name, namefile, square, correctionvalue); 

default(realprecision , 38); 
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];


/* initialization */
 
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 problems with the last interval


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

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

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

m=128; \\ binary precision required for our algorithm
coeff = (m+2)*log(2)*1.0; \\ decimal precision  
\\print(maxindex); 


maxindex = m+10;

\\ REMARK: in fact here we could generate just the values at even integers ; but this is performed just once !
harmonicvector = vector(maxindex);
harmonicvector[1]=1.0; \\psi(2)+gamma
for (i=2, maxindex,harmonicvector[i]=harmonicvector[i-1]+(1.0/i));
\\ harmonicvector = vector(maxindex,i, psi(i+1)+Euler);
datavector = vector(maxindex,i, (zeta(i+1)*harmonicvector[i] + zetahurwitz(i+1,1,1))/(i+1)) ;  \\initialisation of the needed coefficients

elaptimeprecomp = 0;
gettime();

a= g^x%q;	   
\\valueatonehalf = sqr(log(Pi))+Pi^2/12 - gamma1 - sqr(Euler);   \\  = 2S(1/2) REMARK: \\  q is prime => a/q <> 1/2 
for(k=x, y,  \\ S = - (S(a/q) + S(1-a/q))   
        S=0;
		aoverq = a/q*1.0; \\ already sorted for the final summation  in external program
        
		\\ if (aoverq == 0.5,  print("1/2");  S = valueatonehalf, \\ S(1/2) + S(1/2) 
		
		\\  q is prime => a/q <> 1/2 
		oneminusaoverq  = 1 - aoverq;
		logaoverq = log(aoverq);
		logoneminusaoverq = log(oneminusaoverq);
		if (aoverq > 0.5,  \\ a/q >1/2, (1-a/q) <1/2 \\ equations (45)-(46) of the paper
 
			r1 = ceil( ( coeff + abs(logaoverq) ) / abs(logoneminusaoverq) ); \\ number of summands  
			square = sqr(oneminusaoverq); \\ (1-a/q)^2; starting point for computing (1-a/q)^k
			correctionvalue =  sqr(logoneminusaoverq), \\ correction value for this case
			
			\\ a/q <1/2, (1-a/q) >1/2  \\ equations (45)-(46) of the paper
 
			r1 = ceil( ( coeff + abs(logoneminusaoverq) ) / abs(logaoverq) ); \\ number of summands  
			square = sqr(aoverq);  \\ (a/q)^2; starting point for computing (a/q)^k
			correctionvalue = sqr(logaoverq) ; \\ correctionvalue for this case
			);  
           
           u= square;
           forstep (k=2, r1, 2, \\ computation of S(a/q) + S(1-a/q) with equations (45)-(46) of the paper;  
 			  S += u  * datavector[k-1];
 			  u *= square; 
 			 );
        S = 2*S; 
		S = 2*S + correctionvalue;      
		S = -S;  \\ change sign to compare with previous implementations
		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");	

} 
 
/*******************
 
gp2c-run -pmy_ -g -W  DIFprecS-new-v5.gp 

? DIFprecS_new(0,10007)
Precomputation (decimated in frequency) S-values for q = 10007 and saving on file
Precomputation time (I/O included): 0 min, 0 sec, 53 millisec

---
? DIFprecS_new(0,305741)
Precomputation (decimated in frequency) S-values for q = 305741 and saving on file
Precomputation time (I/O included): 0 min, 1 sec, 466 millisec


---
? DIFprecS_new(0,6766811)
Precomputation (decimated in frequency) S-values for q = 6766811 and saving on file
Precomputation time (I/O included): 0 min, 32 sec, 660 millisec

---
? DIFprecS_new(0,10000019)
Precomputation (decimated in frequency) S-values for q = 10000019 and saving on file
Precomputation time (I/O included): 0 min, 48 sec, 908 millisec


---
? DIFprecS_new(0,28227761)
Precomputation (decimated in frequency) S-values for q = 28227761 and saving on file
Precomputation time (I/O included): 2 min, 18 sec, 815 millisec

---
? DIFprecS_new(0,193894451)
Precomputation (decimated in frequency) S-values for q = 193894451 and saving on file
Precomputation time (I/O included): 15 min, 41 sec, 47 millisec

precendente versione; 
? DIFprecS_new(0,193894451)
Precomputation (decimated in frequency) S-values for q = 193894451 and saving on file
Precomputation time (I/O included): 246 min, 0 sec, 507 millisec
 
 quella con serie/integrale era molto piu' lenta 
 
 -----
 
? DIFprecS_new(0,212634221)
Precomputation (decimated in frequency) S-values for q = 212634221 and saving on file
Precomputation time (I/O included): 17 min, 33 sec, 951 millisec


precendente versione; 
? DIFprecS_new(0,212634221)
Precomputation (decimated in frequency) S-values for q = 212634221 and saving on file
 Precomputation time (I/O included): 275 min, 5 sec, 520 millisec
 

versione con integrale/serie impiegava circa 21200 minuti; circa 1206 volte piu' veloce;
per 500 000 valori versione con integrale/serie: circa 100 minuti
ossia circa 5000 valori al minuto
questa versione fa circa 6 075 000 valori al minuto !
 
*******************/