/*
    Copyright (C) 2021 Alessandro Languasco
*/

/****************   	A. LANGUASCO  		******************** 
************* COMPUTATION OF Sr - LvR  *******/  

\\ Global variables:

      
global(tab);   \\global variable used to initialise sumnum
global(pir); 
global(defaultprecision);

 /************* COMPUTATION OF Sr - LvR; q(PRIME) ******* 
*************    FOR q1<=q<=q2  *********/

  
{Srall(rstart,rstop,q1,q2,Pbound,defaultprecision)=
local(minutes, millisec, seconds, elaptimefinalcomp, 
s1, s2, s3, s4, errs1, errs2, errs3, errs4, S, errS, s3expos, s4expos,
qminusone, v, qfrac, recqminusone, gpminusone,
logp, qpow, gp, fp,  qpowminusone, gppowminusone, gppow, gp1pow,
Pbound1, emptyset, V, d, divbound, err1, err2, flaggp1, flaggp2,
recomp, recompfile, Srqfile, S1qfile, S2qfile, S3qfile, S4qfile, S5qfile, S6qfile
);

 
	 
\\ minutes,millisec,seconds: used just to compute the elapsed computation time; local variables

\\ defaultprecision: used to fix the precision used in the computations; global variable


print("************ A. LANGUASCO *************");
print("********* COMPUTATION of Sr - LvR  **********");
print("******* VALUES IN ONE INTERVAL ********"); 
 
\\ precision setting

default(realprecision,defaultprecision);       
default(format,f);   
default(primelimit, Pbound);  

Pbound1=Pbound+1; \\ for the errors

if (rstop < rstart, print("error: 1 <= r <= 6. END PROGRAM");return);
if (rstop > 6 || rstart < 1,  print("error: 1 <= r <= 6. END PROGRAM");return);
print("**********  IT WORKS ONLY for q = 1 mod 2r **********");

\\print("case r = ", r);

q1=nextprime(q1);
q2=precprime(q2);
if (q1>q2, print("error: no odd primes in this interval. END PROGRAM");return);

if (q1 < 3, q1=3);
print("q-Interval=[",q1,",",q2,"]");
print("r-Interval=[",rstart,",",rstop,"]");

pir=primepi(q2)-primepi(q1)+1; 

tab = sumnuminit(); \\ needed for the sumnum functions in the errors	

gettime(); 
 

\\ uncomment to obtain the outputfile
S1qfile = fileopen("S1qfile.csv", "w");filewrite(S1qfile,"q;r;Sr;errSr");
S2qfile = fileopen("S2qfile.csv", "w");filewrite(S2qfile,"q;r;Sr;errSr");
S3qfile = fileopen("S3qfile.csv", "w");filewrite(S3qfile,"q;r;Sr;errSr");
S4qfile = fileopen("S4qfile.csv", "w");filewrite(S4qfile,"q;r;Sr;errSr");
S5qfile = fileopen("S5qfile.csv", "w");filewrite(S5qfile,"q;r;Sr;errSr");
S6qfile = fileopen("S6qfile.csv", "w");filewrite(S6qfile,"q;r;Sr;errSr");
recompfile = fileopen("Srecomp.txt", "w"); 


gettime();
print("Starting computation of Sr(q)");
print("internal sums over primes p up to ", Pbound);
print("---------");
print("*********** THE RESULTS ARE IN THE OUPUT FILES *********** ");
print("---------");

  
emptyset = Set([]);
recomp = [emptyset, emptyset, emptyset, emptyset, emptyset, emptyset];

forprime(q=q1,q2,
	print("q = ", q);
	
    errs1=vector(rstop,rows,0);
    errs2=vector(rstop,rows,0);
	errs3=vector(rstop,rows,0);
	errs4=vector(rstop,rows,0);
	err1=0;
	err2=0;
	flaggp1=vector(rstop,rows,0);
	flaggp2=vector(rstop,rows,0); 
  
	s1=vector(rstop,rows,0);	
	s2=vector(rstop,rows,0);
	s3=vector(rstop,rows,0);
	s4=vector(rstop,rows,0);
	S = vector(rstop,rows,0);
	\\print(S);
	errS = vector(rstop,rows,0);  
  
  
	qminusone = q-1;	
	recqminusone = (1.0)/qminusone;
	qfrac = q*recqminusone;    
	v=[qminusone, factor(qminusone)];

	divbound = numdiv(qminusone);
	\\print(divbound); 
	 
	s3expos = [emptyset, emptyset, emptyset, emptyset, emptyset, emptyset];
	s4expos = [emptyset, emptyset, emptyset, emptyset, emptyset, emptyset];
	
 		
	forprime(p=2, qminusone,  \\ to keep q=3 in this case
					logp = log(p);
					qpowminusone = (1.0)*p^qminusone;   
					qpow = p^qminusone * p; 
					fp = znorder(Mod(p,q),v); 
					
				    for (r = rstart, rstop,	
				    \\print("r=",r);	
				        if(qminusone%(2*r) == 0,
				           if( r == 1, gp = fp, gp = fp/gcd(fp,r)); 				
					       if( gp == 1, flaggp1[r] = 1;
					       				s1[r] += ( 1 / (qpowminusone -1) - qfrac/(qpow - 1)) * logp , \\ gp =  1; if not 
					     	   gp == 2, flaggp2[r] = 1;
					     	   			s2[r] += logp / (p^2 -1),  \\ gp =  2; if not  
					     	    gp >= 3, 
								 if ( setsearch(s3expos[r], gp) == 0, s3expos[r]=concat(s3expos[r],gp); s3expos[r]=Set(s3expos[r]));
								 gpminusone = (1.0)*(gp-1) ; gppowminusone = (1.0)*p^gpminusone; gppow = gppowminusone * p;
								 s3[r] += ( gpminusone / (gppowminusone -1) - gp/(gppow - 1))*logp;
					   			 if ( (gp >= 4) && (gp % 2 == 0),  \\ if gp is also gp >=4 and even				   			 		
								       if ( setsearch(s4expos[r], gp) == 0, s4expos[r]=concat(s4expos[r],gp); s4expos[r]=Set(s4expos[r]));
					   			 	   gp1pow = sqrt(gppow); 
					   			 	   s4[r] +=  gp1pow * logp / (gppow -1) 
					   			 	);
					   			); 	
					       );
				       );	   	 
		);

	forprime(p=qminusone+3, Pbound, \\ from q+2
					logp = log(p);
					qpowminusone = (1.0)*p^qminusone;   
					qpow = p^qminusone * p; 
					fp = znorder(Mod(p,q),v); 
					
				    for (r = rstart, rstop,	
				    \\print("r=",r);	
				        if(qminusone%(2*r) == 0,
				           if( r == 1, gp = fp, gp = fp/gcd(fp,r)); 					
					       if( gp == 1, flaggp1[r] = 1;
					       				s1[r] += ( 1 / (qpowminusone -1) - qfrac/(qpow - 1)) * logp , \\ gp =  1; if not 
					     	   gp == 2, flaggp2[r] = 1;
					     	   			s2[r] += logp / (p^2 -1),  \\ gp =  2; if not 
					     	   gp >= 3, 
								 if ( setsearch(s3expos[r], gp) == 0, s3expos[r]=concat(s3expos[r],gp); s3expos[r]=Set(s3expos[r]));
								 gpminusone = (1.0)*(gp-1) ; gppowminusone = (1.0)*p^gpminusone; gppow = gppowminusone * p;
								 s3[r] += (gpminusone / (gppowminusone -1) - gp/(gppow - 1))*logp;
					   			 if ( (gp >= 4) && (gp % 2 == 0),  \\ if gp is also gp >=4 and even				   			 		
								       if ( setsearch(s4expos[r], gp) == 0, s4expos[r]=concat(s4expos[r],gp); s4expos[r]=Set(s4expos[r]));
					   			 	   gp1pow = sqrt(gppow); 
					   			 	   s4[r] +=  gp1pow * logp / (gppow -1) 
					   			 	);
					   			); 	
					       );
				       ); 	 
 			);
 
 	err1 = sumnum(n=Pbound1, log(n)*( qminusone / (n^qminusone*1.0 -1) - q/(n^q*1.0 - 1) ), tab);
	err2 = sumnum(n=Pbound1, log(n)/(n^2*1.0 -1), tab);
	
	for (r = rstart, rstop,	
	 if(qminusone%(2*r) == 0,
		S[r] = qminusone*s1[r] - s2[r] + s3[r] - s4[r];    
	 	
	 	if (flaggp1[r]  == 1, errs1[r] = err1);
	 	if (flaggp2[r]  == 1, errs2[r] = err2); 
	 		
    	for (j = 1, #s3expos[r], 
     				V = s3expos[r];
     				d = V[j];
      				errs3[r] +=  sumnum(n=Pbound1, log(n)*( (d-1)/(n^(d-1)*1.0 -1) - d/(n^d*1.0 -1) ), tab);
					);			
	    for (j = 1, #s4expos[r], 	
    				V = s4expos[r];	
    				d = V[j];	
					errs4[r] += sumnum(n=Pbound1, n^(d/2)*log(n) / (n^d*1.0 -1) , tab)
					);			
					 
		errS[r] = abs(errs1[r]) + abs(errs2[r]) + abs(errs3[r]) + abs(errs4[r]) ;	 
		   
 		\\print("S",r,"(",q,") = ", S[r] );  
 		\\print("errS",r,"(",q,") = ", errS[r]); 
 		\\if (abs(S[r]) < 10*errS[r], recomp[r]= concat(recomp[r],q); print("*********  recompute this case !!!"));
 		\\print("---------");
 		
 		if (abs(S[r]) < 10*errS[r], recomp[r]= concat(recomp[r],q);
 									filewrite1(recompfile, r);
									filewrite1(recompfile,"; "); 
							 		filewrite(recompfile, q); 
							 		fileflush(recompfile); 	
 									);
 
 		if(r==1, Srqfile = S1qfile,
 		  r==2, Srqfile = S2qfile,		
 		  r==3, Srqfile = S3qfile,		
 		  r==4, Srqfile = S4qfile,		
 		  r==5, Srqfile = S5qfile,	
 		  r==6, Srqfile = S6qfile	 					    					    					   
 		  );
 		filewrite1(Srqfile,q);
 		filewrite1(Srqfile,";"); 
 		filewrite1(Srqfile,r);
 		filewrite1(Srqfile,";");  
 		filewrite1(Srqfile, S[r]); 
 		filewrite1(Srqfile,";");
 		filewrite(Srqfile, errS[r]); 
 		fileflush(Srqfile); 						
		);
	); 
	 
 
);

print("---------");
fileclose(S1qfile); 
fileclose(S2qfile); 
fileclose(S3qfile); 
fileclose(S4qfile); 
fileclose(S5qfile);  
fileclose(S6qfile);  
elaptimefinalcomp=gettime();

for (r = rstart, rstop,	
 filewrite(recompfile, "------------------");
 filewrite(recompfile, "you should recompute these cases: ");
 filewrite1(recompfile, "r = ");
 filewrite1(recompfile, r);
 filewrite1(recompfile,"; ");
 filewrite1(recompfile,"Pbound = 10^");
 filewrite(recompfile, floor(log(Pbound)/log(10)));
 filewrite(recompfile, recomp[r]); 
 filewrite(recompfile, "----------");
 fileflush(recompfile); 
 );	
fileclose(recompfile);  
 		 
seconds=floor(elaptimefinalcomp/1000)%60;
minutes=floor(elaptimefinalcomp/60000);
millisec=elaptimefinalcomp- minutes*60000 - seconds*1000;
 
print("Sr(q) computation time (output time included): ", minutes, " min, ", seconds, " sec, ", millisec, " millisec");

print("******  END PROGRAM   ********"); 
}
 


 
 
/**** RESULTS ****

 gp2.13.1 and gp2c0.0.12 compiled by myself, see below 
 

------- 
 
 optiplex
------- 

 
 ? init_Srall(); Srall(1,6,3,3000,10^8,19)
************ A. LANGUASCO *************
********* COMPUTATION of Sr - LvR  **********
******* VALUES IN ONE INTERVAL ********
**********  IT WORKS ONLY for q = 1 mod 2r **********
q-Interval=[3,2999]
r-Interval=[1,6]
Starting computation of Sr(q)
internal sums over primes p up to 100000000
---------
*********** THE RESULTS ARE IN THE OUPUT FILES ***********
---------
q = 3
.....
q = 2887
q = 2897
q = 2903
q = 2909
q = 2917
q = 2927
q = 2939
q = 2953
q = 2957
q = 2963
q = 2969
q = 2971
q = 2999
---------
Sr(q) computation time (output time included): 2342 min, 17 sec, 797 millisec
******  END PROGRAM   ********

---------
*******/