/*
    Copyright (C) 2020 Alessandro Languasco
*/

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

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


/*****************   Precomputations for the external C program   ****************/
 \\ \p needed to decide the precision  default (format, f)  for the floating point 
 
{precT_new(x,y)=  local(coeff, r1, maxindex, oneminusaoverq, m, aux, harmonicvector,
datavector, v, T, vec, aoverq, q, g, a, minutes, millisec, minusaoverq,
logaoverq, logoneminusaoverq, factor, correction, abslogaoverq, abslogoneminusaoverq,
seconds, elaptimeprecomp, Tfile, 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("Precomputation T-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_T-new";
namefile = Strprintf("%s-%012d-%012d.res", name, x, y); 

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

\\ print(q);print(g);
 
m=128; \\ binary precision required for our algorithm
coeff = (m+2)*log(2)*1.0; \\ decimal precision  
maxindex = m+20;
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));  \\initialisation of the needed coefficients

\\ valueatonehalf = logtwo * (logtwo + twicegamma); \\  q is prime => a/q <> 1/2 

gettime();

a=g^x%q;	 
for(k=x, y,  
		aoverq = (a/q)*1.0; \\ already sorted for the final summation  in external       
        T=0;  \\ T will contain :  - T(a/q) 
        \\ if (aoverq == 0.5, T = valueatonehalf, \\  q is prime => a/q <> 1/2 
	      
	      \\ a/q <>1/2
        minusaoverq = -aoverq;
        oneminusaoverq  = 1 - aoverq;
		logaoverq = log(aoverq);
		abslogaoverq= abs(logaoverq);
		logoneminusaoverq = log(oneminusaoverq); 
		abslogoneminusaoverq = abs(log(oneminusaoverq));
        if (aoverq > 0.5,  
		                 \\ a/q >1/2
 			aux =   (coeff - log(abslogoneminusaoverq))/abslogoneminusaoverq;  \\ for the number of summands  
            factor = oneminusaoverq;
            correction = 0.,  	    \\ correction values  
            
		                    \\ a/q <1/2 
 			aux =   (coeff -log(abslogaoverq))/abslogaoverq; \\ for the number of summands   
            factor = minusaoverq;
            correction =  - logaoverq/aoverq; \\ correction values  
           );   
        
        r1 = 1+ ceil ( aux + log(log(aux))); \\ number of summands  
        v = factor;   
        for (k=2, r1, \\ computation of T(a/q);
 			  T+= v * datavector[k-1]; 
 			  v*= factor;
             );   	  
        T= correction + T;
        T= -T; \\changing sign (to be coherent with other definition of T previously used)
		filewrite(Tfile, T);
		a=(a*g)%q;  \\a=g^k%q;     
      );

elaptimeprecomp=gettime();
fileclose(Tfile); 

/* 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  precT-new-v3.gp
  
? precT_new(0,10007)
Precomputation T-values for q = 10007 and saving on file
Precomputation time (I/O included): 0 min, 0 sec, 232 millisec

? precT_new(0,305741)
Precomputation T-values for q = 305741 and saving on file
Precomputation time (I/O included): 0 min, 6 sec, 897 millisec

? precT_new(0,6766811)
Precomputation T-values for q = 6766811 and saving on file
Precomputation time (I/O included): 2 min, 37 sec, 531 millisec

? precT_new(0,10000019)
Precomputation T-values for q = 10000019 and saving on file
Precomputation time (I/O included): 3 min, 37 sec, 546 millisec

? precT_new(0,28227761)
Precomputation T-values for q = 28227761 and saving on file
Precomputation time (I/O included): 10 min, 15 sec, 839 millisec

?  precT_new(0,193894451)
Precomputation T-values for q = 193894451 and saving on file
Precomputation time (I/O included): 69 min, 27 sec, 456 millisec

? precT_new(0,212634221)
Precomputation T-values for q = 212634221 and saving on file
Precomputation time (I/O included): 76 min, 46 sec, 848 millisec

 
*******/