/****************   	A. LANGUASCO & A. ZACCAGNINI 	*********************** 
************* COMPUTATION OF THE MERTENS CONSTANTS MOD q; 3 <= q<= 100 ******* 
************* WITH A PRECISION OF AT LEAST 100 DECIMAL DIGITS ***************** 
************* Verification of eq. (25) of the paper ***************************
*******************************************************************************/

/******************************************************************************   
*****   The results are collected in a matrix named MCmatrix.gp    ************
*****   every row contains q, a, C(q,a), number of correct digits  ************
*****   for a fixed q>=3, the starting row for the C(q,a)s is detected ********
*****   with the q-component of the vector position which is precomputed ******
*****           The value corresponding to q=2 is inserted as      ************
*****           C(2,1) = 2 exp(-gamma); gamma = Euler constant	   ************
*****       The results of this control are printed on a file      ************ 
*****                named MCcheckresults.txt                      ************
*******************************************************************************/

{checksMC() = local(MC, position, div, q2, q1, j, n, ndiv,a, factors,nfactors, 
l, i0, phiq1,phiq2,digitprecision,flagerror,counter);

position=vector(100);
position[1]=0;
position[2]=1;

\\ REMARK: C(2,1) = 2 exp(-gamma); gamma=Euler constant

position[3]=2;
for (q=4,100, position[q]=position[q-1]+eulerphi(q-1));
\\write("~/MCcheckresults.txt",position);

\\ input of the constants
MC=read("~/MCmatrix.gp");



default(realprecision,120);   
\\ REMARK: in this verification we admit a precision of 100 decimal digits  
digitprecision=100;

print("****** A. LANGUASCO & A. ZACCAGNINI *******");
write("~/MCcheckresults.txt","****** A. LANGUASCO & A. ZACCAGNINI *******");

print("****** VERIFICATION OF EQ. (25) FOR 4 <= q <= 100 *******");
write("~/MCcheckresults.txt","****** VERIFICATION OF EQ. (25) FOR 4 <= q <= 100 *******");
print("****** REMARK: C(2,1) = 2 exp(-gamma); gamma=Euler constant *******");
write("~/MCcheckresults.txt","****** REMARK: C(2,1) = 2 exp(-gamma); gamma=Euler constant *******");

print("******************* START CHECK ********************");
write("~/MCcheckresults.txt","******************* START CHECK ********************");

write("~/MCcheckresults.txt","**************************");
flagerror=0;
counter=0;
for (q2=4,100,   
		phiq2=eulerphi(q2);
		div=divisors(q2);
		factors=factor(q2);
		nfactors=length(factors[,1]);
		ndiv=length(div);
		\\if (div[2]==2,i0=1,i0=0);
		i0=0;
		for (i=2+i0, ndiv-1,
				q1=div[i]; 
				phiq1=eulerphi(q1);
				n=q2\q1;
				for (a = 1, q1, 
						if (gcd(a,q1)==1, counter=counter+1;
								write("~/MCcheckresults.txt","The constant C(",q1,",",a,")", " has to be equal to the product of: ");
								P=1;
								for(j=0,n-1,
									if (gcd(a+j*q1,q2)==1,  
											for (i=1, phiq2,
													if (MC[position[q2]+i-1,2]==(a+j*q1)%q2,  
														P=P*MC[position[q2]+i-1,3];
														write("~/MCcheckresults.txt","C(",q2,",",(a+j*q1)%q2,") = ",MC[position[q2]+i-1,3]);
														);
													);
											);
									);
								for (l=1, nfactors,
									p=factors[l,1];
									if ((p-a)%q1==0,
										P=P*(1-1/p);
										write("~/MCcheckresults.txt","And of the factor: (1-",1/p,")");
										)
									);
							write("~/MCcheckresults.txt","Their product is = ",P);
							u=1;
							for(i=1, phiq1,
							   if (MC[position[q1]+i-1,2]==a,
							      write("~/MCcheckresults.txt","The constant C(",q1,",",a,")"," is = ",  MC[position[q1]+i-1,3]);
							      diff = P-MC[position[q1]+i-1,3];
							      write("~/MCcheckresults.txt","The difference (constant-product) is = ", diff);
								  write("~/MCcheckresults.txt","**************************");
								  if (abs(diff)>10^(-digitprecision), flagerror=1; write("~/MCcheckresults.txt","**************************");
								  				write("~/MCcheckresults.txt","ERROR: less than ",digitprecision," digits are correct");
								  				write("~/MCcheckresults.txt","**************************"));
								);
							);
						);
				);	
		);
	);
write("~/MCcheckresults.txt","**************************************************");
if (flagerror==0, print("We built and checked n. ", counter, " possible relationships between the given constants");
				print("All the results are correct with a precision of at least ",digitprecision," decimal digits");
				write("~/MCcheckresults.txt","We built and checked n. ", counter, " possible relationships between the given constants");
				write("~/MCcheckresults.txt","All the results are correct with a precision of at least ",digitprecision," decimal digits"),
				print("ERROR: At least one result has not the prescribed precision of more than ",digitprecision," decimal digits");
				write("~/MCcheckresults.txt","ERROR: At least one result has not the prescribed precision of more than ", digitprecision," decimal digits");
				);
print("******************* END CHECK ********************");
write("~/MCcheckresults.txt","******************* END CHECK ********************");
return();
}


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

languasco1:~ languasc$ gp
                  GP/PARI CALCULATOR Version 2.3.2 (released)
           i386 running darwin (ix86/GMP-4.2.2 kernel) 32-bit version
       compiled: Oct  8 2007, gcc-4.0.1 (Apple Computer, Inc. build 5367)
                (readline v5.0 enabled, extended help available)

                     Copyright (C) 2000-2006 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 ?12 for how to get moral (and possibly technical) support.

parisize = 4000000, primelimit = 500000
? \r /Users/languasc/ChecksMC.gp 
? checksMC
****** A. LANGUASCO & A. ZACCAGNINI *******
****** VERIFICATION OF EQ. (25) FOR 4 <= q <= 100 *******
****** REMARK: C(2,1) = 2 exp(-gamma); gamma=Euler constant *******
******************* START CHECK ********************
We built and checked n. 1907 possible relationships between the given constants
All the results are correct with a precision of at least 100 decimal digits
******************* END CHECK ********************
? ##
  ***   last result computed in 2,428 ms.
  
**********************/