/*
Copyright (C) 2024 Alessandro Languasco
*/

/**********
****** COMPUTATION OF Euler's constants in AP (L.-Moree)
**********/

\\ Global variables
global(CHI_VECTOR_NON_PRINC, DFACTORS_POW_MATR, PPOWERS_MATR, NUM_PRIME_FACTORS_D, VALCHI_CONJ_A_MATRIX, VALCHI_P_MATRIX);
global(DFACTORS, TWOPII, S5VECTOR, S2VECTOR, GAMMA1VECTOR, LOG_DFACTORS, LOG_P, LOG_P_OVER_P); 
global(SQUARE_FREE_UP_TO_J, LOG_DER_AT_1, SJCHI_VECTOR, NUMPRIMES_UP_TO_P, ONE_OVER_P_MINUS_ONE);
global(ONE_OVER_P_D_MINUS_ONE, ORD_CHI_J_MATRIX, CHI_J_MATRIX, ISPRIMEDIVISOR_OF_D, NUMPRIMES_UP_TO_D, PRIMES_UP_TO_D);

\\ Routine to compute the value of L'/L(s,chi) when chi is a Dirichlet character (also imprimitive)

{L_derlog(L, chi, G, d, s) = local(conductor, Gstar, chistar, corr_chi_imprim, val_chistar_d1,
ud1, L_value, L_prime_value, L_derlog_value, d1);

L_value = lfun(L, s); \\ value of L at s for the attached primitive character		
L_prime_value = lfun(L, s, 1); \\ value of L'at s for the attached primitive character		
L_derlog_value = L_prime_value / L_value; \\ value of logarithmic derivatie at s for the attached primitive character	

conductor = zncharconductor(G, chi);
				
if (conductor < d, \\ chi is imprimitive
			\\ computing the correction sum for the logarithmic derivative
			[Gstar,chistar] = znchartoprimitive(G, chi); 
			\\ chistar is the primitive character such that chi is induced by chistar
	 											  
			corr_chi_imprim = 0.;
			for(i=1, NUM_PRIME_FACTORS_D, 
						 d1 = DFACTORS[i];
	 					 if(conductor % d1 <> 0, \\ d1 not divides both d and conductor of chi and conductor <> 1
				 				ud1 = chareval(Gstar, chistar, d1); 
	 							val_chistar_d1 = exp(TWOPII * ud1); \\ value of chistar(d1)								
		 						corr_chi_imprim += val_chistar_d1 * LOG_DFACTORS[i] / (DFACTORS_POW_MATR[i,s] - val_chistar_d1) 
		 						\\ we can use precomputed powers of d1^s since s is k*j, 2 <= k < = K, 1 <= j <= J
	 							);
	 						);

						L_derlog_value += corr_chi_imprim ;
				);

return(L_derlog_value);
}

\\ Routine to compute the sum p <= P of chi(p) log(p) / (p^s - chi(p))); to get the truncated log derivative L(s, chi)
\\ Computes directly the difference of the finite sums for chi^j and chi^{kj}

{chi_sum_fin_delta(P, k, chi_kj, chi_j, G, d, s) = local(fin_sum_chi, up_j, val_chi_j_p, 
val_chi_kj_p, den, p);
 
fin_sum_chi = 0.; \\ computation of the truncated sum over chi(p) log(p) / (p^s - chi(p)

\\ we can use precomputed powers of p^s since s is k*j, 2 <= k < = K, 1 <= j <= J

\\ working on (p,d) == 1 
for (i = 1, NUMPRIMES_UP_TO_D,
			\\if ( d % p <> 0,	\\ p does not divide d	
			if ( ISPRIMEDIVISOR_OF_D[i] == 0,
						p = PRIMES_UP_TO_D[i];
						up_j = chareval(G, chi_j , p ); 
						val_chi_j_p = exp(TWOPII * up_j);

						den = denominator (up_j); 
						val_chi_kj_p = (val_chi_j_p)^(k%den);

						fin_sum_chi += (val_chi_kj_p  /(PPOWERS_MATR[i,s] - val_chi_kj_p)  - val_chi_j_p   /(PPOWERS_MATR[i,s] - val_chi_j_p)) * LOG_P[i];	

						\\ we can use precomputed powers of p^s since s is k*j, 2 <= k < = K, 1 <= j <= J
						);
		); 		 
\\ working on (p,d) == 1 
for (i = NUMPRIMES_UP_TO_D + 1, NUMPRIMES_UP_TO_P,
						p = PPOWERS_MATR[i,1];
						up_j = chareval(G, chi_j , p ); 
						val_chi_j_p = exp(TWOPII * up_j); 
						
						den = denominator (up_j); 
						val_chi_kj_p = (val_chi_j_p)^(k%den); 
						
						fin_sum_chi += (val_chi_kj_p  /(PPOWERS_MATR[i,s] - val_chi_kj_p)  - val_chi_j_p   /(PPOWERS_MATR[i,s] - val_chi_j_p)) * LOG_P[i];		
						\\ we can use precomputed powers of p^s since s is k*j, 2 <= k < = K, 1 <= j <= J									
		); 		

return(fin_sum_chi);
}

{AP_gamma_all(D, prec=19) = local(results, logtwo, phid, resfile_total,
toterrs, dhalf, i, output);

resfile_total = fileopen("gamma_AP-results_TOTAL.csv", "a"); 
\\filewrite(resfile,"d; a; gamma_da; accuracy"); 

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

D = floor(D);
if (D<3, error("D must be >=3"));
results = matrix(D,D-1);
toterrs = matrix(D,D-1);

logtwo = log(2);
results[1,1] = Euler;
results[2,1] = Euler + logtwo;
filewrite1(resfile_total,1); 
filewrite1(resfile_total,";");
filewrite1(resfile_total,1); 
filewrite1(resfile_total,";");	
filewrite1(resfile_total,results[1,1]); 
filewrite(resfile_total,";");	

filewrite1(resfile_total,2); 
filewrite1(resfile_total,";");
filewrite1(resfile_total,1); 
filewrite1(resfile_total,";");	
filewrite1(resfile_total,results[2,1]); 
filewrite(resfile_total,";");	
fileflush(resfile_total);


for (d=3, D,
	phid = eulerphi(d);
	dhalf = d\2;
	if (d%2 == 0 && (dhalf % 2) == 1,
			print("----------------------------");
			print("----------------------------");
			print("For d = ", d, " we can re-use previous computations");
			print("----------------------------");
			print("----------------------------");
		 	results[d,1] = results[dhalf,1];
			toterrs[d,1] = toterrs[dhalf,1];
			results[d, 2+ dhalf] = results[dhalf, 2] + logtwo;
			toterrs[d, 2+ dhalf] = toterrs[dhalf, 2] ;
			for (a = 3, dhalf,
				if (gcd(a,dhalf) == 1,
					if (a % 2 ==1,
						results[d,a] = results[dhalf, a];
						toterrs[d,a] = toterrs[dhalf, a],
						\\ else
						results[d,a+dhalf] = results[dhalf, a];
						toterrs[d,a+dhalf] = toterrs[dhalf, a]
						)
					)
				)		
			,
		\\else	dhalf is even or d is odd
		i = 0;
		output = matrix(phid, 2);
		output = AP_gamma(d,prec);
		for (a = 1, d,	
			if (gcd(a,d) ==1,	
				 i+=1;	
				 results[d,a] = output[i,1];
				 toterrs[d,a] = output[i,2];
				 );
			);
		);
	\\saving on file	
	for (a = 1, d,	
		if (gcd(a,d) ==1,	
			filewrite1(resfile_total,d); 
			filewrite1(resfile_total,";");
			filewrite1(resfile_total,a); 
			filewrite1(resfile_total,";");	
			filewrite1(resfile_total,results[d,a]); 
			filewrite1(resfile_total,";");	
			default(format,e); 
			default(realprecision,5);  
			filewrite(resfile_total,toterrs[d,a]); 
			default(format,f); 
			default(realprecision, prec+10);   
			fileflush(resfile_total);
			);
		);	
	);
fileclose(resfile_total); 
}

{AP_gamma(d, prec=19) = local(dminusone, A, K, P, err1, G, moebius_vector, 
Smallprimes_sum, chi, chi_to_kj, dcoprimes, toterr, phid,
L1, L2, elaptimecomp, seconds, minutes, millisec, ua, Sjchi, ord_chi_j,
Sjkchi, err2, val_chi_p, up, factorp, i, accuracy, J, KtimesJ, ktimesj, 
phidminusone, aux, test_value, sum_chars, numdcoprimes, p, chi_to_j, L, j,
sum_constants, L1_derlog_value, L2_derlog_value, num_square_free,
fin_sum_chi_delta, val_div_d, firstterm, Lderlog_at_1, results);

if (d > 1000, error("d must be <= 1000"));
if (d <= 2, error("d must be >= 3"));

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

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

dminusone = d-1;
phid = eulerphi(d);

print("****** COMPUTATION OF Euler's constants in AP *******");
print("****** WITH A PRECISION OF AT LEAST ", prec + 10," DECIMALS *******");

\\ settings of the main parameters for a precision of about 100 decimals 

gettime();  
A = 1000;  


P = A*d; \\ is P in the file
if (P > 9600, P = 9600);

accuracy = 10^(-prec-10);

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

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

print("-------------------------- ");
print ("Accuracy parameters: ");
if (prec == 19 , print("default precision is with prec = 19"));
print("prec = ", prec, "; A = ", A, "; 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("-------------------------- ");

\\ setting the group
G = znstar(d, 1); 
TWOPII = 2*Pi*I;

\\ PRECOMPUTATIONS 

KtimesJ = floor(K*J);
moebius_vector = vector(J);

moebius_vector[1] = 1;
for (j = 2, J, 
			moebius_vector[j] = moebius(j) * 1.0; 
			);
			
\\ needed since gp2c does not support forsquarefree
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
					);
	);


NUMPRIMES_UP_TO_P = primepi(P);
PPOWERS_MATR = matrix(NUMPRIMES_UP_TO_P, KtimesJ); 
LOG_P = vector(NUMPRIMES_UP_TO_P);
LOG_P_OVER_P = vector(NUMPRIMES_UP_TO_P);
ONE_OVER_P_MINUS_ONE = vector(NUMPRIMES_UP_TO_P);

i=1;
forprime ( p=2, P, 
			factorp = p;
			PPOWERS_MATR[i,1]= factorp;
			LOG_P[i] = log(factorp);
			LOG_P_OVER_P[i] = LOG_P[i]/factorp;
			ONE_OVER_P_MINUS_ONE[i] = 1.0/(factorp - 1.0);
			for(k=2, KtimesJ, 
					PPOWERS_MATR[i,k] = PPOWERS_MATR[i,k-1]*factorp );
			i += 1					
	);

\\ generating the list of integers <= d-1 coprime with d

dcoprimes = vector(phid);
i=1;
dcoprimes[1]=1;
for(m=2, dminusone, 
	if (gcd(m,d) == 1, 
		i+=1;
		dcoprimes[i]=m
		)
	);
numdcoprimes = #dcoprimes;	

\\ generating the non-principal Dirichlet characters mod d
\\ and the orders of his powers
phidminusone = phid - 1;
CHI_VECTOR_NON_PRINC = vector(phidminusone);
CHI_J_MATRIX = matrix(phidminusone, J);
ORD_CHI_J_MATRIX = matrix(phidminusone, J);
i=1;
foreach(dcoprimes, m,
	\\ m = 1 corresponds to the principal character mod d	
	if( m > 1, 
		CHI_VECTOR_NON_PRINC[i]=znconreychar(G,m);
		CHI_J_MATRIX[i,1] = CHI_VECTOR_NON_PRINC[i];
		ORD_CHI_J_MATRIX[i,1] = charorder(G,CHI_J_MATRIX[i,1]);
		foreach(SQUARE_FREE_UP_TO_J,j, \\ needed since gp2c does not support forsquarefree
			CHI_J_MATRIX[i,j] =  charpow(G, CHI_J_MATRIX[i,1], j);

			ORD_CHI_J_MATRIX[i,j] =  ORD_CHI_J_MATRIX[i,1] / gcd(j, ORD_CHI_J_MATRIX[i,1]);
			);
		i+=1;
		);
	);

\\ DFACTORS contains the list of prime factors of d
\\ it is needed to handle the L functions associated to imprimitive characters

DFACTORS=factor(d)[,1]; 
NUM_PRIME_FACTORS_D = #DFACTORS;

DFACTORS_POW_MATR = matrix(NUM_PRIME_FACTORS_D, KtimesJ); 
LOG_DFACTORS = vector(NUM_PRIME_FACTORS_D);
ONE_OVER_P_D_MINUS_ONE = vector(NUM_PRIME_FACTORS_D);

for(i=1, NUM_PRIME_FACTORS_D, 
			p = DFACTORS[i];
			factorp = p*1.0;
			ONE_OVER_P_D_MINUS_ONE[i] = 1.0/(factorp - 1.0);
			DFACTORS_POW_MATR[i,1]= factorp;
			LOG_DFACTORS[i] = log(factorp);
			for(k=2, KtimesJ, 
					DFACTORS_POW_MATR[i,k] = DFACTORS_POW_MATR[i,k-1]*factorp )	
	);
	
\\ ISPRIMEDIVISOR_OF_D contains 1 or 0 depending if p<=d divides d or not	
NUMPRIMES_UP_TO_D = primepi(d);
ISPRIMEDIVISOR_OF_D = vector(NUMPRIMES_UP_TO_D);
PRIMES_UP_TO_D = vector(NUMPRIMES_UP_TO_D);
i = 1 ;
forprime (p = 2, d, 
			PRIMES_UP_TO_D[i] = p;
			if ( d % p == 0,
					ISPRIMEDIVISOR_OF_D[i] = 1 
					);
			i+=1
			);

\\ S5VECTOR stores the value of the sum over each non trivial character
S5VECTOR = vector(numdcoprimes); \\ initialised with zeros
SJCHI_VECTOR = vector(phidminusone);

for(m=1, phidminusone, \\ we run over the non principal characters
		chi = CHI_VECTOR_NON_PRINC[m];
		Sjchi = 0;    

		foreach(SQUARE_FREE_UP_TO_J,j, \\ needed since gp2c does not support forsquarefree

				chi_to_j = CHI_J_MATRIX[m,j];
				L2 = lfuncreate([G, chi_to_j]);
				ord_chi_j = ORD_CHI_J_MATRIX[m,j];
				Sjkchi = 0;
				
				for(k = 2, K,	

					if ( k % ord_chi_j <> 1,
						\\ for k % ord_chi_j == 1 we have chi^(kj) = chi^j and the summand is zero 
						ktimesj = k*j; 

						chi_to_kj = charpow(G, chi_to_j, k);
						
						L1 = lfuncreate([G, chi_to_kj]);
						L1_derlog_value = L_derlog(L1, chi_to_kj, G, d, ktimesj);
						 
						L2_derlog_value = L_derlog(L2, chi_to_j, G, d, ktimesj);
						 
						fin_sum_chi_delta = chi_sum_fin_delta(P, k, chi_to_kj, chi_to_j, G, d, ktimesj);
						 
						Sjkchi += L1_derlog_value - L2_derlog_value + fin_sum_chi_delta;
						);
					);
			 Sjchi += moebius_vector[j]* Sjkchi;
			);
		SJCHI_VECTOR[m] = Sjchi;
	);
	
VALCHI_CONJ_A_MATRIX = matrix(numdcoprimes,phidminusone);

i=0;
foreach(dcoprimes, a,
	i += 1 ;
	sum_chars = 0; 
	for(m=1, phidminusone, \\ we run over the non principal characters
		chi = CHI_VECTOR_NON_PRINC[m];
		ua = chareval(G, chi , a ); 
		
		VALCHI_CONJ_A_MATRIX[i,m] = exp(-TWOPII * ua); \\ conjugate character
		sum_chars += VALCHI_CONJ_A_MATRIX[i,m]* SJCHI_VECTOR[m];
	);
	S5VECTOR[i] = -sum_chars;
);	


S2VECTOR = vector(numdcoprimes); 


VALCHI_P_MATRIX = matrix(NUMPRIMES_UP_TO_P,phidminusone);
for(m=1, phidminusone, \\ we run over the nonprincipal characters
	chi = CHI_VECTOR_NON_PRINC[m];
 	j = 0;
 	\\ working on (p,d) == 1
	forprime (p = 2, d, 
			j+=1;
			if ( d % p  == 0 ,		
						VALCHI_P_MATRIX[j,m] = 0.
						, \\else
						up = chareval(G, chi , p ); 
						VALCHI_P_MATRIX[j,m] = exp(TWOPII * up)
				);
			);	
	forprime (p = d+1, P, 
			j+=1;
			up = chareval(G, chi , p ); 
			VALCHI_P_MATRIX[j,m] = exp(TWOPII * up);
			);
);

i=0;
foreach(dcoprimes, a,
	i += 1 ;
	Smallprimes_sum = 0;
	for(m=1, phidminusone, \\ we run over the nonprincipal characters
		aux = 0;
		j = 0;
		forprime(p = 2, P,
			j+=1;			
			val_chi_p = VALCHI_P_MATRIX[j,m];
			aux += val_chi_p * LOG_P_OVER_P[j] * ( val_chi_p/(p-val_chi_p) - ONE_OVER_P_MINUS_ONE[j] ) ;		
		    );
	   	Smallprimes_sum += VALCHI_CONJ_A_MATRIX[i,m] * aux;
	);
	S2VECTOR[i] = Smallprimes_sum;
);

GAMMA1VECTOR = vector(numdcoprimes); 

val_div_d = 0;
for(i=1, NUM_PRIME_FACTORS_D, 
			val_div_d += LOG_DFACTORS[i]*ONE_OVER_P_D_MINUS_ONE[i];
	);

LOG_DER_AT_1 = vector(phidminusone);
for(m=1, phidminusone,
	chi = CHI_VECTOR_NON_PRINC[m];
	L = lfuncreate([G, chi]); 
	LOG_DER_AT_1[m] = L_derlog(L, chi, G, d, 1); 
	);

i=0;
foreach(dcoprimes, a,
	i += 1 ;
	firstterm = 0;
	for(m=1, phidminusone, \\ we run over the nonprincipal characters

		Lderlog_at_1 = LOG_DER_AT_1[m];
	   	firstterm += VALCHI_CONJ_A_MATRIX[i,m]* Lderlog_at_1;	
	);
	GAMMA1VECTOR[i] = Euler + val_div_d + firstterm;
);


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

sum_constants = 0;
print("-------------------------- ");
i=0;
results = matrix(numdcoprimes,2);
foreach(dcoprimes, a,
	i += 1 ;
	aux = real(GAMMA1VECTOR[i] + S2VECTOR[i] + S5VECTOR[i])/phid;
	results[i,1] = aux;
	results[i,2] = toterr;
	print("gamma(",d,";",a,") = ", aux);
	sum_constants += aux;
	);
	


print("-------------------------- ");
test_value = Euler + val_div_d;
print("Consistency test = ", sum_constants - test_value );

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);
}
						
/************

Last login: Sun May  5 13:02:57 on ttys007
languasc@zygalski ~ % cd /Users/languasc/Documents/LANGUASCO/matematica/lavori/\[80\]EK_for_AP/programs
languasc@zygalski programs % gp2c-run -pmy_ -g -W gamma_arprog_final.gp
Warning:gamma_arprog_final.gp:164: variable undeclared
e
Warning:gamma_arprog_final.gp:167: variable undeclared
f
Reading GPRC: /Users/languasc/.gprc
GPRC Done.

                  GP/PARI CALCULATOR Version 2.15.5 (released)
          i386 running darwin (x86-64/GMP-6.3.0 kernel) 64-bit version
     compiled: Feb 24 2024, Apple clang version 15.0.0 (clang-1500.1.0.2.5)
                           threading engine: pthread
                 (readline v8.0 enabled, extended help enabled)

                     Copyright (C) 2000-2022 The PARI Group

PARI/GP is free software, covered by the GNU General Public License, and comes 
WITHOUT ANY WARRANTY WHATSOEVER.

Type ? for help, \q to quit.
Type ?18 for how to get moral (and possibly technical) support.

parisizemax = 2048000000, primelimit = 500000, nbthreads = 12
? AP_gamma_all(10,50)
****** COMPUTATION OF Euler's constants in AP *******
****** WITH A PRECISION OF AT LEAST 60 DECIMALS *******
-------------------------- 
Accuracy parameters: 
prec = 50; A = 1000; P = 3000; K = 18; J = 9
Accuracy is 10^(-60)
err1 = 2.43208529966564416997966094010728207007655846855676182216720 E-63
err2 = 5.16234958342639436398006301623350651924394121773093271293970 E-63
toterr = 7.59443488309203853395972395634078858932049968628769453510689 E-63
-------------------------- 
          RESULTS          
-------------------------- 
gamma(3;1) = 1.09904952586667653048446536830561011726877035116603294279006
gamma(3;2) = 0.0274722833689111758196693402380551660971342636852653818830593
-------------------------- 
Consistency test = -3.454467422037777850 E-77
-------------------------- 
******  COMPUTATION TIME   ********
Computation time: 0 min, 1 sec, 20 millisec
******  END PROGRAM   ********
****** COMPUTATION OF Euler's constants in AP *******
****** WITH A PRECISION OF AT LEAST 60 DECIMALS *******
-------------------------- 
Accuracy parameters: 
prec = 50; A = 1000; P = 4000; K = 17; J = 9
Accuracy is 10^(-60)
err1 = 0.0000000000000000000000000000000000000000000000000000000000000603622441923637424518716161439957706799684599998793940586093
err2 = 0.0000000000000000000000000000000000000000000000000000000000000000291038304567337036132812500000000000277555756156289135105908
toterr = 0.0000000000000000000000000000000000000000000000000000000000000603913480228204761554848973939957706799962155754950229721199
-------------------------- 
          RESULTS          
-------------------------- 
gamma(4;1) = 0.986722568313428628851628485266200068085358865038395262040658
gamma(4;3) = 0.283640277148049541172115726274378931032300605261783590885789
-------------------------- 
Consistency test = 0.000000000000000000000000000000000000000000000000000000000000
-------------------------- 
******  COMPUTATION TIME   ********
Computation time: 0 min, 1 sec, 22 millisec
******  END PROGRAM   ********
****** COMPUTATION OF Euler's constants in AP *******
****** WITH A PRECISION OF AT LEAST 60 DECIMALS *******
-------------------------- 
Accuracy parameters: 
prec = 50; A = 1000; P = 5000; K = 17; J = 9
Accuracy is 10^(-60)
err1 = 0.00000000000000000000000000000000000000000000000000000000000000139573607969757585515414239744172562330462885534523748423492
err2 = 0.000000000000000000000000000000000000000000000000000000000000000000524288000000000000000000000000000000053687091200000000000000
toterr = 0.00000000000000000000000000000000000000000000000000000000000000139626036769757585515414239744172562330468254243643748423492
-------------------------- 
          RESULTS          
-------------------------- 
gamma(5;1) = 0.608523951539242472079793258724297784485288468537460426713149
gamma(5;2) = -0.352151568727300949661662936463231782979398446526485211013255
gamma(5;3) = -0.0277297956309668265257700614505971705577472592709054495812843
gamma(5;4) = 0.750932555829083258364341662578480509975416911766983263165319
-------------------------- 
Consistency test = -0.000000000000000000000000000000000000000000000000000000000000000000000000000008636168555094444625
-------------------------- 
******  COMPUTATION TIME   ********
Computation time: 0 min, 3 sec, 735 millisec
******  END PROGRAM   ********
----------------------------
----------------------------
For d = 6 we can re-use previous computations
----------------------------
----------------------------
****** COMPUTATION OF Euler's constants in AP *******
****** WITH A PRECISION OF AT LEAST 60 DECIMALS *******
-------------------------- 
Accuracy parameters: 
prec = 50; A = 1000; P = 7000; K = 16; J = 8
Accuracy is 10^(-60)
err1 = 0.0000000000000000000000000000000000000000000000000000000000000355267311321500157764873863491435651956486076557726607870119
err2 = 0.0000000000000000000000000000000000000000000000000000000000000601812709779332688901710474186200635595013025160464830267913
toterr = 0.0000000000000000000000000000000000000000000000000000000000000957080021100832846666584337677636287551499101718191438138032
-------------------------- 
          RESULTS          
-------------------------- 
gamma(7;1) = 0.524815040692497356262522341149533289906194861404909678348647
gamma(7;2) = -0.205605438156889246453478698566563348079177033238215998096621
gamma(7;3) = -0.207379226204487017674963327989581998053903224541037016890610
gamma(7;4) = 0.394136290298642495869491773880949339848115670869076082320908
gamma(7;5) = 0.00220396249714958098054633164538511148407692480072827852585877
gamma(7;6) = 0.393363393950505242473285793869876657543032924908106106007483
-------------------------- 
Consistency test = -0.000000000000000000000000000000000000000000000000000000000000000000000000000025908505665283333876
-------------------------- 
******  COMPUTATION TIME   ********
Computation time: 0 min, 8 sec, 726 millisec
******  END PROGRAM   ********
****** COMPUTATION OF Euler's constants in AP *******
****** WITH A PRECISION OF AT LEAST 60 DECIMALS *******
-------------------------- 
Accuracy parameters: 
prec = 50; A = 1000; P = 8000; K = 16; J = 8
Accuracy is 10^(-60)
err1 = 0.00000000000000000000000000000000000000000000000000000000000000425772382599000392766964880487485037243624016100796507190878
err2 = 0.00000000000000000000000000000000000000000000000000000000000000710542735760100185871124267578125005293955920339377119177016
toterr = 0.0000000000000000000000000000000000000000000000000000000000000113631511835910057863808914806561004253757993644017362636789
-------------------------- 
          RESULTS          
-------------------------- 
gamma(8;1) = 0.854935464772769233543358720452450064927724612173143130188155
gamma(8;3) = -0.0798291608610312540562783915706500134207138863064139763276721
gamma(8;5) = 0.131787103540659395308269764813750003157634252865252131852503
gamma(8;7) = 0.363469438009080795228394117845028944453014491568197567213461
-------------------------- 
Consistency test = 0.000000000000000000000000000000000000000000000000000000000000
-------------------------- 
******  COMPUTATION TIME   ********
Computation time: 0 min, 3 sec, 180 millisec
******  END PROGRAM   ********
****** COMPUTATION OF Euler's constants in AP *******
****** WITH A PRECISION OF AT LEAST 60 DECIMALS *******
-------------------------- 
Accuracy parameters: 
prec = 50; A = 1000; P = 9000; K = 16; J = 8
Accuracy is 10^(-60)
err1 = 0.000000000000000000000000000000000000000000000000000000000000000655218019242319177505869043711966067036573789502399157961946
err2 = 0.00000000000000000000000000000000000000000000000000000000000000107931905547085803064583170332768337932861875862937140974660
toterr = 0.00000000000000000000000000000000000000000000000000000000000000173453707471317720815170074703964944636519254813177056770854
-------------------------- 
          RESULTS          
-------------------------- 
gamma(9;1) = 0.458393411675144887073401432170122456606603184021781783568418
gamma(9;2) = -0.439091624987433551772260723033844635656382105053552664408795
gamma(9;4) = 0.376643569772260748501338788635142398143701665457401484371480
gamma(9;5) = 0.0310377589755237254080223102800284801451355305044033869513626
gamma(9;7) = 0.264012544419270894909725147500345262518465501686849674850158
gamma(9;8) = 0.435526149380821002183907752991871321608380838234414659340492
-------------------------- 
Consistency test = 0.000000000000000000000000000000000000000000000000000000000000
-------------------------- 
******  COMPUTATION TIME   ********
Computation time: 0 min, 8 sec, 171 millisec
******  END PROGRAM   ********
----------------------------
----------------------------
For d = 10 we can re-use previous computations
----------------------------
----------------------------
? 

 


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