/*
    Copyright (C) 2021 Alessandro Languasco
    
    Computation of the following functions using the algorithms in the paper
    "Efficient computation some special functions"
    
    1) loggamma(x), x>0;
    2) psi(x), x>0;
    3) Bateman G function, x>0;  
    
    4) zeta Hurwitz (s,x), x>0, s>1; 
    	(in another file named hurwitz_forfft.gp there's version with x in (0,1) for FFT applications);
    4*) ONLY in another file: zetaprime Hurwitz (s,x),x in (0,1), s>1;
	    (ONLY in another file named hurwitz_forfft.gp there's version with x in (0,1) for FFT applications);
    
    5) Dirichlet beta (s), s>1;
    6) Dirichlet betaprime (s), s>1;
    7) Dirichlet beta values (s), s>1  (it computes betaprime/beta(s), beta(s) and betaprime(s) at the same time);
    
    8) Catalan constant computation.
    
    REMARKS:
    1) the values of the Riemann zeta function, of zetaprime, of Euler's Beta
    function, and of the differences of the digamma function are precomputed at
    runtime each time the main functions are called; another important
    gain in the running time would follow from having such values stored
    once for all (for example, they could be computed and stored when
	the first time this package is called)
	
	2) the variable LOOP (as default initialized as 1) is used to test
	the performances by repeating LOOP times the main computations each
	function performs
	
	3) Some of these functions requires a precision parameter that could be
	32, 64, 80, 128 bits and uses the optimal m-value pre-determined in 
	Tables 1 and 2 of the paper
	
	4) toward the bottom of this file you will also find some examples about
	how to use these functions and a comparison of their running times with 
	the analogue ones obtained with the internal PARI/GP functions.
     

*/

/****************   	A. LANGUASCO  		***************************/


/********************************************
 Here we compute the values of loggamma(x); x > 0
 we split in three parts: 
 1) the tail: x > myshift (with Euler-Maclaurin; myshift is the horizontal shift mentioned in the paper) 
 2) intermediate: 1<= x <= myshift (trivially using g_Gamma of the paper)
 3) {x} in (0,1) 		

 IMPORTANT: to have a correct representation of
 the result, the order of magnitude of the result
 has to be considered
 		
******************************************* */
  

\\  Euler-Maclaurin formula error function for loggamma
{gGamma(x,y,m) = local(res);   
	if (y>0 && y<3, res = 0 ,
  		res = x^(2*m+1)/(2*Pi^2) * factorial(2*m)/(2*Pi)^(2*m) * (1+3/4^(m+1))* (1-1/(1+x*(y-1))^(2*m+1))
  		);
   return(res);
}
  

\\ Looking for the optimal m (Euler-MacLaurin for loggamma) 
{opterrgamma(x)= local(m_opt, y, maxprec, toterr, toterr1, logten, logtwo, m);

\\choosing precision
default(realprecision, 100);  
print("internal precision: 100 decimal digits");

y=10^(10); 

toterr = gGamma(x,y,1);   
toterr1 = gGamma(x,y,2);  
m=2; 
\\print(toterr);
\\print(toterr1);
while( toterr1 < toterr,  
			toterr = toterr1;  
			m= m+1;  
			toterr1 = gGamma(x,y,m); 
			\\print(toterr);
			\\print(toterr1);
		);  
m_opt=m-1; \\ optimal m for this precision

print("optimal tilde-m = ", m_opt);
default(format,e); 
default(realprecision,10); 
print("minimal error = ", toterr);   

logten = log(10);
logtwo = log(2);
maxprec = floor(abs(log(toterr)/logten))-1;
print("total precision: at most ", maxprec, " decimal digits");
logtwo = log(2);
maxprec = floor(abs(log(toterr)/logtwo))-1;
print("total precision: at most ", maxprec, " binary digits");
print("******");   
return();
}

 
{loggamma(x, prec, LOOP=1)=  local(S, oneminusv, minusv, u , coeff, r1, m, datavector_zeta, abslogv, m_opt, maxindex,
abslogoneminusv, elaptime, valueatonehalf, starttime, start, correction, datavector_bern, quot, ok,
starttimeprecomp, elaptimeprecomp, y, y1, v, v1, myshift, magnitude, decprec, w, starttimehorshift, 
startv1, stepv1, fattorev1, startquot, fattorequot, stepquot, logten, logtwo, quad, cube, fourPifourth,
factstart, twoPi, twoPisquared, rec, stepfirstpart, stepsecondpart,firstpartstart, secondpartstart,
toterr, firstpart, secondpart, toterr1, elaptimehorshift, endEMtime, endshifttime,
endseriestime, startEMtime, startshifttime, startseriestime); 

default(format , f);  \\ print numbers with floating point notation (not E notation)

if (x<=0, error(" input has to be > 0 !!" ));


/*if (prec == 32, m = 32; myshift = 4; m_opt = 12,
	if (prec == 64, m = 64; myshift = 8; m_opt = 25,
		if (prec == 80, m = 80; myshift = 10; m_opt = 31,
			if (prec == 128, m = 128; myshift = 15; m_opt = 47, 
			 		error(" prec = 32, 64, 80, 128 (bits) !!!!")
			 	))));
*/	 	
	 	
if (prec <> 32 && prec <> 64 && prec <> 80 && prec <> 128, error(" prec = 32, 64, 80, 128 (bits) !!!!") );		

logtwo = log(2); 
logten = log(10);

decprec =  ceil(prec*logtwo/logten); \\ required decimal precision
\\ print("precision: at least ", decprec , " decimal digits");

y = floor(x); \\ y = integral part of x
v = frac(x); \\ v = fractional part of x
if ( v==0, x = floor(x) );

\\ maximal order of magnitude for  x->+infty and v->0+ of the function in base 10
magnitude = ceil ( log( abs(x*log(x)) ) / logten );
if (v <> 0, magnitude = max( magnitude , ceil ( abs(log(abs(log(v)))/logten) ) )); 

magnitude = ceil ( 1.1*magnitude );
\\ 1.1 just to be safer

\\print(magnitude);
decprec = magnitude + decprec + 16; \\ decimal precision  
\\ +16 just to be safer

default( realprecision, decprec );  
print("Internal precision: # decimal digits: ", decprec);

 
print("Evaluation of the horizontal shift");
starttimehorshift = getwalltime();

/* start computation horizontal shift */
\\m = prec ; 
m = ceil(decprec*logten/logtwo);	
	 	 
ok = 0; 
factstart = 2;
twoPi = 2*Pi;
twoPisquared = (twoPi)^2;
fourPifourth = (twoPisquared)^2;

for(j=1,10000, rec=1.0/j;
	  \\w=1;
	  quad =  rec^2;
	  cube = quad*rec;	
	  stepfirstpart = quad/twoPisquared;
	  stepsecondpart = 1/4;
	  
	  firstpartstart = 2* cube * factstart/fourPifourth;
	  secondpartstart = 1/16;
	  toterr = firstpartstart * (1+secondpartstart*5/3); 
	  
	  w=2; 	  
	  firstpart = firstpartstart* stepfirstpart * 3 * 4;
	  secondpart = secondpartstart * stepsecondpart;
	  toterr1=  firstpart *(1+secondpart*7.0/5);

	  while( toterr1 < toterr,  
			 toterr = toterr1;  
			 w = w+1;  
			 firstpart = firstpart * stepfirstpart * (2*w-1)*(2*w);
			 secondpart = secondpartstart * stepsecondpart;
			 toterr1 = firstpart *(1+secondpart*(2*w+3)/(2*w+1));
		  ); 
     if ( toterr < 10^(-decprec), ok = 1; m_opt= w - 1; myshift = j; break);
  	);
if (ok == 0, error(" Horizontal shift too large for the required precision; decrease the required accuracy"));  

/* end computation horizontal shift */

elaptimehorshift=getwalltime()-starttimehorshift; \\ elapsed time
print("horizontal shift = ", myshift);
print("Length of the EM-sum = ", m_opt);


\\ meaning of the parameters
\\ m = length of the sum for the computations in (0,1)
\\ m_opt = length of the Euler-MacLaurin sum
\\ myshift; horizontal shift to have the error in E-MC sufficiently small; precomputed to have 32, 64, 80, 128, bits 

\\ initialisation of the coefficients of the sums
starttimeprecomp = getwalltime();

maxindex = 2*m+10;
datavector_zeta = vector(maxindex,i, zeta(i+1)/(i+1)) ;  \\initialisation of the needed coefficients
maxindex = ceil(maxindex/2);
datavector_bern = vector(maxindex,i,bernreal(2*i)/(2*i*(2*i-1)));  \\precomputation of the first (weighted) even-index Bernoulli numbers
 
elaptimeprecomp = getwalltime()-starttimeprecomp; \\ elapsed time

starttime = getwalltime();
endEMtime = 0;
endshifttime = 0;
endseriestime = 0;

for (j=1, LOOP, 
	if (x == 1,     elaptime=getwalltime()-starttime;
			  		print("Raw time: ", elaptime, " millisec");
			  		return(0));  \\ special value  \log\Gamma(1) = 0

 
	S=0;
	
	startEMtime = getwalltime();
	\\ tail computation with Euler-MacLaurin formula; you have the tail if x > myshift
	if (x >= myshift, 
			v1 = v + myshift - 1; \\ change the real parameter based on the fractional part so that it is > myshift
			y1 = y - myshift + 1; \\ re-scaling the integral part to compensate
			\\ tail computation with Euler-MacLaurin formula (with myshift)  
			
			quot = 1/(1+(y1-1)/v1);
			
			\\ repeated products algorithm
			startv1 = 1.0/v1;
			startquot = quot;
			stepv1 = startv1 * startv1;
			stepquot = startquot * startquot;
		
			\\res = 0;
			fattorev1 = startv1;
			fattorequot = startquot;
			
			S += datavector_bern[1] * fattorev1 * ( - 1 + fattorequot );
			for (n=2,m_opt, 
				fattorev1 = fattorev1 * stepv1;
				fattorequot = fattorequot * stepquot;
				S += datavector_bern[n] * fattorev1 * ( - 1 + fattorequot );
				);
			
			S += - y1 * log(startv1) - (y1 -1) - (v1+y1-1/2) * log(quot);
		
		\\ the repeated products algorithm computed the following quantity
		\\	S= - y1 * log(1/v1) - (y1 -1) - (v1+y1-1/2) * log(quot) + sum(n=1, m_opt, datavector_bern[n]*v1^(-2*n+1)*(-1 + quot^(2*n-1)) , 0.);
	);
	endEMtime+= getwalltime()-startEMtime;
 

	\\compensate with the first summands (at most #=myshift summands)
	startshifttime = getwalltime();
	
	if(v<>0,  \\S += sum(k=0, min(myshift-2,y-1), log(v+k), 0.), \\ depending on the size of myshift this might be faster
		  S += log( prod(k=0, min(myshift-2,y-1), v+k,1.) ),
		  \\S += sum(k=2, min(myshift-2,y-1), log(k), 0.);  \\ depending on the size of myshift this might be faster
		  S += log( prod(k=2, min(myshift-2,y-1), k ,1.) ); 
		  elaptime=getwalltime()-starttime;
		  print("Raw time: ", elaptime, " millisec");
		  return(S);
		  ); 
    endshifttime += getwalltime()-startshifttime;

	\\ computing lngamma(v), with v={x} in (0,1);
	  
	startseriestime = getwalltime();
	coeff = (m+1)*logtwo; \\ needed for the length of sums
	valueatonehalf = 1/2*log(Pi);  \\ special value \log\Gamma(1/2) = 0.5 log (pi)

	if (v == 0.5,  S += valueatonehalf,  \\ v=1/2
	      \\ v <>1/2 -- initialisation common parameters
        minusv = -v;
        oneminusv  = 1 - v;
        abslogv = abs(log(v));	
        abslogoneminusv = abs(log(oneminusv)); 
        if (v > 0.5,  
		                 \\ 1 > v > 1/2 
 			r1 = ceil( ( coeff + abslogv) / abslogoneminusv )-1; \\number of summands 
 			if (r1<3, r1=3);  
 			start = oneminusv;
 			correction =Euler*oneminusv
 			,
		                 \\ 0 < v < 1/2 
 			r1 = ceil( ( coeff + abslogoneminusv ) / abslogv )-1; \\number of summands 
 			if (r1<3, r1=3);    
            start = minusv;
            correction = abslogv + Euler * minusv;
        );   
 
        u = sqr(start); \\ starting point for computing using a repeated product strategy
        for (k=2, r1,  
 			  S+= u * datavector_zeta[k-1]; 
 			  u*= start);   	
        S = correction + S; 	    \\ correction values
	);
    endseriestime += getwalltime()-startseriestime;
); \\end of the LOOP iteration

elaptime=getwalltime()-starttime; \\ elapsed time

print("Raw horizontal shift determination time: ", elaptimehorshift, " millisec"); 
print("Raw precomputation time: ", elaptimeprecomp, " millisec");
print("Raw EM comp time: ", endEMtime, " millisec");
print("Raw shift comp time: ", endshifttime, " millisec");
print("Raw series comp time: ", endseriestime, " millisec");
print("Raw total time (no precomp, no horshift det): ", elaptime, " millisec");
return(S);
}



/********************************************
 Here we compute the values of digamma(x); x > 0
 we split in three parts: 
 1) the tail: x > myshift (with Euler-Maclaurin; myshift is the horizontal shift mentioned in the paper)  
 2) intermediate: 1<= x <= myshift (trivially)
 3) {x} in (0,1) 	
 
 IMPORTANT: to have a correct representation of
 the result, the order of magnitude of the result
 has to be considered
 			
******************************************* */
 
 \\  Euler-Maclaurin formula error function for digamma
{gdigamma(x,y,m) = local(res);   
	if (y>0 && y<3, res = 0 ,
  		res = x^(2*m+2)/(2*Pi^2) * factorial(2*m+1)/(2*Pi)^(2*m) * (1+3/4^(m+1))* (1-1/(1+x*(y-1))^(2*m+1))
  		);
   return(res);
}

\\ Looking for the optimal m (Euler-MacLaurin for  digamma) 
{opterrdigamma(x)= local(m_opt, y, maxprec, toterr, toterr1, m, logten, logtwo);

\\choosing precision
default(realprecision, 100);  
print("internal precision: 100 decimal digits");

y=10^(10); 

toterr = gdigamma(x,y,1);   
toterr1 = gdigamma(x,y,2);  
m=2; 
\\print(toterr);
\\print(toterr1);
while( toterr1 < toterr,  
			toterr = toterr1;  
			m= m+1;  
			toterr1 = gdigamma(x,y,m); 
			\\print(toterr);
			\\print(toterr1);
		);  
m_opt=m-1; \\ optimal m for this precision

print("optimal tilde-m = ", m_opt);
print("minimal error = ", toterr);   
logten = log(10);
logtwo = log(2);

maxprec = floor(abs(log(toterr)/logten))-1;
print("total precision: at most ", maxprec, " decimal digits");
maxprec = floor(abs(log(toterr)/logtwo))-1;
print("total precision: at most ", maxprec, " binary digits");
print("******");   
return();
}


{digamma(x, prec, LOOP=1)=  local(S, oneminusv, minusv, u , coeff, r1, m, datavector_zeta, abslogv, S1, 
m_opt, maxindex, starttimeprecomp, elaptimeprecomp, abslogoneminusv, elaptime, decprec, 
valueatonehalf, starttime, start, correction, quot, magnitude, datavector_bern, y, y1, v, 
v1, myshift, startv1, stepv1, fattorev1, startquot, fattorequot, stepquot, logten, logtwo,
factstart, twoPi, twoPisquared, rec, stepfirstpart, stepsecondpart,firstpartstart, secondpartstart,
toterr, firstpart, secondpart, toterr1, elaptimehorshift, ok, starttimehorshift, fourPifourth,
fourth, quad, endEMtime, endshifttime,
endseriestime, startEMtime, startshifttime, startseriestime); 

default(format , f);  \\ print numbers with floating point notation (not E notation)

if (x<=0, error(" input has to be > 0 !!" ));

/*
if (prec == 32, m = 32; myshift = 4; m_opt = 12,
	if (prec == 64, m = 64; myshift = 8; m_opt = 24,
		if (prec == 80, m = 80; myshift = 10; m_opt = 31,
			if (prec == 128, m = 128; myshift = 15; m_opt = 46, 
			 		error(" prec = 32, 64, 80, 128 (bits) !!!!")
			 	))));		 
*/
if (prec <> 32 && prec <> 64 && prec <> 80 && prec <> 128, error(" prec = 32, 64, 80, 128 (bits) !!!!") );		


logtwo = log(2); 
logten = log(10);

decprec =  ceil(prec*logtwo/logten); \\ required decimal precision
\\ print("precision: at least ", decprec , " decimal digits");

y = floor(x); \\ y = integral part of x
v = frac(x); \\ v = fractional part of x
if ( v==0, x = floor(x) );

\\ maximal order of magnitude for  x->+infty and v->0+ of the function in base 10
magnitude =  ceil ( log( abs( log(x)) ) / logten );
if (v <> 0, magnitude = max( magnitude , ceil ( abs(log(v)/logten) ))); 

magnitude = ceil ( 1.1*magnitude );
\\ 1.1 just to be safer

decprec = magnitude + decprec + 16; \\ decimal precision  
\\ +16 just to be safer

default( realprecision, decprec );  
print("Internal precision: # decimal digits: ", decprec);


print("Evaluation of the horizontal shift");

starttimehorshift = getwalltime();

/* start computation horizontal shift */
\\m = prec ; 
m = ceil(decprec*logten/logtwo);	
	 	 
ok = 0; 
factstart = 6;
twoPi = 2*Pi;
twoPisquared = (twoPi)^2;
fourPifourth = (twoPisquared)^2;

for(j=1,10000, rec=1.0/j;
	  \\w=1;
	  quad = rec^2;
	  fourth = quad*quad;	
	  stepfirstpart = quad/twoPisquared;
	  stepsecondpart = 1/4;
	  
	  firstpartstart = 2* fourth * factstart/fourPifourth;
	  secondpartstart = 1/16;
	  toterr = firstpartstart * (1+secondpartstart*5/3); 
	  
	  w=2; 	  
	  firstpart = firstpartstart* stepfirstpart * 4 * 5;
	  secondpart = secondpartstart * stepsecondpart;
	  toterr1=  firstpart *(1+secondpart*7/5);

	  while( toterr1 < toterr,  
			 toterr = toterr1;  
			 w = w+1;  
			 firstpart = firstpart * stepfirstpart * (2*w)*(2*w+1);
			 secondpart = secondpartstart * stepsecondpart;
			 toterr1 = firstpart *(1+secondpart*(2*w+3)/(2*w+1));
		  ); 
     if ( toterr < 10^(-decprec), ok = 1; m_opt= w - 1; myshift = j; break);
  );
if (ok == 0, error(" Horizontal shift too large for the required precision; decrease the required accuracy"));  

/* end computation horizontal shift */

elaptimehorshift=getwalltime()-starttimehorshift; \\ elapsed time

print("horizontal shift = ", myshift);
print("Length of the EM-sum = ", m_opt);


\\ meaning of the parameters
\\ m = length of the sum for the computations in (0,1)
\\ m_opt = length of the Euler-MacLaurin sum
\\ myshift; myshift to have the error in E-MC sufficiently small; precomputed to have 32, 64, 80, 128, bits 


\\ initialisation of the coefficients of the sums
starttimeprecomp = getwalltime();

maxindex = 2*m+10;
datavector_zeta = vector(maxindex,i, zeta(i+1)) ;  \\initialisation of the needed coefficients
maxindex = ceil(maxindex/2);
datavector_bern = vector(maxindex,i,bernreal(2*i)/(2*i));  \\precomputation of the first (weighted) even-index Bernoulli numbers
 
elaptimeprecomp=getwalltime()-starttimeprecomp; \\ elapsed time

starttime = getwalltime();
endEMtime = 0;
endshifttime = 0;
endseriestime = 0;

for (j=1, LOOP, 

	if (x == 1,     elaptime=getwalltime()-starttime;
		  		print("Raw time: ", elaptime, " millisec");
		  		return(-Euler));  \\ special value  \digamma(1) = -Euler

	S=0;
	startEMtime = getwalltime();
	\\ tail computation with Euler-MacLaurin formula; you have the tail if x > myshift
	if (x >= myshift, 
			v1 = v + myshift - 1 ; \\ change the real parameter based on the fractional part so that it is > myshift
			y1 = y - myshift + 1;  \\ re-scaling the integral part to compensate
			\\ tail computation with Euler-MacLaurin formula (with myshift)  
			
			quot = 1/(1+(y1-1)/v1);
			
			\\ repeated products algorithm
			startv1 = v1^(-2);
			startquot = quot*quot;
			stepv1 = startv1;
			stepquot = startquot;
		
			fattorev1 = startv1;
			fattorequot = startquot;
			S = datavector_bern[1] * fattorev1 * ( 1 - fattorequot );		
			
			for (n=2,m_opt, 
				fattorev1 = fattorev1 * stepv1;
				fattorequot = fattorequot * stepquot;
				S += datavector_bern[n] * fattorev1 * ( 1 - fattorequot );
				);
			
			S += -log(quot) + (1+quot)/(2*v1) ;

		   \\ the repeated products algorithm computed the following quantity			
		   \\ S = -log(quot)  + 0.5/v1*(1+quot)+ sum(n=1, m_opt, datavector_bern[n]*v1^(-2*n)*(1-quot^(2*n)) , 0.)
			
	);
	endEMtime+= getwalltime()-startEMtime;
 
	\\compensate with the first summands (at most #=myshift summands)
  	startshifttime = getwalltime();
  	
	if(v<>0,  S += sum(k=0, min(myshift-2,y-1), 1/(v+k),0.),
		  S += -Euler + sum(k=1, min(myshift-2,y-1), 1/k,0.); 
		  elaptime=getwalltime()-starttime;
		  print("Raw time: ", elaptime, " millisec");
		  return(S);
		  ); \\ sum of the first terms
    endshifttime += getwalltime()-startshifttime;
    
	\\ computing digamma(v), with v={x} in (0,1); 

	startseriestime = getwalltime();
	coeff = (m+1)*logtwo; \\ needed for the length of sums
	valueatonehalf = -2*logtwo-Euler;  \\ special value digamma(1/2) = -2 log(2) - gamma

	if (v == 0.5,  S += valueatonehalf,  \\ v=1/2
	      \\ v <>1/2 -- initialisation common parameters
	    S1=0;
        minusv = -v;
        oneminusv  = 1 - v;
        abslogv = abs(log(v));	
        abslogoneminusv = abs(log(oneminusv)); 
        if (v > 0.5,  
		                 \\ 1 > v > 1/2 
 			r1 = ceil( ( coeff + abslogv + 0.2 ) / abslogoneminusv )- 1; \\number of summands  
 			if (r1<3, r1=3);     
 			start = oneminusv;
 			correction = - Euler
 			,
		                 \\ 0 < v < 1/2 
 			r1 = ceil( ( coeff + abslogoneminusv + 0.2) / abslogv ) - 1; \\number of summands   
 			if (r1<2, r1=3);    
            start = minusv;
            correction = - Euler + 1/minusv;
        );   
 
        u = start; \\ starting point for computing using a repeated product strategy
        for (k=2, r1,  
 			  S1+= u * datavector_zeta[k-1]; 
 			  u*= start);    
        S1 = correction - S1; 	    \\ correction values
	);

	S = S + S1;
    endseriestime += getwalltime()-startseriestime;
); \\end of the LOOP iteration	

elaptime=getwalltime()-starttime; \\ elapsed time

print("Raw horizontal shift determination time: ", elaptimehorshift, " millisec"); 
print("Raw precomputation time: ", elaptimeprecomp, " millisec");
print("Raw EM comp time: ", endEMtime, " millisec");
print("Raw shift comp time: ", endshifttime, " millisec");
print("Raw series comp time: ", endseriestime, " millisec");
print("Raw total time (no precomp, no horshift det): ", elaptime, " millisec");
return(S);
}

/********************************************
 Here we compute the values of batemanG(x); x > 0
 we split in three parts: see section 8.5 of the paper for the cases
 we do about {x/2}
 1) the tail: x > myshift (with Euler-Maclaurin; myshift is the horizontal shift mentioned in the paper) 
 2) intermediate: 1<= x <= myshift (trivially)
 3) {x} in (0,1) 				

 IMPORTANT: to have a correct representation of
 the result, the order of magnitude of the result
 has to be considered
 
******************************************* */

{verifbatemanG(x) = psi((x+1)/2) - psi(x/2);
}

\\\\  BATEMAN G  

{batemanG(x, prec, LOOP=1)=  local(S, oneminusv2, minusv2, u , coeff, r1, m, datavector_zeta, abslogv2, 
S1, m_opt, xhalf, w1, fracx, oneplusv2, magnitude, decprec, twoplusv2,
maxindex, starttimeprecomp, elaptimeprecomp, logtwo, logten, twicelogtwo, multiplier,
abslogoneminusv2, elaptime, valueatonehalf, starttime, start, correction, valueatone,
datavector_powers, datavector_bern, y, y1, v, v1, v2, myshift, coeffvector,
factstart, twoPi, twoPisquared, rec, stepfirstpart, stepsecondpart,firstpartstart, secondpartstart,
toterr, firstpart, secondpart, toterr1, elaptimehorshift, ok, starttimehorshift, fourPifourth,
fourth, quad, w, endEMtime, endshifttime, endseriestime, startEMtime, startshifttime, startseriestime,
quotv, quotw, startv1, startw1, startquotv, startquotw, stepv1, stepw1, stepquotv, stepquotw,
fattorev1, fattorequotv, fattorew1, fattorequotw); 

default(format , f);  \\ print numbers with floating point notation (not E notation)

if (x<=0, error(" input has to be > 0 !!" ));

/*  done with the parameters for psi(x) in the EM formula */

/*
if (prec == 32, m = 32; myshift = 4; m_opt = 12,
	if (prec == 64, m = 64; myshift = 8; m_opt = 24,
		if (prec == 80, m = 80; myshift = 10; m_opt = 31,
			if (prec == 128, m = 128; myshift = 15; m_opt = 46, 
			 		error(" prec = 32, 64, 80, 128 (bits) !!!!")
			 	))));	
*/			 	
			 	
if (prec <> 32 && prec <> 64 && prec <> 80 && prec <> 128, error(" prec = 32, 64, 80, 128 (bits) !!!!") );		
			 			 	 			 	

logtwo = log(2);
logten = log(10);

twicelogtwo = 2 * logtwo;
valueatone = twicelogtwo;  \\ special value G(1) = 2*log(2)
valueatonehalf = Pi;  \\ special value G(1/2) = Pi


decprec =  ceil(prec*logtwo/logten); \\ required decimal precision
\\ print("precision: at least ", decprec , " decimal digits");

y = floor(x); \\ y = integral part of x
v = frac(x); \\ v = fractional part of x
if ( v==0, x = floor(x) );

\\ maximal order of magnitude for  x->+infty and v->0+ of the function in base 10 (as for digamma)
magnitude = ceil ( log( abs(log(x))));
if (v <> 0, magnitude = max( magnitude , ceil ( abs(log(v) / logten )))); 
 
magnitude = ceil ( 1.1*magnitude );
\\ 1.1 just to be safer

decprec = magnitude + decprec + 16; \\ decimal precision  
\\ +16 just to be safer

default( realprecision, decprec );  
print("Internal precision: # decimal digits: ", decprec);



print("Evaluation of the horizontal shift");

starttimehorshift = getwalltime();

/* start computation horizontal shift */
\\m = prec ; 
m = ceil(decprec*logten/logtwo);	
	 	 
ok = 0; 
factstart = 6;
twoPi = 2*Pi;
twoPisquared = (twoPi)^2;
fourPifourth = (twoPisquared)^2;

for(j=1,10000, rec=1.0/j;
	  \\w=1;
	  quad = rec^2;
	  fourth = quad*quad;	
	  stepfirstpart = quad/twoPisquared;
	  stepsecondpart = 1/4;
	  
	  firstpartstart = 2* fourth * factstart/fourPifourth;
	  secondpartstart = 1/16;
	  toterr = firstpartstart * (1+secondpartstart*5/3); 
	  toterr = 2*toterr; \\ change for the G-function
	  
	  w=2; 	  
	  firstpart = firstpartstart* stepfirstpart * 4 * 5;
	  secondpart = secondpartstart * stepsecondpart;
	  toterr1=  firstpart *(1+secondpart*7/5);
	  toterr1 = 2*toterr1; \\ change for the G-function
	  
	  while( toterr1 < toterr,  
			 toterr = toterr1;  
			 w = w+1;  
			 firstpart = firstpart * stepfirstpart * (2*w)*(2*w+1);
			 secondpart = secondpartstart * stepsecondpart;
			 toterr1 = firstpart *(1+secondpart*(2*w+3)/(2*w+1));
			 toterr1 = 2*toterr1; \\ change for the G-function
		  ); 
     if ( toterr < 10^(-decprec), ok = 1; m_opt= w - 1; myshift = j; break);
  );
if (ok == 0, error(" Horizontal shift too large for the required precision; decrease the required accuracy"));  

/* end computation horizontal shift */

elaptimehorshift=getwalltime()-starttimehorshift; \\ elapsed time
print("horizontal shift = ", myshift);
print("Length of the EM-sum = ", m_opt);

\\ meaning of the parameters
\\ m = length of the sum for the computations in (0,1)
\\ m_opt = length of the Euler-MacLaurin sum
\\ myshift; myshift to have the error in E-MC sufficiently small; precomputed to have 32, 64, 80, 128, bits 



\\ initialisation of the coefficients of the sums
starttimeprecomp = getwalltime();

maxindex = 2*m+10;
datavector_zeta = vector(maxindex,i, zeta(i+1)-1) ;  \\initialisation of the needed coefficients
datavector_powers = vector(maxindex,i, (1 - 2^(-i))*1.0 );  
coeffvector = vector(maxindex,i, 0);  
for ( i = 1, maxindex, coeffvector[i] =  datavector_zeta[i] * datavector_powers[i]);


\\print(datavector_powers);
maxindex = ceil(maxindex/2);
datavector_bern = vector(maxindex,i,bernreal(2*i)/(2*i));  \\precomputation of the first (weighted) even-index Bernoulli numbers
 
elaptimeprecomp=getwalltime()-starttimeprecomp; \\ elapsed time

starttime = getwalltime();
endEMtime = 0;
endshifttime = 0;
endseriestime = 0;

for(j=1, LOOP,

	if (x == 1,     elaptime=getwalltime()-starttime;
		  		print("Raw time: ", elaptime, " millisec");
		  		return( valueatone ));  \\ special value  G(1) = 2*log(2)
		  		
	if (x == 0.5,     elaptime=getwalltime()-starttime;
		  		print("Raw time: ", elaptime, " millisec");
		  		return( valueatonehalf ));  \\ special value  G(1) = 2*log(2)

	xhalf=x/2;
	y=floor(xhalf); \\ y = integral part of x/2
	v = frac(xhalf); \\ v = fractional part of x/2
	fracx = frac(x); \\ fracx = fractional part of x


	S=0;
	
	startEMtime = getwalltime();
	\\ tail computation with Euler-MacLaurin formula; you have the tail if x > myshift
	if (xhalf >= myshift,  
			if(fracx == 0 && x <> 2*y, v = v-0.5); \\ x odd integer
			
			v1 = v + myshift - 1 ; \\ change the real parameter based on the fractional part so that it is > myshift
			w1 = v + myshift - 1/2;
			y1 = y - myshift + 1;  \\ re-scaling the integral part to compensate
			\\ tail computation with Euler-MacLaurin formula (with myshift)  

		  	quotv = 1/(v1+y1-1);
		  	quotw = 1/(w1+y1-1);
			
			\\ repeated products algorithm
			startv1 = v1^(-2);
			startw1 = w1^(-2);
			startquotv = quotv*quotv;
			startquotw = quotw*quotw;
			stepv1 = startv1;
			stepw1 = startw1;
			stepquotv = startquotv;
			stepquotw = startquotw;
		
			fattorev1 = startv1;
			fattorequotv = startquotv;
			fattorew1 = startw1;
			fattorequotw = startquotw;		
				
			S = datavector_bern[1] * (fattorew1 - fattorev1 - fattorequotw + fattorequotv );		
			
			for (n=2,m_opt, 
				fattorev1 = fattorev1 * stepv1;
				fattorew1 = fattorew1 * stepw1;				
				fattorequotv = fattorequotv * stepquotv;
				fattorequotw = fattorequotw * stepquotw;
				S += datavector_bern[n] * (fattorew1 - fattorev1 - fattorequotw + fattorequotv );
				);
		  	
		  	\\S += log(v1/w1) + log( (w1+y1-1) / (v1+y1-1) ) + 1/2*(1/w1 -1/v1) + 1/2*( 1/(w1+y1-1) - 1/(v1+y1-1));
		  	S += log( v1 / w1 * quotv / quotw ) + 1/2*(1/w1 -1/v1) + 1/2*( quotw - quotv);
		  	
		  	\\ the repeated products algorithm computed the following quantity			
		  	\\S = log(v1/w1) + log( (w1+(y1-1)) / (v1+(y1-1)) )
			\\	+ 0.5*(1/w1 -1/v1) + 0.5*( 1/(w1+y1-1) - 1/(v1+y1-1))
		  	\\	+ sum(n=1, m_opt, datavector_bern[n] * ( w1^(-2*n) - v1^(-2*n) - 1/(w1+y1-1)^(2*n) + 1/(v1+y1-1)^(2*n) ) , 0.); 
 
			
			if(fracx == 0 && x <> 2*y, S = - S); \\ x odd integer
 
	);
	endEMtime+= getwalltime()-startEMtime;
 
	\\compensate with the first summands (at most #=myshift summands)
	startshifttime = getwalltime();
	
	v = frac(xhalf); \\ v = fractional part of x/2 (reset)
  
	if(fracx <> 0,   S = S - 1/2 * sum(k=0, min(myshift-2,y-1), 1/((v+k)*(v+k+1/2)),0. ), \\ 1/(v+k+0.5) - 1/(v+k) = - 0.5/((v+k)*(v+k+0.5))
		if( x == 2*y, 
			\\print("even");
		  		S =  S  + 2 - valueatone  - 1/2 * sum(k=1, min(myshift-2,y-1), 1/(k*(k+1/2)),0. ); \\ 1/(k+0.5)-1/k = - 0.5/(k*(k+0.5))
		  		elaptime=getwalltime()-starttime;
		  		print("Raw time: ", elaptime, " millisec");
		  		return(S),  \\ psi(3/2)-psi(1) = 2 - 2*log(2)
          	\\print("odd"); 
          		S =  S  + 1/y - 2 + valueatone  + 1/2 * sum(k=1, min(myshift-2,y-1), 1/(k*(k+1/2)),0. ); 
          		elaptime=getwalltime()-starttime;
		  		print("Raw time: ", elaptime, " millisec");
		  		return(S);
		  )); 
    endshifttime += getwalltime()-startshifttime;

	\\computing batemanG(v), with v={x/2} in (0,1); 

	startseriestime = getwalltime();
	coeff = (m+1)*logtwo; \\ needed for the length of sums
 
	S1=0; 
	\\print(S);
	if (v == 0.25 ,  S += valueatonehalf, \\print(S),  \\ {x}=1/2
	 
		if ( x < 1, v2 = x , \\compute G(v2) with the zeta series; READ THE PAPER !!
			if (v > 0.5, v2=2*v-1, v2=2*v); \\ compute G(v2) with the zeta series; READ THE PAPER !!
			); 
        
        minusv2 = -v2;
        oneplusv2 = 1 + v2;
        twoplusv2 = 2 + v2;        
        oneminusv2 = 1 - v2; 
        
        abslogv2 = abs(log(v2));	
        abslogoneminusv2 = abs(log(oneminusv2));
         
        if (v2 > 0.5,  
		                 \\ 1 > v2 > 1/2  
            
 			r1 = ceil( ( coeff + abslogv2 + 1.2 ) / abslogoneminusv2 ) - 1 ; \\number of summands   
  			if (r1<2, r1=2);  
 			start = oneminusv2;
 			correction = 2*(logtwo + oneminusv2/(v2*oneplusv2));
 			multiplier = 2;
 			,
		                 \\ 0 < v2 < 1/2 
 			r1 = ceil( ( coeff + abslogoneminusv2 + 1.2 ) / abslogv2 ) -1 ; \\number of summands   
  			if (r1<2, r1=2);  
            start = minusv2;
            correction = 2*(-logtwo + v2/(oneplusv2*twoplusv2) + 1/v2) ;
            multiplier = -2;
        );   
 
        u = start; \\ starting point for computing using a repeated product strategy
 
        for (k=1, r1,  
 			  \\S1 += u * datavector_zeta[k] * datavector_powers[k]; 
 			  S1 += u * coeffvector[k]; 
 			  u *= start;  
 			  );
         S1 = correction + multiplier * S1; 	    \\ correction values

		if ( x < 1, S = S + S1 , if (v > 0.5,  S = S - S1 + 1/(v-0.5), S = S + S1 )); 
	);
	endseriestime += getwalltime()-startseriestime;
); \\end of the LOOP iteration	

elaptime=getwalltime()-starttime; \\ elapsed time
 
print("Raw horizontal shift determination time: ", elaptimehorshift, " millisec"); 
print("Raw precomputation time: ", elaptimeprecomp, " millisec");
print("Raw EM comp time: ", endEMtime, " millisec");
print("Raw shift comp time: ", endshifttime, " millisec");
print("Raw series comp time: ", endseriestime, " millisec");
print("Raw total time (no precomp, no horshift det): ", elaptime, " millisec");

return(S);
}

 
  

/********************************************
 Here we compute the values of hurwitzzeta(s,x); s > 1, x > 0  
 we split in three parts: 
 1) the tail: x > myshift (with Euler-Maclaurin; myshift is the horizontal shift mentioned in the paper)  
 2) intermediate: 1<= x <= myshift (trivially)
 3) {x} in (0,1) 				
 
 REMARK: horizontal shift evaluation performed at runtime		
 
 IMPORTANT: to have a correct representation of
 the result, the order of magnitude of the result
 has to be considered
 ******************* */
 
 
{hurwitzzeta(s, x, prec, LOOP=1)=  local(S, oneminusv, minusv, u , coeff, r1, m, datavector_zeta, maxindex, abslogv, S1, start, 
starttimeprecomp, elaptimeprecomp, rec, ok, w, elaptimehorshift, starttimehorshift, toterr, toterr1, logten, logtwo,
m_opt, abslogoneminusv, elaptime, decprec, valueatonehalf, gammaratiosvector, factstart, twoPi, twoPisquared, y, v,
stepsecondpart, stepfirstpart, firstpart, secondpart, firstpartstart, secondpartstart, sminusone, minuss, quot,magnitude,
myshift, datavector_bern, coeffvector, starttime, v1, y1, ceils, ceilsplusone, ceilsplustwo, correction, EMcoeffvector,
oneplusv, fourPifourth, quad, cube, endEMtime, endshifttime, endseriestime, startEMtime, startshifttime, startseriestime,
quotsminusone, v1tominuss, v1tooneminuss, startv1, done, startquot, stepv1, stepquot, fattorev1, fattorequot, zetas); 

default(format , f);  \\ print numbers with floating point notation (not E notation)
 

if (x <= 0, error(" x input has to be > 0 !!" )); 
if (s <= 1, error(" s input has to be > 1 !!" ));
if (frac(s) == 0, s=floor(s));
			 	
if (prec <> 32 && prec <> 64 && prec <> 80 && prec <> 128, error(" prec = 32, 64, 80, 128 (bits) !!!!") );		

logten = log(10);
logtwo = log(2);

decprec =  ceil(prec*logtwo/logten); \\ required decimal precision
\\ print("precision: at least ", decprec , " decimal digits");
 

y = floor(x); \\ y = integral part of x
v = frac(x); \\ v = fractional part of x
if ( v==0, x = floor(x) );

\\ max order of magnitude as x->+infty and v->0+ of the function (base 10)

magnitude = ceil(s * abs(log(x)) / logten );
if (v <> 0, magnitude = max( magnitude , ceil ( s*abs(log(v))/logten ) ) ); 

magnitude = ceil ( 1.1*magnitude );

\\print(magnitude);
decprec = magnitude + decprec + 16; \\ decimal precision  
\\ +16 solo per "starci largo"; per avere sempre altre 16 cifre decimali dopo la virgola

default( realprecision, decprec );  
print("Internal precision: # decimal digits: ", decprec);

if (x == 1, return(zeta(s)));  \\ special value  zetahurwitz(s,1) =  zeta(s)

valueatonehalf =  (2^s-1)*zeta(s);  
if (x == 0.5, return(valueatonehalf));


print("Evaluation of the horizontal shift");
starttimehorshift = getwalltime();
 
/* start computation horizontal shift */
\\m = prec ; 
m = ceil(decprec*logten/logtwo);
	 	 
ok = 0; 
factstart=s*(s+1)*(s+2);
twoPi = 2*Pi;
twoPisquared = (twoPi)^2;
fourPifourth = (twoPisquared)^2;

for(j=1,10000, rec=1.0/j;
	  \\w=1;
	  quad =  rec^2;
	  cube = quad*rec;	
	  stepfirstpart = quad/twoPisquared;
	  stepsecondpart = 1/4;
	  
	  firstpartstart = 2*rec^s*cube * factstart/fourPifourth;
	  secondpartstart = 1/16;
	  toterr = firstpartstart * (1+secondpartstart*5/3); 
	  
	  w=2; 	  
	  firstpart = firstpartstart* stepfirstpart * (s+3)*(s+4);
	  secondpart = secondpartstart * stepsecondpart;
	  toterr1=  firstpart *(1+secondpart*7/5);

	  while( toterr1 < toterr,  
			 toterr = toterr1;  
			 w = w+1;  
			 firstpart = firstpart * stepfirstpart * (s+2*w-1)*(s+2*w);
			 secondpart = secondpartstart * stepsecondpart;
			 toterr1 = firstpart *(1+secondpart*(2*w+3)/(2*w+1));
		  ); 
     if ( toterr < 10^(-decprec), ok = 1; m_opt= w - 1; myshift = j; break);
  );
if (ok == 0, error(" Horizontal shift too large for the required precision; decrease the required accuracy"));  

/* end computation horizontal shift */

elaptimehorshift=getwalltime()-starttimehorshift; \\ elapsed time
print("horizontal shift = ", myshift);
print("Length of the EM-sum = ", m_opt);

\\ meaning of the parameters
\\ m = length of the sum for the computations in (0,1)
\\ m_opt = length of the Euler-MacLaurin sum
\\ myshift; myshift to have the error in E-MC sufficiently small; computed at runtime 

 
\\ initialisation of the coefficients of the sums
starttimeprecomp = getwalltime();


maxindex = ceil(3*(m+2*s)); 
 
zetas=zeta(s);
datavector_zeta = vector(maxindex,i, zeta(s+i-1)-1) ;  \\initialisation of the needed coefficients
 
gammaratiosvector = vector(maxindex,i, 0) ;  
gammaratiosvector[1] = s;
for ( i = 2, maxindex, gammaratiosvector[i] =  gammaratiosvector[i-1]*(s+i-1)/i);
coeffvector = vector(maxindex,i, 0);  
for ( i = 1, maxindex-1, coeffvector[i] =  gammaratiosvector[i]* datavector_zeta[i+1]);

maxindex = ceil(maxindex/2);
datavector_bern = vector(maxindex,i,bernreal(2*i)/(2*i));   

EMcoeffvector = vector(maxindex,i, 0);  
for ( i = 1, maxindex, EMcoeffvector[i] =  datavector_bern[i] * gammaratiosvector[2*i-1]);


elaptimeprecomp=getwalltime()-starttimeprecomp; \\ elapsed time

starttime = getwalltime(); 
endEMtime = 0;
endshifttime = 0;
endseriestime = 0;


for(j=1, LOOP,
	done = 0;
	y = floor(x); \\ y = integral part of x
	v = frac(x); \\ v = fractional part of x
	if ( v==0, x = floor(x) );


	S=0;

	sminusone = s-1;
	minuss = -s;

	startEMtime = getwalltime();
	\\ tail computation with Euler-MacLaurin formula; you have the tail if x > myshift
	if (x >= myshift, 
			v1 = v + myshift - 1 ; \\ change the real parameter based on the fractional part so that it is > myshift
			y1 = y - myshift + 1;  \\ re-scaling the integral part to compensate

			\\ tail computation with Euler-MacLaurin formula (with myshift) 
		  	quot = 1/(1+(y1-1)/v1);
		  	quotsminusone = quot^sminusone;
		  	v1tominuss = v1^(minuss);
		  	v1tooneminuss = v1 * v1tominuss;
		  	 	
		  	
		  	\\ repeated products algorithm
			startv1 = v1^(-2);
			startquot = quot^2;
			stepv1 = startv1;
			stepquot = startquot;

			\\res = 0;
			fattorev1 = startv1;
			fattorequot = startquot;
			
			\\S += datavector_bern[1] * gammaratiosvector[1] * fattorev1 * ( 1 - fattorequot*quotsminusone );
			S += EMcoeffvector[1] * fattorev1 * ( 1 - fattorequot*quotsminusone );
			
			for (n=2,m_opt, 
				fattorev1 = fattorev1 * stepv1;
				fattorequot = fattorequot * stepquot;
				\\S += datavector_bern[n] * gammaratiosvector[2*n-1] * fattorev1 * ( 1 - fattorequot*quotsminusone );
				S += EMcoeffvector[n] * fattorev1 * ( 1 - fattorequot*quotsminusone );
				);  	
		  	
		  	S = v1tooneminuss * S;
		  	S += v1tooneminuss/sminusone*(1 - quotsminusone) + v1tominuss * 1/2 *(1+quotsminusone*quot);
		  				
		  	\\ the repeated products algorithm computed the following quantity
		  	\\S = v1^(-sminusone)/sminusone*(1 - quot^sminusone) + v1^(minuss)*0.5*(1+quot^s)
			\\	+ v1^(-sminusone) * sum(n=1, m_opt, datavector_bern[n]*gammaratiosvector[2*n-1]/(2*n)*v1^(-2*n)*(1 - quot^(2*n+sminusone)) , 0.); 
  			
  			S = -S;	

	);
 	endEMtime+= getwalltime()-startEMtime;	 
 
	\\compensate with the first summands (at most #=myshift summands) 
	startshifttime = getwalltime();

	if(v<>0,  S += (-1)*sum(k=0, min(myshift-2,y-1), (v+k)^(minuss),0.), 
		  S += zetas - 1 + (-1)*sum(k=2, min(myshift-2,y-1), k^(minuss),0.);  
		  done = 1;
		  \\return(S);
		  );  
    endshifttime += getwalltime()-startshifttime;

	\\ computing ZETAHURWITZ(s,v), with v={x} in (0,1)  
	if ( done == 0, 		
		startseriestime = getwalltime();
		coeff = (m+1)*logtwo; \\ needed for the length of sums 
		ceils = ceil(s);
		ceilsplusone = ceils+1;
		ceilsplustwo = ceils+2; 

		if (v == 0.5,  S += valueatonehalf,  \\ v=1/2 
	      	\\ v <>1/2 -- initialisation common parameters
	    	S1=0;
        	minusv = -v;
        	oneminusv  = 1 - v;
        	oneplusv = 1+v;
        	abslogv = abs(log(v));	
        	abslogoneminusv = abs(log(oneminusv));  
        	if (v > 0.5,  
		                 \\ 1 > v >= 1/2 
 				r1 = ceil( (coeff + abslogv +  abs(log (gammaratiosvector[ceilsplusone])) + abs(log (datavector_zeta[ceilsplustwo])) ) / abslogoneminusv );   \\number of summands  
 				if (r1 < 3, r1 = 3);  
 				if (r1 < ceils , r1 = ceils);    
 				start = oneminusv;
 				correction = datavector_zeta[1] + v^minuss ;
 				,
		                 \\ 0 < v < 1/2 
 				r1 = ceil( (coeff + abslogoneminusv + abs(log (gammaratiosvector[ceilsplusone])) + abs(log (datavector_zeta[ceilsplustwo])) )  / abslogv );   \\number of summands  
 				if (r1 < 3, r1 = 3);  
 				if (r1 < ceils , r1 = ceils);   
            	start = minusv;
            	correction =  datavector_zeta[1] + v^minuss + (oneplusv)^minuss  ;
       		 );    

        	u = start; \\ starting point for computing  using a repeated product strategy
        	for (k=1, r1,  
 			  \\S1+= u * gammaratiosvector[k]* datavector_zeta[k+1]; \\ zeta(s) per k=1 
 			  S1+= u * coeffvector[k];
 			  u*= start
 			  );    
       	 S1 = correction + S1 ; 	    \\ correction values 
		); 
 
		S = S + S1;
    	endseriestime += getwalltime()-startseriestime;
    );	
);

elaptime=getwalltime()-starttime; \\ elapsed time 

print("Raw horizontal shift determination time: ", elaptimehorshift, " millisec"); 
print("Raw precomputation time: ", elaptimeprecomp, " millisec");
print("Raw EM comp time: ", endEMtime, " millisec");
print("Raw shift comp time: ", endshifttime, " millisec");
print("Raw series comp time: ", endseriestime, " millisec");
print("Raw total time (no precomp, no horshift det): ", elaptime, " millisec");

return(S);
print("----- ");
}
 
/* *******************************************
 Here we compute the values of Dirichletbeta(s); s > 1,   
 (based on HURWITZ ZETA (reflection formulae) for the length of the sums)
 
 IMPORTANT: to have a correct representation of
 the result, the order of magnitude of the result
 has to be considered
******************************************* */

 \\  dirichlet beta definition for verification
{verifdirichletbeta(s) =  (zetahurwitz(s, 0.25) - zetahurwitz(s, 0.75)) / (4^s);}

 
{dirichletbeta(s, prec, LOOP=1)=  local(S, u , coeff1, aoverq, r1, m, 
datavector_zeta, maxindex, S1, start, aux, starttimeprecomp, elaptimeprecomp, 
oneminusaoverq, abslogaoverq, abslogoneminusaoverq, elaptime, decprec, gammaratiosvector, 
minuss, ceils, starttime, magnitude, logten, logtwo, coeffvector); 

default(format , f);  \\ print numbers with floating point notation (not E notation)
  
if (s <= 1, error(" s input has to be > 1 !!" ));

if (frac(s) == 0, s=floor(s));
			 	
\\if (prec <> 32 && prec <> 64 && prec <> 80 && prec <> 128, error(" prec = 32, 64, 80, 128 (bits) !!!!") );		

logten = log(10);
logtwo = log(2);
decprec =  ceil(prec*logtwo/logten); \\ required decimal precision
\\print("precision: at least ", decprec , " decimal digits");

magnitude = 2*ceil ( s*abs(log(5))/logten ); \\ maximal order of magnitude of the function in base 10 ?? CHECK !!!
magnitude = ceil ( 1.1*magnitude );
\\ 1.1 just to be safer

\\print(magnitude);
decprec = magnitude + decprec + 16; \\ decimal precision  
\\ +16 just to be safer

default( realprecision, decprec );  
print("Internal precision: # decimal digits: ", decprec);  

\\m = prec;
m = ceil(decprec*logten/logtwo);
 
\\ initialisation of the coefficients of the sums

starttimeprecomp = getwalltime();

maxindex = ceil(m+2*s); 
 
datavector_zeta = vector(maxindex,i, zeta(s+2*i-1)-1) ;   
gammaratiosvector = vector(maxindex,i, 0) ;   
gammaratiosvector[1] = s;
for ( i = 2, maxindex, gammaratiosvector[i] =  gammaratiosvector[i-1]*(s+2*i-2)/(2*i-1) *(s+2*i-3)/(2*i-2));
coeffvector = vector(maxindex,i, 0);  
for ( i = 1, maxindex, coeffvector[i] =  gammaratiosvector[i]* datavector_zeta[i] );



elaptimeprecomp=getwalltime()-starttimeprecomp; \\ elapsed time

\\ END initialisation of the coefficients of the sums

starttime = getwalltime(); 

\\sminusone = s-1;
minuss = -s;
 
\\ computing with riemann zeta values 
coeff1 = (m+1)*logtwo; \\ needed for the length of sums 

ceils = ceil(s);
\\ceilsplusone = ceils+1;
\\ceilsplustwo = ceils+2; 
  
for(j=1, LOOP, 
	S1=0;   
	start=1/16;
	u = start; \\ starting point for computing using a repeated product strategy

	aoverq = 1/4;   
	oneminusaoverq  = 3/4;
	abslogaoverq = abs(log(aoverq));	 
	abslogoneminusaoverq = abs(log(oneminusaoverq));  
	aux = abs(log(4*(zeta(5)-1)));

	r1 = ceil ( (coeff1 + abslogoneminusaoverq + aux)/ abslogaoverq); \\ number of summands
	r1 = ceil(r1/2);
	if (r1 < 3, r1 = 3);  
	if (r1 < ceils , r1 = ceils);     

	for (k=1, r1,  
		\\S1+= u * gammaratiosvector[k]* datavector_zeta[k];  
		S1+= u * coeffvector[k];
		u *= start
		);    
	
	S1 = 8 * 4^minuss * S1 ; 

	S = 1 + 5^minuss - 3^minuss - S1;
);

elaptime=getwalltime()-starttime; \\ elapsed time
 
print("Raw precomputation time: ", elaptimeprecomp, " millisec");
print("Raw series computation time: ", elaptime, " millisec");
print("Raw computation time: ", elaptime+elaptimeprecomp, " millisec");
print("----- ");
return(S);
}

/* *******************************************
 Here we compute the values of Dirichletbetaprime(s); s > 1,   
 (based on HURWITZ ZETA (reflection formulae) for the length of the sums)
 
 IMPORTANT: to have a correct representation of
 the result, the order of magnitude of the result
 has to be considered
******************************************* */

 \\ first derivative of dirichlet beta definition for verification
{verifdirichletbetaprime(s) =  (zetahurwitz(s, 0.25, 1) - zetahurwitz(s, 0.75,1)) / (4^s)
							  - 2*log(2)*verifdirichletbeta(s)
\\							  - 2*log(2)*((zetahurwitz(s, 0.25) - zetahurwitz(s, 0.75)) / (4^s))
									;}

 
{dirichletbetaprime(s, prec, LOOP=1)=  local(S, u , coeff1, aoverq, r1, m, betacoeff, aux1, coeff, psi_vector, 
coeffprimevector, datavector_zeta, datavector_zeta_prime, maxindex, S1, start, aux, starttimeprecomp, 
elaptimeprecomp, boundindex, oneminusaoverq, abslogaoverq, abslogoneminusaoverq, elaptime, decprec, 
gammaratiosvector, minuss, ceilsplusone, ceils, starttime, magnitude, logten, logtwo); 

default(format , f);  \\ print numbers with floating point notation (not E notation)
  
if (s <= 1, error(" s input has to be > 1 !!" ));
if (frac(s) == 0, s=floor(s));
			 	
if (prec <> 32 && prec <> 64 && prec <> 80 && prec <> 128, error(" prec = 32, 64, 80, 128 (bits) !!!!") );		

logten = log(10);
logtwo = log(2);
decprec =  ceil(prec*logtwo/logten); \\ required decimal precision
\\ print("precision: at least ", decprec , " decimal digits");

magnitude = ceil ( s*abs(log(5))/logten ); \\ maximal order of magnitude in base 10
magnitude = ceil ( 1.1*magnitude );
\\ 1.1 just to be safer

\\print(magnitude);
decprec = magnitude + decprec + 16; \\ decimal precision  
\\ +16 just to be safer

default( realprecision, decprec );  
print("Internal precision: # decimal digits: ", decprec);  

\\m = prec;
m = ceil(decprec*logten/logtwo);	
 
\\ initialisation of the coefficients of the sums

starttimeprecomp = getwalltime();

coeff1 = (m+1)*logtwo; \\ needed for the length of sums 
ceils = ceil(s);
ceilsplusone = ceils+1;
\\ceilsplustwo = ceils+2; 

coeff = (m+1+2*ceils)*logtwo;  

betacoeff = abs(lngamma(ceilsplusone+s)-lngamma(s)-lngamma(ceilsplusone+1));

aux1 = betacoeff + abs(log(zeta(ceilsplusone+s)-1));  
boundindex = aux1 / logtwo ;
boundindex = ceil(boundindex); 
maxindex = m + 10 + boundindex;
 
 
datavector_zeta = vector(maxindex,i, zeta(s+2*i-1)-1) ;   
datavector_zeta_prime = vector(maxindex,i, zetahurwitz(s+2*i-1,1,1)) ;   
gammaratiosvector = vector(maxindex,i, 0) ;   
gammaratiosvector[1] = s;
for ( i = 2, maxindex, gammaratiosvector[i] =  gammaratiosvector[i-1]*(s+2*i-2)/(2*i-1) *(s+2*i-3)/(2*i-2)); 

psi_vector = vector(maxindex,i, 0) ;   
psi_vector[1] = 1/s - 2*log(2); \\psi(s+1) - psi(s) - 2*log(2)
for ( i = 2, maxindex, psi_vector[i] =  psi_vector[i-1] + 1/(s+2*i-3) + 1/(s+2*i-2) ); \\ psi(s+2*i-1) - psi(s) - 2*log(2)
 
coeffprimevector = vector(maxindex,i, 0) ; 
for ( i = 1, maxindex, coeffprimevector[i]=gammaratiosvector[i]* (datavector_zeta[i] * psi_vector[i] + datavector_zeta_prime[i]) );

elaptimeprecomp=getwalltime()-starttimeprecomp; \\ elapsed time

\\ END initialisation of the coefficients of the sums

starttime = getwalltime(); 

\\sminusone = s-1;
minuss = -s;
 
\\ computing with riemann zeta values

for(j=1, LOOP,   
	S1=0;   
	start=1./16;
	u = start; \\ starting point for computing  using a repeated product strategy

	aoverq = 0.25;   
	oneminusaoverq  = 0.75;
	abslogaoverq = abs(log(aoverq)); 
	abslogoneminusaoverq = abs(log(oneminusaoverq));  
	aux = coeff + betacoeff;

	r1 = ceil ( (coeff1 + abslogoneminusaoverq + aux)/ abslogaoverq); \\ number of summands
	r1 = ceil(r1/2);
	if (r1 < 3, r1 = 3);  
	if (r1 < ceils , r1 = ceils);     
 
	for (k=1, r1,  
		\\S1+= u * gammaratiosvector[k]* (datavector_zeta[k] * psi_vector[k] + datavector_zeta_prime[k]); 
		S1+= u * coeffprimevector[k];  
		u *= start
	);    
	
	S1 = 8 * 4^minuss * S1 ; 

	S = - log(5)*5^minuss + log(3)*3^minuss - S1;
);

elaptime=getwalltime()-starttime; \\ elapsed time
 
print("Raw precomputation time: ", elaptimeprecomp, " millisec");
print("Raw series computation time: ", elaptime, " millisec");
print("Raw computation time: ", elaptime+elaptimeprecomp, " millisec");
print("----- ");
return(S);
}

/* *******************************************
 Here we compute at the same time the values of 
 beta(s), beta'(s), beta'(s)/beta(s); s > 1,   
 beta = Dirichlet beta
 (based on HURWITZ ZETA (reflection formulae) for the length of the sums)
 
 IMPORTANT: to have a correct representation of
 the result, the order of magnitude of the result
 has to be considered
******************************************* */

 \\ first derivative of dirichlet beta definition for verification
{verifdirichletbetalogder(s) = -2*log(2) + (zetahurwitz(s, 0.25, 1) - zetahurwitz(s, 0.75,1)) /(zetahurwitz(s, 0.25) - zetahurwitz(s, 0.75))
									;}

{dirichletbetavalues(s, prec, LOOP=1)=  local(S, u , coeff1, aoverq, r1, m, betacoeff, aux1, coeff, psi_vector, num,denom,
datavector_zeta, datavector_zeta_prime, maxindex, S1num, S1denom, start, aux, starttimeprecomp, elaptimeprecomp, boundindex,
oneminusaoverq, abslogaoverq, abslogoneminusaoverq, elaptime, decprec, gammaratiosvector, logtwo, logfour,
minuss, ceilsplusone, ceils, starttime, magnitude, corrpowfour, powfive, coeffvector, coeffprimevector,
powthree, powfour, sumpowers, logten, betas, betaprimes); 

default(format , f);  \\ print numbers with floating point notation (not E notation)
  
if (s <= 1, error(" s input has to be > 1 !!" ));

if (frac(s) == 0, s=floor(s));
			 	
if (prec <> 32 && prec <> 64 && prec <> 80 && prec <> 128, error(" prec = 32, 64, 80, 128 (bits) !!!!") );		

logten = log(10);
logtwo = log(2);
decprec =  ceil(prec*logtwo/logten); \\ required decimal precision
\\ print("precision: at least ", decprec , " decimal digits");

magnitude = ceil ( s*abs(log(5))/logten ); \\ maximal order of magnitude in base 10
magnitude = ceil ( 1.1*magnitude );
\\ 1.1 just to be safer

\\print(magnitude);
decprec = magnitude + decprec + 16; \\ decimal precision  
\\ +16 just to be safer

default( realprecision, decprec );  
print("Internal precision: # decimal digits: ", decprec);  

\\m = prec;
m = ceil(decprec*logten/logtwo);	

\\ initialisation of the coefficients of the sums

starttimeprecomp = getwalltime();

logtwo = log(2);
logfour = 2*logtwo;

coeff1 = (m+1)*logtwo; \\ needed for the length of sums 
ceils = ceil(s);
ceilsplusone = ceils+1;
\\ceilsplustwo = ceils+2; 

coeff = (m+1+2*ceils)*logtwo;  

betacoeff = abs(lngamma(ceilsplusone+s)-lngamma(s)-lngamma(ceilsplusone+1));

aux1 = betacoeff + abs(log(zeta(ceilsplusone+s)-1));  
boundindex = aux1 / logtwo ;
boundindex = ceil(boundindex); 
maxindex = m + 10 + boundindex;
 
 
datavector_zeta = vector(maxindex,i, zeta(s+2*i-1)-1) ;  
datavector_zeta_prime = vector(maxindex,i, zetahurwitz(s+2*i-1,1,1)) ;    

gammaratiosvector = vector(maxindex,i, 0) ;   
gammaratiosvector[1] = s;
for ( i = 2, maxindex, gammaratiosvector[i] =  gammaratiosvector[i-1]*(s+2*i-2)/(2*i-1) *(s+2*i-3)/(2*i-2)); 

psi_vector = vector(maxindex,i, 0) ;  
psi_vector[1] = 1/s; \\psi(s+1) - psi(s)
for ( i = 2, maxindex, psi_vector[i] =  psi_vector[i-1] + 1/(s+2*i-3) + 1/(s+2*i-2) ); \\ psi(s+2*i-1) - psi(s) 

coeffvector = vector(maxindex,i, 0);  
for ( i = 1, maxindex, coeffvector[i] =  gammaratiosvector[i]* datavector_zeta[i] );

coeffprimevector = vector(maxindex,i, 0) ; 
for ( i = 1, maxindex, coeffprimevector[i] = gammaratiosvector[i]* (datavector_zeta[i] * psi_vector[i] + datavector_zeta_prime[i]) ); 
elaptimeprecomp=getwalltime()-starttimeprecomp; \\ elapsed time

\\ END initialisation of the coefficients of the sums

starttime = getwalltime(); 

\\sminusone = s-1;
minuss = -s;

powfive = 5^minuss;
powthree = 3^minuss;
powfour = 4^minuss;

corrpowfour = 8 * powfour;
sumpowers = 1 - powthree + powfive;

\\ computing with riemann zeta values

for(j=1, LOOP, 
  
	S1num=0;   
	S1denom=0;   
	start=1./16;
	u = start; \\ starting point for computing using a repeated product strategy

	aoverq = 0.25;   
	oneminusaoverq  = 0.75;
	abslogaoverq = abs(log(aoverq));	 
	abslogoneminusaoverq = abs(log(oneminusaoverq));  
	aux = coeff + betacoeff;

	r1 = ceil ( (coeff1 + abslogoneminusaoverq + aux)/ abslogaoverq); \\ number of summands
	r1 = ceil(r1/2);
	if (r1 < 3, r1 = 3);  
	if (r1 < ceils , r1 = ceils);     
 
	for (k=1, r1,  
		\\S1num += u * gammaratiosvector[k] * (datavector_zeta[k] * psi_vector[k] + datavector_zeta_prime[k]); 
		\\S1denom += u * gammaratiosvector[k] * datavector_zeta[k];
		S1num += u * coeffprimevector[k]; 
		S1denom += u * coeffvector[k];
		u *= start
		);    
	
	num = logfour*sumpowers + log(3)*powthree - log(5)*powfive - corrpowfour * S1num ;  
	denom = sumpowers - corrpowfour * S1denom ; 
	
	S = -logfour + num/denom; \\ betaprimes/betas

	betas = denom;	
	betaprimes =  S * betas;
);

elaptime=getwalltime()-starttime; \\ elapsed time
 
print("Raw precomputation time: ", elaptimeprecomp, " millisec");
print("Raw series computation time: ", elaptime, " millisec");
print("Raw computation time: ", elaptime+elaptimeprecomp, " millisec");
print("----- ");
print("beta(s) = ", betas);
print("betaprime(s) = ", betaprimes);
print("betaprime/beta(s) = ", S);
\\return(S);
}


/* Catalan constant computation */
 
{my_catalan(decprec)=
 local(aoverq, S, oneminusaoverq, maxr1,
aux, u , coeff, coeff1, r1, m, minutes, millisec, abslogaoverq, abslogoneminusaoverq,
seconds, step, elaptimecatalan, elaptimestart, start, correction); 

default(format , f);  \\ print numbers with floating point notation (not E notation)

\\s=2;
coeff = decprec + 10 ; \\ decimal precision  
m = ceil(coeff*log(10)/log(2));\\ binary precision  
coeff1 = (m+1)*log(2);

default(realprecision, ceil(coeff));   

\\zetas = Pi^2/6; \\ zeta(2);

elaptimestart=gettime();

 
maxr1 = 0;

print("Catalan constant computation"); 

aoverq = 0.25; 
S=0;  
oneminusaoverq  = 0.75;
abslogaoverq = abs(log(aoverq));
abslogoneminusaoverq = abs(log(oneminusaoverq));   
aux = abs(log(4*(zeta(5)-1)));
r1 = ceil ( (coeff1 + abslogoneminusaoverq + aux)/ abslogaoverq); \\ number of summands
r1 = ceil(r1/2);
start = 1/16;
step = 1/16;
correction = 209/225;

if (r1>maxr1, maxr1 = r1);
u = start; \\ starting point for computing using a repeated product strategy
for(k=1, r1, 
      S+= u * (zeta(2*k+1)-1) * k; 
      u*= step;
 		  );   	

S = correction - S; 	    \\ correction values

elaptimecatalan=gettime()-elaptimestart; 

/* print computation time */
seconds=floor(elaptimecatalan/1000)%60;
minutes=floor(elaptimecatalan/60000);
millisec=elaptimecatalan- minutes*60000 - seconds*1000;
print("Catalan constant computation time: ", minutes, " min, ", seconds, " sec, ", millisec, " millisec");	
return(S);
 }
 
/************************************  


OPTIPLEX

gp2c-run -pmy_ -g -W  specialfunct_package-final-v3.gp   

GP/PARI CALCULATOR Version 2.13.3 (released)
amd64 running linux (x86-64/GMP-6.2.0 kernel) 64-bit version
compiled: Oct 25 2021, gcc version 9.3.0 (Ubuntu 9.3.0-17ubuntu1~20.04)
threading engine: single
(readline v7.0 enabled, extended help enabled)
Copyright (C) 2000-2020 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 ?17 for how to get moral (and possibly technical) support.

parisize = 8000000, primelimit = 500000


-------------------------------------------------------------------
loggamma verification

? v=[32,64,80,128]
%1 = [32, 64, 80, 128]
? for (i=1,4, print("----- "); print("binary precision = ",v[i]); print(loggamma(1345.1234,v[i])-lngamma(1345.1234)))
----- 
binary precision = 32
Internal precision: # decimal digits: 31
Evaluation of the horizontal shift
horizontal shift = 11
Length of the EM-sum = 34
Raw horizontal shift determination time: 0 millisec
Raw precomputation time: 22 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.000000000000000000000000000013226275744830732540
----- 
binary precision = 64
Internal precision: # decimal digits: 41
Evaluation of the horizontal shift
horizontal shift = 15
Length of the EM-sum = 47
Raw horizontal shift determination time: 0 millisec
Raw precomputation time: 1 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
-0.000000000000000000000000000000000003636523115995667317
----- 
binary precision = 80
Internal precision: # decimal digits: 46
Evaluation of the horizontal shift
horizontal shift = 17
Length of the EM-sum = 53
Raw horizontal shift determination time: 0 millisec
Raw precomputation time: 1 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
-0.0000000000000000000000000000000000000000000004325558347957441508
----- 
binary precision = 128
Internal precision: # decimal digits: 60
Evaluation of the horizontal shift
horizontal shift = 22
Length of the EM-sum = 69
Raw horizontal shift determination time: 1 millisec
Raw precomputation time: 1 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
-0.00000000000000000000000000000000000000000000000000000008670992117257616652



performances

? loggamma(0.01,128,10000)
Internal precision: # decimal digits: 57
Evaluation of the horizontal shift
horizontal shift = 21
Length of the EM-sum = 66
Raw horizontal shift determination time: 1 millisec
Raw precomputation time: 3 millisec
Raw EM comp time: 2 millisec
Raw shift comp time: 2 millisec
Raw series comp time: 144 millisec
Raw total time (no precomp, no horshift det): 149 millisec
%5 = 4.59947987804202172251394541100874808726100141338528965242


? for(j=1,10000,lngamma(0.01))
? ##
  ***   last result computed in 289 ms.
? for(j=1,10000,log(gamma(0.01)))
time = 292 ms.
?  lngamma(0.01)
%5 = 4.59947987804202172251394541100874808726100141338528965242 + 0.000000000000000000000000000000000000000000000000000000000*I


? loggamma(100.01,128,10000)
Internal precision: # decimal digits: 59
Evaluation of the horizontal shift
horizontal shift = 22
Length of the EM-sum = 69
Raw horizontal shift determination time: 1 millisec
Raw precomputation time: 0 millisec
Raw EM comp time: 333 millisec
Raw shift comp time: 60 millisec
Raw series comp time: 108 millisec
Raw total time (no precomp, no horshift det): 501 millisec
%4 = 359.1802074905942794959513321973192358382494167931552488267


?  for(j=1,10000,log(gamma(100.01)))
time = 244 ms.
?  for(j=1,10000,lngamma(100.01))
time = 256 ms.
? log(gamma(100.01))
%13 = 359.18020749059427949595133219731923583824941679315524882672
-----------------------------------------------------------------
digamma verification

? v=[32,64,80,128]
%1 = [32, 64, 80, 128]
? for (i=1,4, print("----- "); print("binary precision = ",v[i]); print(digamma(1345.1234,v[i])-psi(1345.1234)))
----- 
binary precision = 32
Internal precision: # decimal digits: 28
Evaluation of the horizontal shift
horizontal shift = 11
Length of the EM-sum = 34
Raw horizontal shift determination time: 0 millisec
Raw precomputation time: 1 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
-0.00000000000000000000000000014905274838347512527
----- 
binary precision = 64
Internal precision: # decimal digits: 38
Evaluation of the horizontal shift
horizontal shift = 14
Length of the EM-sum = 43
Raw horizontal shift determination time: 0 millisec
Raw precomputation time: 0 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.0000000000000000000000000000000000001880790961315660013
----- 
binary precision = 80
Internal precision: # decimal digits: 43
Evaluation of the horizontal shift
horizontal shift = 16
Length of the EM-sum = 50
Raw horizontal shift determination time: 1 millisec
Raw precomputation time: 0 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
-0.0000000000000000000000000000000000003291384182302405022
----- 
binary precision = 128
Internal precision: # decimal digits: 57
Evaluation of the horizontal shift
horizontal shift = 21
Length of the EM-sum = 65
Raw horizontal shift determination time: 1 millisec
Raw precomputation time: 1 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
-0.0000000000000000000000000000000000000000000000000000000012744735289059618216


performances

? digamma(0.01,128,10000)
Internal precision: # decimal digits: 58
Evaluation of the horizontal shift
horizontal shift = 22
Length of the EM-sum = 68
Raw horizontal shift determination time: 2 millisec
Raw precomputation time: 3 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 115 millisec
Raw total time (no precomp, no horshift det): 119 millisec
%13 = -100.5608854578686744974809613039329488963941105697188339255

? for(j=1,10000,psi(0.01))
time = 255 ms.

? psi(0.01)
%17 = -100.5608854578686744974809613039329488963941105697188339255


?  digamma(100.01,128,10000)
Internal precision: # decimal digits: 59
Evaluation of the horizontal shift
horizontal shift = 22
Length of the EM-sum = 68
Raw horizontal shift determination time: 1 millisec
Raw precomputation time: 0 millisec
Raw EM comp time: 237 millisec
Raw shift comp time: 58 millisec
Raw series comp time: 134 millisec
Raw total time (no precomp, no horshift det): 431 millisec
%4 = 4.6002623493548090854079423675842515513857094155170245178891

 
? for(j=1,10000,psi(100.01))
time = 172 ms.
? psi(100.01)
%20 = 4.6002623493548090854079423675842515513857094155170245178894


-------------------------------------------------------------------
BATEMAN G verification

? v=[32,64,80,128]
%1 = [32, 64, 80, 128]

for (i=1,4, print("----- "); print("binary precision = ",v[i]); print(batemanG(1345.1234,v[i])-verifbatemanG(1345.1234)))
----- 
binary precision = 32
Internal precision: # decimal digits: 29
Evaluation of the horizontal shift
horizontal shift = 11
Length of the EM-sum = 34
Raw horizontal shift determination time: 0 millisec
Raw precomputation time: 0 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.000000000000000000000000000004579690678953333430
----- 
binary precision = 64
Internal precision: # decimal digits: 39
Evaluation of the horizontal shift
horizontal shift = 15
Length of the EM-sum = 46
Raw horizontal shift determination time: 0 millisec
Raw precomputation time: 0 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.00000000000000000000000000000000000015196191740818102010
----- 
binary precision = 80
Internal precision: # decimal digits: 44
Evaluation of the horizontal shift
horizontal shift = 17
Length of the EM-sum = 53
Raw horizontal shift determination time: 0 millisec
Raw precomputation time: 1 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
-0.00000000000000000000000000000000000000000000015462531523940791522
----- 
binary precision = 128
Internal precision: # decimal digits: 58
Evaluation of the horizontal shift
horizontal shift = 22
Length of the EM-sum = 68
Raw horizontal shift determination time: 1 millisec
Raw precomputation time: 1 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.000000000000000000000000000000000000000000000000000000005440626682052049666



? for (i=1,4, print("----- "); print("binary precision = ",v[i]); print(batemanG(1345.53674,v[i])-verifbatemanG(1345.53674)))
----- 
binary precision = 32
Internal precision: # decimal digits: 29
Evaluation of the horizontal shift
horizontal shift = 11
Length of the EM-sum = 34
Raw horizontal shift determination time: 0 millisec
Raw precomputation time: 0 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.0000000000000000000000000000011999680256569724517
----- 
binary precision = 64
Internal precision: # decimal digits: 39
Evaluation of the horizontal shift
horizontal shift = 15
Length of the EM-sum = 46
Raw horizontal shift determination time: 0 millisec
Raw precomputation time: 1 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.00000000000000000000000000000000000008196318032100715632
----- 
binary precision = 80
Internal precision: # decimal digits: 44
Evaluation of the horizontal shift
horizontal shift = 17
Length of the EM-sum = 53
Raw horizontal shift determination time: 0 millisec
Raw precomputation time: 1 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
-0.00000000000000000000000000000000000000000000004095822787236731653
----- 
binary precision = 128
Internal precision: # decimal digits: 58
Evaluation of the horizontal shift
horizontal shift = 22
Length of the EM-sum = 68
Raw horizontal shift determination time: 1 millisec
Raw precomputation time: 1 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.000000000000000000000000000000000000000000000000000000006734858284465220613


? for (i=1,4, print("----- "); print("binary precision = ",v[i]); print(batemanG(1344.53674,v[i])-verifbatemanG(1344.53674))) 
----- 
binary precision = 32
Internal precision: # decimal digits: 29
Evaluation of the horizontal shift
horizontal shift = 11
Length of the EM-sum = 34
Raw horizontal shift determination time: 0 millisec
Raw precomputation time: 0 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.00000000000000000000000000003205396359116731903
----- 
binary precision = 64
Internal precision: # decimal digits: 39
Evaluation of the horizontal shift
horizontal shift = 15
Length of the EM-sum = 46
Raw horizontal shift determination time: 0 millisec
Raw precomputation time: 0 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.00000000000000000000000000000000000002519736425834883867
----- 
binary precision = 80
Internal precision: # decimal digits: 44
Evaluation of the horizontal shift
horizontal shift = 17
Length of the EM-sum = 53
Raw horizontal shift determination time: 1 millisec
Raw precomputation time: 0 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
-0.0000000000000000000000000000000000000000000010482656037223930473
----- 
binary precision = 128
Internal precision: # decimal digits: 58
Evaluation of the horizontal shift
horizontal shift = 22
Length of the EM-sum = 68
Raw horizontal shift determination time: 1 millisec
Raw precomputation time: 1 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.0000000000000000000000000000000000000000000000000000000023633456271497150918

? for (i=1,4, print("----- "); print("binary precision = ",v[i]); print(batemanG(1344.2345,v[i])-verifbatemanG(1344.2345)))  
----- 
binary precision = 32
Internal precision: # decimal digits: 29
Evaluation of the horizontal shift
horizontal shift = 11
Length of the EM-sum = 34
Raw horizontal shift determination time: 0 millisec
Raw precomputation time: 0 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 1 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 1 millisec
0.00000000000000000000000000008923493696518732259
----- 
binary precision = 64
Internal precision: # decimal digits: 39
Evaluation of the horizontal shift
horizontal shift = 15
Length of the EM-sum = 46
Raw horizontal shift determination time: 0 millisec
Raw precomputation time: 0 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.0000000000000000000000000000000000002208184505118898682
----- 
binary precision = 80
Internal precision: # decimal digits: 44
Evaluation of the horizontal shift
horizontal shift = 17
Length of the EM-sum = 53
Raw horizontal shift determination time: 0 millisec
Raw precomputation time: 0 millisec
Raw EM comp time: 1 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 1 millisec
-0.000000000000000000000000000000000000000000002847660703916066829
----- 
binary precision = 128
Internal precision: # decimal digits: 58
Evaluation of the horizontal shift
horizontal shift = 22
Length of the EM-sum = 68
Raw horizontal shift determination time: 0 millisec
Raw precomputation time: 1 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.000000000000000000000000000000000000000000000000000000015655861836043060307

performances

? batemanG(1344.2345,128, 10000)
Internal precision: # decimal digits: 58
Evaluation of the horizontal shift
horizontal shift = 22
Length of the EM-sum = 68
Raw horizontal shift determination time: 2 millisec
Raw precomputation time: 3 millisec
Raw EM comp time: 324 millisec
Raw shift comp time: 104 millisec
Raw series comp time: 211 millisec
Raw total time (no precomp, no horshift det): 652 millisec
%18 = 0.0007441945276717117527836797619885823490444664194851457096284

? for(j=1,10000,verifbatemanG(1344.2345))
time = 345 ms.

? verifbatemanG(1344.2345)
%27 = 0.0007441945276717117527836797619885823490444664194851457095762

? batemanG(0.2345,128, 10000)
Internal precision: # decimal digits: 57
Evaluation of the horizontal shift
horizontal shift = 21
Length of the EM-sum = 65
Raw horizontal shift determination time: 2 millisec
Raw precomputation time: 2 millisec
Raw EM comp time: 1 millisec
Raw shift comp time: 2 millisec
Raw series comp time: 183 millisec
Raw total time (no precomp, no horshift det): 186 millisec
%21 = 7.44875088476837717985432385159721060725696200745027313898


? for(j=1,10000,verifbatemanG(0.2345))
time = 352 ms.

? verifbatemanG(0.2345)
%32 = 7.44875088476837717985432385159721060725696200745027313898

? batemanG(0.7345,128, 10000)
Internal precision: # decimal digits: 57
Evaluation of the horizontal shift
horizontal shift = 21
Length of the EM-sum = 65
Raw horizontal shift determination time: 1 millisec
Raw precomputation time: 1 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 2 millisec
Raw series comp time: 185 millisec
Raw total time (no precomp, no horshift det): 191 millisec
%14 = 1.99878080384105851086077727818910041831892148122283220219

? for(j=1,10000,verifbatemanG(0.7345))
time = 354 ms.
? verifbatemanG(0.7345)
%6 = 1.99878080384105851086077727818910041831892148122283220219

-------------------------------------------------------------------
zeta hurwitz verification

? v=[32,64,80,128]
%1 = [32, 64, 80, 128]
? for (i=1,4, print("----- "); print("binary precision = ",v[i]); print(hurwitzzeta(1.1,0.1234,v[i])-zetahurwitz(1.1,0.1234)))
----- 
binary precision = 32
Internal precision: # decimal digits: 28
Evaluation of the horizontal shift
horizontal shift = 11
Length of the EM-sum = 34
Raw horizontal shift determination time: 1 millisec
Raw precomputation time: 13 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.0000000000000000000000000000000000001880790961315660013
----- 
binary precision = 64
Internal precision: # decimal digits: 38
Evaluation of the horizontal shift
horizontal shift = 14
Length of the EM-sum = 43
Raw horizontal shift determination time: 0 millisec
Raw precomputation time: 8 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.0000000000000000000000000000000000005642372883946980038
----- 
binary precision = 80
Internal precision: # decimal digits: 43
Evaluation of the horizontal shift
horizontal shift = 16
Length of the EM-sum = 49
Raw horizontal shift determination time: 0 millisec
Raw precomputation time: 8 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.0000000000000000000000000000000000003761581922631320025
----- 
binary precision = 128
Internal precision: # decimal digits: 57
Evaluation of the horizontal shift
horizontal shift = 21
Length of the EM-sum = 65
Raw horizontal shift determination time: 1 millisec
Raw precomputation time: 14 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.00000000000000000000000000000000000000000000000000000003058736469374308372


? for (i=1,4, print("----- "); print("binary precision = ",v[i]); print(hurwitzzeta(1.1,0.6234,v[i])-zetahurwitz(1.1,0.6234)))
----- 
binary precision = 32
Internal precision: # decimal digits: 28
Evaluation of the horizontal shift
horizontal shift = 11
Length of the EM-sum = 34
Raw horizontal shift determination time: 1 millisec
Raw precomputation time: 13 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.00000000000000000000000000000000000014105932209867450096
----- 
binary precision = 64
Internal precision: # decimal digits: 38
Evaluation of the horizontal shift
horizontal shift = 14
Length of the EM-sum = 43
Raw horizontal shift determination time: 0 millisec
Raw precomputation time: 8 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.0000000000000000000000000000000000005172175143618065035
----- 
binary precision = 80
Internal precision: # decimal digits: 43
Evaluation of the horizontal shift
horizontal shift = 16
Length of the EM-sum = 49
Raw horizontal shift determination time: 0 millisec
Raw precomputation time: 8 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.0000000000000000000000000000000000003761581922631320025
----- 
binary precision = 128
Internal precision: # decimal digits: 57
Evaluation of the horizontal shift
horizontal shift = 21
Length of the EM-sum = 65
Raw horizontal shift determination time: 1 millisec
Raw precomputation time: 14 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.00000000000000000000000000000000000000000000000000000002548947057811923643


? for (i=1,4, print("----- "); print("binary precision = ",v[i]); print(hurwitzzeta(2,0.6234,v[i])-zetahurwitz(2,0.6234)))
----- 
binary precision = 32
Internal precision: # decimal digits: 28
Evaluation of the horizontal shift
horizontal shift = 11
Length of the EM-sum = 33
Raw horizontal shift determination time: 1 millisec
Raw precomputation time: 0 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
-0.00000000000000000000000000000000000002350988701644575016
----- 
binary precision = 64
Internal precision: # decimal digits: 38
Evaluation of the horizontal shift
horizontal shift = 15
Length of the EM-sum = 46
Raw horizontal shift determination time: 0 millisec
Raw precomputation time: 1 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
-0.00000000000000000000000000000000000002350988701644575016
----- 
binary precision = 80
Internal precision: # decimal digits: 43
Evaluation of the horizontal shift
horizontal shift = 16
Length of the EM-sum = 49
Raw horizontal shift determination time: 0 millisec
Raw precomputation time: 1 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
-0.00000000000000000000000000000000000002350988701644575016
----- 
binary precision = 128
Internal precision: # decimal digits: 57
Evaluation of the horizontal shift
horizontal shift = 22
Length of the EM-sum = 68
Raw horizontal shift determination time: 1 millisec
Raw precomputation time: 1 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
-0.0000000000000000000000000000000000000000000000000000000019117102933589427324
? 



? for (i=1,4, print("----- "); print("binary precision = ",v[i]); print(hurwitzzeta(11.6,0.6234,v[i])-zetahurwitz(11.6,0.6234)))
----- 
binary precision = 32
Internal precision: # decimal digits: 30
Evaluation of the horizontal shift
horizontal shift = 12
Length of the EM-sum = 32
Raw horizontal shift determination time: 1 millisec
Raw precomputation time: 9 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.000000000000000000000000000000
----- 
binary precision = 64
Internal precision: # decimal digits: 40
Evaluation of the horizontal shift
horizontal shift = 16
Length of the EM-sum = 44
Raw horizontal shift determination time: 1 millisec
Raw precomputation time: 5 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.0000000000000000000000000000000000000000
----- 
binary precision = 80
Internal precision: # decimal digits: 45
Evaluation of the horizontal shift
horizontal shift = 17
Length of the EM-sum = 47
Raw horizontal shift determination time: 1 millisec
Raw precomputation time: 10 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.000000000000000000000000000000000000000000000
----- 
binary precision = 128
Internal precision: # decimal digits: 59
Evaluation of the horizontal shift
horizontal shift = 22
Length of the EM-sum = 63
Raw horizontal shift determination time: 1 millisec
Raw precomputation time: 10 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.0000000000000000000000000000000000000000000000000000004486146821748985612


? for (i=1,4, print("----- "); print("binary precision = ",v[i]); print(hurwitzzeta(11.6,0.6234,v[i])-zetahurwitz(11.6,0.6234)))
----- 
binary precision = 32
Internal precision: # decimal digits: 30
Evaluation of the horizontal shift
horizontal shift = 12
Length of the EM-sum = 32
Raw horizontal shift determination time: 0 millisec
Raw precomputation time: 16 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.000000000000000000000000000000
----- 
binary precision = 64
Internal precision: # decimal digits: 40
Evaluation of the horizontal shift
horizontal shift = 16
Length of the EM-sum = 44
Raw horizontal shift determination time: 0 millisec
Raw precomputation time: 6 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.0000000000000000000000000000000000000000
----- 
binary precision = 80
Internal precision: # decimal digits: 45
Evaluation of the horizontal shift
horizontal shift = 17
Length of the EM-sum = 47
Raw horizontal shift determination time: 1 millisec
Raw precomputation time: 9 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.000000000000000000000000000000000000000000000
----- 
binary precision = 128
Internal precision: # decimal digits: 59
Evaluation of the horizontal shift
horizontal shift = 22
Length of the EM-sum = 63
Raw horizontal shift determination time: 1 millisec
Raw precomputation time: 10 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.0000000000000000000000000000000000000000000000000000004486146821748985612


? for (i=1,4, print("----- "); print("binary precision = ",v[i]); print(hurwitzzeta(11.6,0.1234,v[i])-zetahurwitz(11.6,0.1234))) 
----- 
binary precision = 32
Internal precision: # decimal digits: 39
Evaluation of the horizontal shift
horizontal shift = 15
Length of the EM-sum = 41
Raw horizontal shift determination time: 1 millisec
Raw precomputation time: 15 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
-0.00000000000000000000000000000000000000000000010947644252537633367
----- 
binary precision = 64
Internal precision: # decimal digits: 49
Evaluation of the horizontal shift
horizontal shift = 19
Length of the EM-sum = 54
Raw horizontal shift determination time: 0 millisec
Raw precomputation time: 10 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 1 millisec
Raw total time (no precomp, no horshift det): 1 millisec
0.0000000000000000000000000000000000000000000000000
----- 
binary precision = 80
Internal precision: # decimal digits: 54
Evaluation of the horizontal shift
horizontal shift = 21
Length of the EM-sum = 60
Raw horizontal shift determination time: 1 millisec
Raw precomputation time: 10 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.000000000000000000000000000000000000000000000000000000
----- 
binary precision = 128
Internal precision: # decimal digits: 68
Evaluation of the horizontal shift
horizontal shift = 26
Length of the EM-sum = 76
Raw horizontal shift determination time: 1 millisec
Raw precomputation time: 11 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.00000000000000000000000000000000000000000000010947644252537633367



?  for (i=1,4, print("----- "); print("binary precision = ",v[i]); print(hurwitzzeta(2,1345.1234,v[i])-zetahurwitz(2,1345.1234)))
----- 
binary precision = 32
Internal precision: # decimal digits: 34
Evaluation of the horizontal shift
horizontal shift = 13
Length of the EM-sum = 40
Raw horizontal shift determination time: 0 millisec
Raw precomputation time: 0 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.0000000000000000000000000000000017776049336060739151
----- 
binary precision = 64
Internal precision: # decimal digits: 44
Evaluation of the horizontal shift
horizontal shift = 17
Length of the EM-sum = 52
Raw horizontal shift determination time: 0 millisec
Raw precomputation time: 1 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.0000000000000000000000000000000000028781101002951698842
----- 
binary precision = 80
Internal precision: # decimal digits: 49
Evaluation of the horizontal shift
horizontal shift = 19
Length of the EM-sum = 58
Raw horizontal shift determination time: 1 millisec
Raw precomputation time: 1 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.00000000000000000000000000000000000000000000000005462801257545052930
----- 
binary precision = 128
Internal precision: # decimal digits: 63
Evaluation of the horizontal shift
horizontal shift = 24
Length of the EM-sum = 74
Raw horizontal shift determination time: 1 millisec
Raw precomputation time: 1 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.00000000000000000000000000000000000000000000000000000022309525328670387673


? for (i=1,4, print("----- "); print("binary precision = ",v[i]); print(hurwitzzeta(2,1345.6234,v[i])-zetahurwitz(2,1345.6234)))
----- 
binary precision = 32
Internal precision: # decimal digits: 34
Evaluation of the horizontal shift
horizontal shift = 13
Length of the EM-sum = 40
Raw horizontal shift determination time: 0 millisec
Raw precomputation time: 0 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.00000000000000000000000000000000006474787040278542003
----- 
binary precision = 64
Internal precision: # decimal digits: 44
Evaluation of the horizontal shift
horizontal shift = 17
Length of the EM-sum = 52
Raw horizontal shift determination time: 0 millisec
Raw precomputation time: 1 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.00000000000000000000000000000000000003648165084876200879
----- 
binary precision = 80
Internal precision: # decimal digits: 49
Evaluation of the horizontal shift
horizontal shift = 19
Length of the EM-sum = 58
Raw horizontal shift determination time: 1 millisec
Raw precomputation time: 1 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.0000000000000000000000000000000000000000000000000022026476856997827534
----- 
binary precision = 128
Internal precision: # decimal digits: 63
Evaluation of the horizontal shift
horizontal shift = 24
Length of the EM-sum = 74
Raw horizontal shift determination time: 1 millisec
Raw precomputation time: 1 millisec
Raw EM comp time: 0 millisec
Raw shift comp time: 0 millisec
Raw series comp time: 0 millisec
Raw total time (no precomp, no horshift det): 0 millisec
0.000000000000000000000000000000000000000000000000000000006287267910605057995


? hurwitzzeta(20,1345.1234,128,10000)
Internal precision: # decimal digits: 125
Evaluation of the horizontal shift
horizontal shift = 45
Length of the EM-sum = 131
Raw horizontal shift determination time: 2 millisec
Raw precomputation time: 6 millisec
Raw EM comp time: 690 millisec
Raw shift comp time: 319 millisec
Raw series comp time: 535 millisec
Raw total time (no precomp, no horshift det): 1553 millisec
%29 = 0.0000000000000000000000000000000000000000000000000000000000001895963461786327125428610368722191005815672395983749858921

?  for(j=1,10000,zetahurwitz(20,1345.1234))
  ***   last result computed in 8,800 ms.
  
? zetahurwitz(20,1345.1234)
%24 = 0.0000000000000000000000000000000000000000000000000000000000001895963461786327125428610368722191005815672395983749858209273492131871103796478617527663555108805323353188177796328019430

? hurwitzzeta(50,1345.1234,128,10000)
Internal precision: # decimal digits: 228
Evaluation of the horizontal shift
horizontal shift = 75
Length of the EM-sum = 210
Raw horizontal shift determination time: 9 millisec
Raw precomputation time: 20 millisec
Raw EM comp time: 1885 millisec
Raw shift comp time: 976 millisec
Raw series comp time: 1290 millisec
Raw total time (no precomp, no horshift det): 4159 millisec
%31 = 0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000101940514555541454287276281790110006444
? zetahurwitz(50,1345.1234)
%25 = 0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001019405145555414542872762817895070526334567852612974058802912348663941543741814659131417041880139913006077614260624765457
? for(j=1,10000,zetahurwitz(50,1345.1234))
  ***   last result computed in 18,503 ms.

performances
 
 
? hurwitzzeta(3,1345.1234,128,10000)
Internal precision: # decimal digits: 66
Evaluation of the horizontal shift
horizontal shift = 25
Length of the EM-sum = 77
Raw horizontal shift determination time: 0 millisec
Raw precomputation time: 2 millisec
Raw EM comp time: 258 millisec
Raw shift comp time: 91 millisec
Raw series comp time: 265 millisec
Raw total time (no precomp, no horshift det): 618 millisec
%29 = 0.000000276546785405675961721446622621541608533156635172023384794409511879
 
? zetahurwitz(3,1345.1234) 
%46 = 0.000000276546785405675961721446622621541608533156635172023384794415969
? for(j=1,10000, zetahurwitz(3,1345.1234))
time = 4,467 ms.
? 4.467/.618
%30 = 7.22815533980582524271844660194174757281553398058252427184466019417
  

?  hurwitzzeta(10,1345.1234,128,10000)
Internal precision: # decimal digits: 91
Evaluation of the horizontal shift
horizontal shift = 34
Length of the EM-sum = 102
Raw horizontal shift determination time: 1 millisec
Raw precomputation time: 3 millisec
Raw EM comp time: 409 millisec
Raw shift comp time: 162 millisec
Raw series comp time: 348 millisec
Raw total time (no precomp, no horshift det): 924 millisec
%32 = 0.0000000000000000000000000000077330462226126189435010377485437297914952570333023366395536649285085852819161
? 

? 6.190/.924
%33 = 6.699134199134199134199134199134199134199134199134199134199134199134199134199134199134199134
? 


? for(j=1,10000, zetahurwitz(10,1345.1234))
time = 6,190 ms.
? zetahurwitz(10,1345.1234)
%50 = 0.00000000000000000000000000000773304622261261894350103774854372979149525703330233663953805731716806977404424017881121

 
 
? hurwitzzeta(20,1345.1234,128,10000)
Internal precision: # decimal digits: 125
Evaluation of the horizontal shift
horizontal shift = 45
Length of the EM-sum = 131
Raw horizontal shift determination time: 2 millisec
Raw precomputation time: 6 millisec
Raw EM comp time: 672 millisec
Raw shift comp time: 317 millisec
Raw series comp time: 548 millisec
Raw total time (no precomp, no horshift det): 1545 millisec
%35 = 0.0000000000000000000000000000000000000000000000000000000000001895963461786327125428610368722191005815672395983749858921

? 9.115/1.545
%38 = 5.8996763754045307443365695792880258899676375404530744336569579288025889967637540453074433656957928802588996763754045307443366

? for(j=1,10000, zetahurwitz(20,1345.1234))
time = 9,115 ms.
? zetahurwitz(20,1345.1234)
time = 1 ms.
%53 = 
0.0000000000000000000000000000000000000000000000000000000000001895963461786327125428610368722191005815672395983749858209273492131871103796478617527663555108805323353188177796328019430

 
? hurwitzzeta(8,1345.1234,128,10000)
Internal precision: # decimal digits: 84
Evaluation of the horizontal shift
horizontal shift = 32
Length of the EM-sum = 96
Raw horizontal shift determination time: 3 millisec
Raw precomputation time: 11 millisec
Raw EM comp time: 389 millisec
Raw shift comp time: 139 millisec
Raw series comp time: 345 millisec
Raw total time (no precomp, no horshift det): 880 millisec
%41 = 0.000000000000000000000017976152573162809213704220728153567202718996334681286706857092892145283352725

?  zetahurwitz(8,1345.1234)
%42 = 0.0000000000000000000000179761525731628092137042207281535672027189963346812867068570912087308024944892378327

? for(j=1,10000, zetahurwitz(8,1345.1234))
time = 5,480 ms. 
? 5.480/.880
%43 = 6.22727272727

? hurwitzzeta(8.3,1345.1234,128,10000)
Internal precision: # decimal digits: 84
Evaluation of the horizontal shift
horizontal shift = 32
Length of the EM-sum = 96
Raw horizontal shift determination time: 2 millisec
Raw precomputation time: 28 millisec
Raw EM comp time: 614 millisec
Raw shift comp time: 3389 millisec
Raw series comp time: 572 millisec
Raw total time (no precomp, no horshift det): 4583 millisec
%45 = 0.0000000000000000000000019855996153015416870238271027874537150218225654217185471133694617832247317290

?  zetahurwitz(8.3,1345.1234)
time = 14 ms.
%55 = 0.00000000000000000000000198559961530154168702382710278745371502182256542171854711336786501313993448235798

? for(j=1,10000, zetahurwitz(8.3,1345.1234))
time = 2min, 20,557 ms.

? 140.557/4.583
%46 = 30.6692123063495526947414357407811477198341697578005673139864717433995199650883700633

---------------------------------------------------
Dirichlet beta

? v=[32,64,80,128]
%1 = [32, 64, 80, 128] 

? for (i=1,4, print("----- "); print("binary precision = ",v[i]); print(dirichletbeta(2,v[i])-verifdirichletbeta(2)))
----- 
binary precision = 32
Internal precision: # decimal digits: 31
Raw precomputation time: 0 millisec
Raw series computation time: 0 millisec
Raw computation time: 0 millisec
----- 
0.0000000000000000000000000000000
----- 
binary precision = 64
Internal precision: # decimal digits: 41
Raw precomputation time: 0 millisec
Raw series computation time: 0 millisec
Raw computation time: 0 millisec
----- 
0.00000000000000000000000000000000000000000
----- 
binary precision = 80
Internal precision: # decimal digits: 46
Raw precomputation time: 1 millisec
Raw series computation time: 0 millisec
Raw computation time: 1 millisec
----- 
0.0000000000000000000000000000000000000000000000
----- 
binary precision = 128
Internal precision: # decimal digits: 60
Raw precomputation time: 0 millisec
Raw series computation time: 0 millisec
Raw computation time: 0 millisec
----- 
0.000000000000000000000000000000000000000000000000000000000000000000000000000008636168555094444625

? dirichletbeta(2,128,10000) 
Internal precision: # decimal digits: 60
Raw precomputation time: 1 millisec
Raw series computation time: 244 millisec
Raw computation time: 245 millisec
----- 
%1 = 0.915965594177219015054603514932384110774149374281672134266498

 
? for(i=1,10000,verifdirichletbeta(2))
time = 882 ms.
  

? for (i=1,4, print("----- "); print("binary precision = ",v[i]); print(dirichletbeta(8,v[i])-verifdirichletbeta(8)))
----- 
binary precision = 32
Internal precision: # decimal digits: 40
Raw precomputation time: 1 millisec
Raw series computation time: 0 millisec
Raw computation time: 1 millisec
----- 
-0.00000000000000000000000000000000000000000000000000000000015930919111324522770
----- 
binary precision = 64
Internal precision: # decimal digits: 50
Raw precomputation time: 0 millisec
Raw series computation time: 0 millisec
Raw computation time: 0 millisec
----- 
-0.00000000000000000000000000000000000000000000000000000000015930919111324522770
----- 
binary precision = 80
Internal precision: # decimal digits: 55
Raw precomputation time: 1 millisec
Raw series computation time: 0 millisec
Raw computation time: 1 millisec
----- 
-0.00000000000000000000000000000000000000000000000000000000015930919111324522770
----- 
binary precision = 128
Internal precision: # decimal digits: 69
Raw precomputation time: 0 millisec
Raw series computation time: 0 millisec
Raw computation time: 0 millisec
----- 
0.000000000000000000000000000000000000000000000000000000000000000000000
time = 2 ms.
 
 

? for(i=1,10000,verifdirichletbeta(8))
time = 961 ms.
? verifdirichletbeta(8)
%5 = 0.99984999024682965633806705924046378147600743300742806972498743


? dirichletbeta(8,128,10000) 
Internal precision: # decimal digits: 69
Raw precomputation time: 1 millisec
Raw series computation time: 333 millisec
Raw computation time: 334 millisec
----- 
time = 334 ms.
%21 = 0.999849990246829656338067059240463781476007433007428069724987429240671


? for (i=1,4, print("----- "); print("binary precision = ",v[i]); print(dirichletbeta(8.3,v[i])-verifdirichletbeta(8.3)))
----- 
binary precision = 32
Internal precision: # decimal digits: 40
Raw precomputation time: 7 millisec
Raw series computation time: 0 millisec
Raw computation time: 7 millisec
----- 
-0.00000000000000000000000000000000000000000000000000000000015930919111324522770
----- 
binary precision = 64
Internal precision: # decimal digits: 50
Raw precomputation time: 5 millisec
Raw series computation time: 0 millisec
Raw computation time: 5 millisec
----- 
-0.00000000000000000000000000000000000000000000000000000000015930919111324522770
----- 
binary precision = 80
Internal precision: # decimal digits: 55
Raw precomputation time: 4 millisec
Raw series computation time: 0 millisec
Raw computation time: 4 millisec
----- 
-0.00000000000000000000000000000000000000000000000000000000015930919111324522770
----- 
binary precision = 128
Internal precision: # decimal digits: 69
Raw precomputation time: 5 millisec
Raw series computation time: 0 millisec
Raw computation time: 5 millisec
----- 
-0.00000000000000000000000000000000000000000000000000000000000000015767972126762518684


? dirichletbeta(8.3,128,10000) 
Internal precision: # decimal digits: 69
Raw precomputation time: 8 millisec
Raw series computation time: 535 millisec
Raw computation time: 543 millisec
----- 
time = 543 ms.
%17 = 0.999891872076192554837935987435697213602603483733314431452489852735161

? for(i=1,10000,verifdirichletbeta(8.3))
time = 7,095 ms.
? verifdirichletbeta(8.3)
time = 1 ms.
%70 = 0.99989187207619255483793598743569721360260348373331443145248985
? 7.095/.535
%18 = 13.2616822429906542056074766355140186915887850467289719626168224299065
? 



? for (i=1,4, print("----- "); print("binary precision = ",v[i]); print(dirichletbeta(2^(7),v[i])-verifdirichletbeta(2^(7))))
----- 
binary precision = 32
Internal precision: # decimal digits: 224
Raw precomputation time: 0 millisec
Raw series computation time: 1 millisec
Raw computation time: 1 millisec
----- 
-0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000006441148769597133308
----- 
binary precision = 64
Internal precision: # decimal digits: 234
Raw precomputation time: 1 millisec
Raw series computation time: 0 millisec
Raw computation time: 1 millisec
----- 
0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
----- 
binary precision = 80
Internal precision: # decimal digits: 239
Raw precomputation time: 2 millisec
Raw series computation time: 0 millisec
Raw computation time: 2 millisec
----- 
0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
----- 
binary precision = 128
Internal precision: # decimal digits: 253
Raw precomputation time: 1 millisec
Raw series computation time: 1 millisec
Raw computation time: 2 millisec
----- 
0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000

? dirichletbeta(2^7,128,10000) 
Internal precision: # decimal digits: 253
Raw precomputation time: 1 millisec
Raw series computation time: 1307 millisec
Raw computation time: 1308 millisec
----- 
time = 1,308 ms.
%53 = 0.9999999999999999999999999999999999999999999999999999999999999151836857678944802122803476734813190482512090889661239183039835241779393226085986016801247031937392944645724927288136731327812710875032577528729926045110192032286586045639850796011060503298666

? for(i=1,10000,verifdirichletbeta(2^7))
time = 5,477 ms.
-------------------------------------------------
Dirichlet betaprime

? dirichletbetaprime(2,128,10000)
Internal precision: # decimal digits: 58
Raw precomputation time: 125 millisec
Raw series computation time: 229 millisec
Raw computation time: 354 millisec
----- 
%56 = 0.08158073611659279510291216978594115145773887519885962618374

? 


? for(i=1,10000,verifdirichletbetaprime(2))
time = 17,407 ms.


? verifdirichletbetaprime(2)
time = 2 ms.
%76 = 0.08158073611659279510291216978594115145773887519885962618374

? dirichletbetaprime(8.3,128,10000) 
Internal precision: # decimal digits: 62
Raw precomputation time: 160 millisec
Raw series computation time: 443 millisec
Raw computation time: 603 millisec
----- 
%58 = 0.00011805486611657653425932052800811664636460897881682366533768402

? verifdirichletbetaprime(8.3)
time = 2 ms.
%20 = 0.00011805486611657653425932052800811664636460897881682366533768402
? for(i=1,10000,verifdirichletbetaprime(8.3))
time = 23,956 ms.

? 23.956/0.603
%59 = 39.728026533996683250414593698175787728026533996683250414593698

? dirichletbetaprime(2^7,128,10000) 
Internal precision: # decimal digits: 154
Raw precomputation time: 1247 millisec
Raw series computation time: 613 millisec
Raw computation time: 1860 millisec
----- 
%61 = 0.00000000000000000000000000000000000000000000000000000000000009318024509492700967831012686159465914022707331503592551218816965162041248130174447467473438297578566735117170435373179718204044154126896220036056162697301
? verifdirichletbetaprime(2^7)
time = 3 ms.
%14 = 0.0000000000000000000000000000000000000000000000000000000000000931802450949270096783101268615946591402270733150359255121881696516204124813017444746747343766749
? for(i=1,10000,verifdirichletbetaprime(2^7))
time = 39,447 ms.

-------------------------------------------------
Dirichlet beta log derivative

? dirichletbetavalues(2,128,10000) 
Internal precision: # decimal digits: 58
Raw precomputation time: 124 millisec
Raw series computation time: 286 millisec
Raw computation time: 410 millisec
----- 
beta(s) = 0.9159655941772190150546035149323841107741493742816721342665
betaprime(s) = 0.08158073611659279510291216978594115145773887519885962618374
betaprime/beta(s) = 0.08906528436788503775577121532873506418642022280689761350881
%2 = 0.08906528436788503775577121532873506418642022280689761350881

? verifdirichletbetalogder(2) 
%26 = 0.08906528436788503775577121532873506418642022280689761350881
? for(i=1,10000,verifdirichletbetalogder(2))
time = 19,097 ms.

? dirichletbetavalues(2^7,128,10000)
Internal precision: # decimal digits: 154
Raw precomputation time: 1269 millisec
Raw series computation time: 802 millisec
Raw computation time: 2071 millisec
----- 
beta(s) = 0.9999999999999999999999999999999999999999999999999999999999999151836857678944802122803476734813190482512090889661239183039835241779393226085986016801247032
betaprime(s) = 0.0000000000000000000000000000000000000000000000000000000000000931802450949270096783101268615946591402270733150359255121881696516204124813017444746747343828790
betaprime/beta(s) = 0.0000000000000000000000000000000000000000000000000000000000000931802450949270096783101268615946591402270733150359255121881775548253606772400261364135620336229
%4 = 0.0000000000000000000000000000000000000000000000000000000000000931802450949270096783101268615946591402270733150359255121881775548253606772400261364135620336229

? for(i=1,10000,verifdirichletbetalogder(2^7))
time = 44,461 ms.

? verifdirichletbetalogder(2^7)
time = 4 ms.
%30 = 0.0000000000000000000000000000000000000000000000000000000000000931802450949270096783101268615946591402270733150359255121881775548253606772400261364135620273579

? dirichletbetavalues(8.3,128,10000) 
Internal precision: # decimal digits: 62
Raw precomputation time: 156 millisec
Raw series computation time: 301 millisec
Raw computation time: 457 millisec
----- 
beta(s) = 0.99989187207619255483793598743569721360260348373331443145248985
betaprime(s) = 0.00011805486611657653425932052800811664636460897881682366533768402
betaprime/beta(s) = 0.00011806763252455027426566425441126211895938523953494237401120170
%9 = 0.00011806763252455027426566425441126211895938523953494237401120170


?  for(i=1,10000,verifdirichletbetalogder(8.3))
time = 23,711 ms.

? verifdirichletbetalogder(8.3)
%23 = 0.00011806763252455027426566425441126211895938523953494237401120170

? 23.711/.457
%10 = 51.884026258205689277899343544857768052516411378555798687089716
-------------------------------------------------
-------------------------------------------------

Catalan constant 
 
? my_catalan(10)
Catalan constant computation
Catalan constant computation time: 0 min, 0 sec, 0 millisec
%1 = 0.91596559417721901505

? my_catalan(100)
Catalan constant computation
Catalan constant computation time: 0 min, 0 sec, 0 millisec
%2 = 0.91596559417721901505460351493238411077414937428167213426649811962176301977625476947935651292611510624857442262


? my_catalan(1000)
Catalan constant computation
Catalan constant computation time: 0 min, 0 sec, 230 millisec
%3 = 0.91596559417721901505460351493238411077414937428167213426649811962176301977625476947935651292611510624857442261919619957903589880332585905943159473748115840699533202877331946051903872747816408786590902470648415216300022872764094238825995774150881639747025248201156070764488380787337048990086477511322599713434074854075532307685653357680958352602193823239508007206803557610482357339423191498298361899770690364041808621794110191753274314997823397610551224779530324875371878665828082360570225594194818097535097113157126158042427236364398500173828759779765306837009298087388749561089365977194096872684444166804621624339864838916280448281506273022742073884311722182721904722558705319086857354234985394983099191159673884645086151524996242370437451777372351775440708538464401321748392999947572446199754961975870640074748707014909376788730458699798606448749746438720623851371239273630499850353922392878797906336440323547845358519277777872709060830319943013323167124761587097924554791190921262018548039639342434956537597

?  my_catalan(2000)
Catalan constant computation
Catalan constant computation time: 0 min, 1 sec, 834 millisec
%4 = 0.915965594177219015054603514932384110774149374281672134266498119621763019776254769479356512926115106248574422619196199579035898803325859059431594737481158406995332028773319460519038727478164087865909024706484152163000228727640942388259957741508816397470252482011560707644883807873370489900864775113225997134340748540755323076856533576809583526021938232395080072068035576104823573394231914982983618997706903640418086217941101917532743149978233976105512247795303248753718786658280823605702255941948180975350971131571261580424272363643985001738287597797653068370092980873887495610893659771940968726844441668046216243398648389162804482815062730227420738843117221827219047225587053190868573542349853949830991911596738846450861515249962423704374517773723517754407085384644013217483929999475724461997549619758706400747487070149093767887304586997986064487497464387206238513712392736304998503539223928787979063364403235478453585192777778727090608303199430133231671247615870979245547911909212620185480396393424349565375967394943547300143851807050512507488613285641293449595022987229831628948164616225739894762318195420066071881427594975599589836373037675338533813545031276817240118140721534688316835681686393272936775866739258395406180333878306870649014334860172981069921799565309581871579115539560366890369904939667538437758104931899553855162621962533168040162737521301209406045387950760538271231974679008823691786155733891244172238339381481207759942984917243976685756327180688082799829793788494327249346576074905438748195268130744370462946358928102765317050765479744948399489594770927885911958487241278660840885545978238124922605056100945844866989585768716111717866623368474099493855413210937552818155258815915022282444544417186099465881517664960782236789705192697113125713754543701243296730572468450158193130160877662156509575546796667866170823476825581335186819377456500145652617040960746889539302347919806000842455621751084234717363878793695778784409337922198945753409616474245546224787880029229148036907115270795546

? my_catalan(3000)
Catalan constant computation
Catalan constant computation time: 0 min, 7 sec, 265 millisec
%1 = 0.9159655941772190150546035149323841107741493742816721342664981196217630197762547694793565129261151062485744226191961995790358988033258590594315947374811584069953320287733194605190387274781640878659090247064841521630002287276409423882599577415088163974702524820115607076448838078733704899008647751132259971343407485407553230768565335768095835260219382323950800720680355761048235733942319149829836189977069036404180862179411019175327431499782339761055122477953032487537187866582808236057022559419481809753509711315712615804242723636439850017382875977976530683700929808738874956108936597719409687268444416680462162433986483891628044828150627302274207388431172218272190472255870531908685735423498539498309919115967388464508615152499624237043745177737235177544070853846440132174839299994757244619975496197587064007474870701490937678873045869979860644874974643872062385137123927363049985035392239287879790633644032354784535851927777787270906083031994301332316712476158709792455479119092126201854803963934243495653759673949435473001438518070505125074886132856412934495950229872298316289481646162257398947623181954200660718814275949755995898363730376753385338135450312768172401181407215346883168356816863932729367758667392583954061803338783068706490143348601729810699217995653095818715791155395603668903699049396675384377581049318995538551626219625331680401627375213012094060453879507605382712319746790088236917861557338912441722383393814812077599429849172439766857563271806880827998297937884943272493465760749054387481952681307443704629463589281027653170507654797449483994895947709278859119584872412786608408855459782381249226050561009458448669895857687161117178666233684740994938554132109375528181552588159150222824445444171860994658815176649607822367897051926971131257137545437012432967305724684501581931301608776621565095755467966678661708234768255813351868193774565001456526170409607468895393023479198060008424556217510842347173638787936957787844093379221989457534096164742455462247878800292291480369071152707955455054147826884981852460058144665178681423154114878554099665167385397276146970169043915114900893330791845746576209967754812313820154360109885272162977010876157478173564163698570355340672649351963169554767211507772315900448338260516116383430865139797225161741385381293248011946362518800840398194553905518210424606292185217560246548601929767239740511039526456924297864212424037518926787296027177337873837997832667620861195206791215126382119252329404069205994386427469321533885667117330827142408332659203260753165928042310230997358400395940342632227688070118681961767809056315815978453763757835637359027716488313102887693795053507320801807581022382308031762504329424722268391229712955351355104314761886655474367692184120188771617992285620563522054703200691808688066121174204060992412348760515406820226255950481248589411873583468229042308361555476947777083194087481249167489290065936961641662343683707543963838945144011955648738134292122982001302107996192242492449305199923585815808260352497998505918669722
? ##
  ***   last result computed in 7,270 ms.
  
  ? my_catalan(10000)
Catalan constant computation
Catalan constant computation time: 7 min, 25 sec, 982 millisec
%2 = 0.91596559417721901505460351493238411077414937428167213426649811962176301977625476947935651292611510624857442261919619957903589880332585905943159473748115840699533202877331946051903872747816408786590902470648415216300022872764094238825995774150881639747025248201156070764488380787337048990086477511322599713434074854075532307685653357680958352602193823239508007206803557610482357339423191498298361899770690364041808621794110191753274314997823397610551224779530324875371878665828082360570225594194818097535097113157126158042427236364398500173828759779765306837009298087388749561089365977194096872684444166804621624339864838916280448281506273022742073884311722182721904722558705319086857354234985394983099191159673884645086151524996242370437451777372351775440708538464401321748392999947572446199754961975870640074748707014909376788730458699798606448749746438720623851371239273630499850353922392878797906336440323547845358519277777872709060830319943013323167124761587097924554791190921262018548039639342434956537596739494354730014385180705051250748861328564129344959502298722983162894816461622573989476231819542006607188142759497559958983637303767533853381354503127681724011814072153468831683568168639327293677586673925839540618033387830687064901433486017298106992179956530958187157911553956036689036990493966753843775810493189955385516262196253316804016273752130120940604538795076053827123197467900882369178615573389124417223833938148120775994298491724397668575632718068808279982979378849432724934657607490543874819526813074437046294635892810276531705076547974494839948959477092788591195848724127866084088554597823812492260505610094584486698958576871611171786662336847409949385541321093755281815525881591502228244454441718609946588151766496078223678970519269711312571375454370124329673057246845015819313016087766215650957554679666786617082347682558133518681937745650014565261704096074688953930234791980600084245562175108423471736387879369577878440933792219894575340961647424554622478788002922914803690711527079554550541478268849818524600581446651786814231541148785540996651673853972761469701690439151149008933307918457465762099677548123138201543601098852721629770108761574781735641636985703553406726493519631695547672115077723159004483382605161163834308651397972251617413853812932480119463625188008403981945539055182104246062921852175602465486019297672397405110395264569242978642124240375189267872960271773378738379978326676208611952067912151263821192523294040692059943864274693215338856671173308271424083326592032607531659280423102309973584003959403426322276880701186819617678090563158159784537637578356373590277164883131028876937950535073208018075810223823080317625043294247222683912297129553513551043147618866554743676921841201887716179922856205635220547032006918086880661211742040609924123487605154068202262559504812485894118735834682290423083615554769477770831940874812491674892900659369616416623436837075439638389451440119556487381342921229820013021079961922424924493051999235858158082603524979985059186697220123164897104830701793528112228966355128317437352393011402792389808744569648309013207877658785362301354280001629055877295006795876178247374871378060042208445346045064702443258085164777173903196028655538328281415915248735263307150513147882844999238663243198106336515243311321463900933362159160744482923457177454817169580181688900175285645046489139090420356029836045652425265797270138586757653899302958449258692189788644388819358114526770563160609737684654083694230203816826392458579107404870879877852426140868715178575801006023681703491797733622196629537718913853116739965565885912164628015582629873541376336076073020045591202946657347571852745311633847776486838248504116301605227086944442703644251242363971814999234960838959168258036164749881042639483890042940550431502193126864230059992926361540649262664186583594904249371523622068403940370108680740098440001512465343535067233845469463576021186762114341424761178341043127306116782248833969915539091310973231066781117485537679027231845076545775699887411395686146631581361573674061881125914620397423401125882131569075175754979658229689846231329257273175338302313533232870056595688534175204573932758183513982347678009261426521074710456668763134325667275929891952548849037809046546488268575204454695053813498309021460489718319387780863409014168285484524248093104343217724788778248739486061800233415225914146138782700545170971410457656614928953108672486080484204376637936230213645817798022720882738071736711299822289069125763027779162651035762577038104288680376054636303337940367377696744757171918712803954370966413877226626889837311111602004518593973174764621542838460162144526553720292552051504941828003032550267579038252786139633572720650890367820176258573636602459644914533528141037251683822090097101943680278336708963314672497329503919259298514966414498521873384370124517467421871213110205726174340134056876555104187866544518902765005382178609412105353899784905982180023067890821606141367018393687028304544346780536499566495053180837980207950365835227622006506786177171095672005629703023553593373869771832835337557262344415664916005762666604199085276789703504193295554568745338842121304879862000928706178007678592735175386523677348535053066125396025536280809350562562821347432394399222442739711562755985244339104126180433506987134104280978456869518977668882650503756167591535473173668135683353168588440266726203196600785194905261819016135540883210564405409027216204498851041761292787884227851835200704439460961571665543448392802592501156306227650740050312351417656526449943042570531502230552233576634208943102385867060630430297719853224212043298619528633162199479803021651170071853216768095061934167286284674753307211005511854225758629292681406381602461376952043278677852351940897487799588262651018857167526448964259516245608164680586662605844328281537669209501700131691093864391470033345906701868799246483109181855848104631118954767258303668922657116990565431759988680286731145873457549777440562265841337924742718870078268554656782290336462515389898481303382848801578806646984480217166949381713998561787717378778712739969834267499971632268327257972572115422822471585175485105077970961560718370771383998265316365376758124751878398350457588311790755545686617395928729558719386219223573876438607017401059359744278581411271395680504961269960048434583896436697014771140329178065084925873008209906179587580402996618292021820696155745628810980223576195163967867626609736795492343789154100185728989816837858427303612448456532426353483148925506480782198270518366562137380923695907752151698346526103237738415089830658136487130918231383360055922540017526278742124582625282370841549068233176525686246245609564332012497970680412465220417099693819728527361263918229564824346904280303582683293573927934144962552827643618433542620664134683156370222632876839087900597166332580664331095881812753248627892980094868158902452714692410818394313034916873369765811519402277339800954019992521514349607341474539039230411990899640390760329165111929551028666741487888146370780055212485563601811272094261309853801454061531585422664625843416142595014823689366366735542832720063760749018108182214340861973911548328544384315811917349722801704172459572971660809528522210471512385830056016372167813179800095725635672059859413601259647770490028817068063068944393809042066274116418284908815132355567686236301149153615835301929518669408266688060717443297219890870195003345427130936321480897947531976335087806006735145801884222411475843568459595642012392746288926431857096931750781433192768351309119874941923765705532160176241226139675957034480490314072757977766284315056551277393854632895972421929264699107855832088971233051922897717260248197053719683239188010267755856584529891014031105506836583574883237481454355676661805833306498942393994388379965296325401019799634146428782283377256296616846911720121828661474404077735591989241075123126002042456108959229931398206017121311234380896737575014319013773810588142008068583268348860835919739558648453632085482935260075146055373154056807915110104867854204502431784812177640696641654466861923921035823089944277756117355726622412507724360772601174557283433060318105899040780191874490094925401279748286924188618885357723618971933738259474905697075534520115315402733311828377976124717015057625706465529679989718223087381422674357748520494951554930115312237349799758257618202831754452121431136790914270888108434238961583657657348848869901897285510569345637556197136235680510433375824604783740671629050969122053094969342273645888616180749732254029065265157374208441696502868219289687780182022916709817519263029166075560544631057607509820784122037469094403481232669125215712433087010353789034900090664329264372577975024182433028437809210844805131449637450450720821172885237993949048626138199252205239306727369305935217637216618890419427844146908993764901848571379429452377041302501047569886854093043365370302936769520956382771770453471932906996474983798241227697609433107940858403105491213134469627575233133732110803099242578856580410400798316619455697015999176218878664471943175469916497415202518064488778879851669666905018069611778953196777655119951365783552962285756790257213860335416088979846955972160571371734597829251501311460097902373072412036354335883448781036365025215337515426356641550685916208348770008127064546983842306741350697392385488411364274333848991803820765123559375590563012417219171980889432932113650166895367035795783718714129279518231967750770527301435566890987728156213938442321621621619714293261707601370746636004397668423210772544037792407419121770198942861545330158425131835182081137819138154235421559608389967085629516895262671993966405734193091489098184941482188411810178142374557350988488825815942045631984223214360317409340713561767809014347986994190579589435555830750700972883787924906533440124683055020087050873645653095557917071622999437962377763190566337687723700887467453957110781672472188004312957875214467796449960511800034246725245859539203204670880889789317541179493224981308995641838295378715037761463626
? ##
  ***   last result computed in 7min, 25,984 ms.


*******/