/*
    Copyright (C) 2021 Alessandro Languasco
*/

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

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


/*****************   Precomputations for the external C program   ****************/
 \\ \p needed to decide the precision  default (format, f)  for the floating point 
 
{ale_loggamma_dif_odd(x,y)=  local(vec, aoverq, S, oneminusaoverq, minusaoverq, q, 
g, a, u , coeff, r1, m, maxindex, datavector, minutes, millisec, abslogaoverq, abslogoneminusaoverq,
seconds, elaptimeprecomp, aleloggammafile, valueatonehalf, logoneminusaoverq, start, nhalf, bound, name, 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 loggamma-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_aleloggamma_dif_odd";
namefile = Strprintf("%s-%012d-%012d.res", name, x, y); 

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

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


maxindex = m+10;
datavector = vector(maxindex,i, zeta(i+1)/(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 : loggamma(a/q) - loggamma(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));
        logoneminusaoverq = log(oneminusaoverq);
        abslogoneminusaoverq = abs(logoneminusaoverq); 
        if (aoverq > 0.5,  
		                 \\ a/q >1/2 
 			r1 = ceil( ( coeff + abslogaoverq ) / abslogoneminusaoverq ); \\ number of summands   
 			start = oneminusaoverq;
 			correction = logoneminusaoverq + 2*Euler*oneminusaoverq
 			,
		                 \\ a/q <1/2 
 			r1 = ceil( ( coeff + abslogoneminusaoverq ) / abslogaoverq ); \\ number of summands   
            start = minusaoverq;
            correction = abslogaoverq + 2*Euler*minusaoverq;
        );   
        
        square=sqr(start);
        u = start*square; \\ starting point for computing 
        forstep (k=3, r1,  2,
 			  S+= u * datavector[k-1]; 
 			  u*= square;
 			  );   	
 		S = 2*S;
        S = correction + S; 	    \\ correction values
 		filewrite(aleloggammafile, S);
		a=(a*g)%q;  \\a=g^k%q;     
      );

elaptimeprecomp=gettime();
fileclose(aleloggammafile); 

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

---------

? ale_loggamma_dif_odd(0,305741)
Precomputation loggamma-values for q = 305741 and saving on file
Precomputation time (I/O included): 0 min, 1 sec, 426 millisec


----

? ale_loggamma(0,6766811) 
 
----

? ale_loggamma(0,10000019) 

-------

 ? ale_loggamma(0,28227761)  
*******/