/*
    Copyright (C) 2021 Alessandro Languasco
*/

/****************   	A. LANGUASCO  		********************  *******/

 /* Here we compute the values of psi(a/q)+psi(1-a/q) */


/*****************   Precomputations for the external C program   ****************/
 \\ \p needed to decide the precision  default (format, f)  for the floating point 
 
{ale_psi_dif_even(x,y)=  local(vec, aoverq, S, oneminusaoverq, minusaoverq, q,
g, a, u , coeff, r1, m, maxindex, datavector, minutes, millisec, abslogaoverq, abslogoneminusaoverq,
seconds, elaptimeprecomp, alepsifile, valueatonehalf, name, nhalf, bound, namefile, square, correction); 



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("Precomputation psi-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_alepsi_dif_even";
namefile = Strprintf("%s-%012d-%012d.res", name, x, y); 

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

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


maxindex = m+10;
datavector = vector(maxindex,i, zeta(i+1)) ;  \\initialisation of the needed coefficients
 
gettime();

a=g^x%q;		
for(k=x, y,  
		aoverq = (a/q)*1.0; \\ already sorted for the final summation  in external       
        S=0;  \\ S will contain :  psi(a/q)+psi(1-a/q) 
        \\if (aoverq == 0.5, print("1/2"); S = valueatonehalf, \\  q is prime => a/q <> 1/2 
	      \\ a/q <>1/2
        minusaoverq = -aoverq;
        oneminusaoverq  = 1 - aoverq;
        abslogaoverq = abs(log(aoverq));	
        abslogoneminusaoverq = abs(log(oneminusaoverq)); 
        if (aoverq > 0.5,  
		                 \\ a/q >1/2 
 			r1 = ceil( ( coeff + abslogaoverq ) / abslogoneminusaoverq ); \\ number of summands     
 			square = sqr(oneminusaoverq);
 			correction = - 2*Euler - 1/oneminusaoverq
 			,
		                 \\ a/q <1/2 
 			r1 = ceil( ( coeff + abslogoneminusaoverq ) / abslogaoverq ); \\ number of summands     
            square = sqr(aoverq);
            correction = - 2*Euler + 1/minusaoverq;
        );   
 
        u = square; \\ squareing point for computing 
        forstep (k=2, r1,  2,
 			  S+= u * datavector[k]; 
 			  u*= square;
 			  );   
 		
 		S = 2*S;	
        S = correction - S; 	    \\ correction values
 		filewrite(alepsifile, S);
		a=(a*g)%q;  \\a=g^k%q;     
      );

elaptimeprecomp=gettime();
fileclose(alepsifile); 

/* 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  ale_psi_dif_even.gp
  
     
  
? ale_psi_dif_even(0,10007)
Precomputation psi-values for q = 10007 and saving on file
Precomputation time (I/O included): 0 min, 0 sec, 51 millisec

----

? ale_psi_dif_even(0,305741)
Precomputation psi-values for q = 305741 and saving on file
Precomputation time (I/O included): 0 min, 1 sec, 418 millisec
 
*******/