/****************   	A. LANGUASCO & A. ZACCAGNINI 	*********************** 
************* COMPUTATION OF THE MEISSEL-MERTENS CONSTANTS MOD q; 3 <= q<= 100 ******* 
************* WITH A PRECISION OF AT LEAST 100 DECIMAL DIGITS ***************** 
************* COMPARISON BETWEEN B(q,a) directly computed and  ***************** 
*************   B(q,a)  computed using M(q,a) and C(q,a)       ***************** 
*******************************************************************************/


{controlB(limit) = local(MBdir, MB, E, i, bound);

default(realprecision,120);   
\\default(realprecision,40);   

\\ input of the constants
MB=read("~/MBmatrix.gp");

MBdir=read("~/Bmatrixdirect.gp");

bound=0;
for(i=2,limit, bound=bound+eulerphi(i));

print("************");
for(i=1, bound, 
                 E = MB[i,3]-MBdir[i,3];
                 print("Via C and M:");
                 print("q = ", MB[i,1], " a = ", MB[i,2], " B(q,a) = ", MB[i,3], " digits ", MB[i,4]);
                 print("Direct:");
                 print("q = ", MB[i,1], " a = ", MB[i,2], " B(q,a) = ", MBdir[i,3], " digits ", MBdir[i,4]);
                 print("q = ", MB[i,1], " a = ", MB[i,2], " Diff = ", E);
                 print("************");
   );
   
}