/*
    Copyright (C) 2021 Alessandro Languasco
*/

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

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


/*****************   Precomputations for the external C program   ****************/
 \\ \p needed to decide the precision  default (format, f)  for the floating point 
 
{ale_psi(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, namefile, start, 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");
 
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>q-2, y=q-2); \\ needed to avoid problem with the last interval


name="./precomp_alepsi";
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*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) 
        \\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     
 			start = oneminusaoverq;
 			correction = - Euler
 			,
		                 \\ a/q <1/2 
 			r1 = ceil( ( coeff + abslogoneminusaoverq ) / abslogaoverq ); \\ number of summands     
            start = minusaoverq;
            correction = - Euler + 1/minusaoverq;
        );   
 
        u = start; \\ starting point for computing 
        for (k=2, r1,  
 			  S+= u * datavector[k-1]; 
 			  u*= start;
 			  );   	
        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.gp
  
    
  Questo script e' un poco piu' lento di psi implementato in parigp
  (il che vuol dire che facendolo per bene, e' perlomeno della stessa
  velocita')
  
  ? precpsi(0,10007)
Precomputation psi-values for q = 10007 and saving on file
Precomputation time (I/O included): 0 min, 0 sec, 128 millisec
 
 ? ale_psi(0,10007)
Precomputation psi-values for q = 10007 and saving on file
Precomputation time (I/O included): 0 min, 0 sec, 134 millisec

----

? ale_psi(0,305741)
Precomputation psi-values for q = 305741 and saving on file
Precomputation time (I/O included): 0 min, 3 sec, 844 millisec

? precpsi(0,305741)
Precomputation psi-values for q = 305741 and saving on file
Precomputation time (I/O included): 0 min, 3 sec, 812 millisec

----

? ale_psi(0,6766811)
Precomputation psi-values for q = 6766811 and saving on file
Precomputation time (I/O included): 1 min, 26 sec, 849 millisec

? precpsi(0,6766811)
Precomputation psi-values for q = 6766811 and saving on file
Precomputation time (I/O included): 1 min, 25 sec, 46 millisec

----

? ale_psi(0,10000019)
Precomputation psi-values for q = 10000019 and saving on file
Precomputation time (I/O included): 2 min, 10 sec, 899 millisec

? precpsi(0,10000019)
Precomputation psi-values for q = 10000019 and saving on file
Precomputation time (I/O included): 2 min, 6 sec, 766 millisec

-------

? ale_psi(0,28227761)
Precomputation psi-values for q = 28227761 and saving on file
Precomputation time (I/O included): 6 min, 7 sec, 542 millisec

? precpsi(0,28227761)
Precomputation psi-values for q = 28227761 and saving on file
Precomputation time (I/O included): 5 min, 53 sec, 579 millisec
 
*******/