/*
Copyright (C) 2024 Alessandro Languasco
*/

/**********
****** COMPUTATION OF sum log p/ (p^s);  Re(s)> 1
****** FOR Euler's constants in AP (L.-Moree)
**********/


global(SQUARE_FREE_UP_TO_J);

\\ Routine to compute the sum p <= P of  log(p) / (p^s - 1)); to get the truncated log derivative zeta(s)
{sum_fin(P, s) = local(fin_sum);
 

fin_sum = 0; \\ computation of the truncated sum over  log(p) / (p^s - 1)
forprime (p = 2, P, 
					fin_sum +=  log(p) /(p^s - 1);	
					
		); 		

return(fin_sum);
}


{sum_primes(s, prec=19) = local( K, P, err1, moebius_vector, 
Smallprimes_sum, fin_sum_stimesj, toterr, num_square_free,
elaptimecomp, seconds, minutes, millisec, Sj, frac_part_sigma,
err2, accuracy, J, stimesj, real_s, int_part_sigma,
aux, zeta_derlog_value, zetatrunc_derlog);


if (floor(prec)<> prec, prec = floor(prec);
						print("prec forced to be an integer");
						);
if (prec < 5, error("prec must be an integer >= 5"));

default(realprecision, prec+10);        \\ set working precision;  

real_s = real(s);
if (real_s <=1, error ("Re(s) must be >1")); 
int_part_sigma = floor(real_s);
frac_part_sigma = real_s - int_part_sigma;

print("****** COMPUTATION OF  sum log p/ (p^s);  Re(s)> 1  *******");
print("****** WITH A PRECISION OF AT LEAST ", prec," DECIMAL DIGITS *******");

\\ settings of the main parameters for a precision of about 100 decimal digits 

gettime();  
P = 3200;  

accuracy = 10^(-prec-10);

\\K = 26;
K = int_part_sigma + 2;
err1 =  log(P) / ((K +frac_part_sigma - 1)*(P^(K-1+frac_part_sigma))*(P-1));	;
until (err1 < accuracy/2, 
		K += 1;
		err1 = log(P) / ((K +frac_part_sigma - 1)*(P^(K-1+frac_part_sigma))*(P-1));	
		);

\\J = 26
J = 2; 
err2 = P^2/(P^((J+1)*int_part_sigma)*(P^(J+1)-1.0)); 
until (err2 < accuracy/2, 
		J += 1;
		err2 = P^2/(P^((J+1)*int_part_sigma)*(P^(J+1)-1.0)); 
		);	

print("-------------------------- ");
print ("Accuracy parameters: ");
if (prec == 19 , print("default precision is with prec = 19"));
print("prec = ", prec, "; P = ", P, "; K = ", K, "; J = ", J);
print("Accuracy is 10^(-",prec+10,")");
print("err1 = ", err1); 
print("err2 = ", err2);  
toterr = err1+err2;
\\print("toterr = ", toterr); 

print("-------------------------- ");


\\ PRECOMPUTATIONS 


moebius_vector = vector(J);
moebius_vector[1] = 1;
for (j = 2, J, moebius_vector[j] = moebius(j) * 1.0);

num_square_free = hammingweight(moebius_vector);			
SQUARE_FREE_UP_TO_J = vector(num_square_free);
SQUARE_FREE_UP_TO_J[1] = 1;
i=2;
for (j = 2, J, 
			if (moebius_vector[j] <> 0,
					SQUARE_FREE_UP_TO_J[i] = j;
					i+=1
					);
	);
	
Sj = 0;    
\\for (j=1, J ,
foreach(SQUARE_FREE_UP_TO_J,j,
			stimesj = j*s; 
			fin_sum_stimesj = sum_fin(P, stimesj);
			zeta_derlog_value = zetahurwitz(stimesj,1,1)/zeta(stimesj);
			zetatrunc_derlog = zeta_derlog_value + fin_sum_stimesj;								
		    Sj += moebius_vector[j]* zetatrunc_derlog;
	);

Sj = -Sj;
	
Smallprimes_sum = 0;
forprime(p = 2, P,
			Smallprimes_sum  +=  log(p)/(p^s);
		);

print("-------------------------- ");
print("          RESULTS          ");
elaptimecomp=gettime();
print("-------------------------- ");

\\results = vector(2);

aux = Sj + Smallprimes_sum;
\\results[1] = aux;
\\results[2] = toterr;
\\	print("Smallprimes_sum = ", Smallprimes_sum);

print("sum_ p log p/(p^s) (s = ", s, ") = ", aux);
print("toterr = ", toterr);	


print("-------------------------- ");
print("******  COMPUTATION TIME   ********");
seconds=floor(elaptimecomp/1000)%60;
minutes=floor(elaptimecomp/60000);
millisec=elaptimecomp - minutes*60000 - seconds*1000; 
print("Computation time: ", minutes, " min, ", seconds, " sec, ", millisec, " millisec");

print("******  END PROGRAM   ********");
\\return(results);
}
						
/************

? sum_primes(2,50)
****** COMPUTATION OF  sum log p/ (p^s);  Re(s)> 1  *******
****** WITH A PRECISION OF AT LEAST 50 DECIMAL DIGITS *******
-------------------------- 
Accuracy parameters: 
prec = 50; P = 3200; K = 18; J = 6
Accuracy is 10^(-60)
err1 = 3.83627296350736212922875764752309417968907465779715945511153 E-64
err2 = 2.52435489670723777731753214358888519888395307280942320034525 E-67
-------------------------- 
-------------------------- 
          RESULTS          
-------------------------- 
sum_ p log p/(p^s) (s = 2) = 0.493091109368764462197826205056491258055588126367846383144333
toterr = 3.83879731840406936700607517966668306488795861086996887831188 E-64
-------------------------- 
******  COMPUTATION TIME   ********
Computation time: 0 min, 0 sec, 17 millisec
******  END PROGRAM   ********


*****************/