Jacques Gélinas on Wed, 28 Nov 2018 19:06:02 +0100

 Learning with GP: Taylor coefficients of xi(s) - Kreminski

• To: "pari-users@pari.math.u-bordeaux.fr" <pari-users@pari.math.u-bordeaux.fr>
• Subject: Learning with GP: Taylor coefficients of xi(s) - Kreminski
• From: Jacques Gélinas <jacquesg00@hotmail.com>
• Date: Wed, 28 Nov 2018 18:05:55 +0000
• Accept-language: fr-CA, en-US
• Delivery-date: Wed, 28 Nov 2018 19:06:02 +0100
• Thread-topic: Learning with GP: Taylor coefficients of xi(s) - Kreminski

```Accurate values of derivatives of xi(s) at s=1/2 and at s=1 are
available in [8]; their precision decreases almost linearly from
3276 digits for xi^(2)(1/2) down to 328 digits for xi^(1000)(1/2).
The format used includes what can be interpreted as an embeded
Mathematica precision estimate, and we run a GNU sed script to
convert the two Mathematica data files for PARI/GP.

##------------------------------------------- Start cut

cat > xid.sed << 'EOF'
#============================================= xid.sed
#
#	Convert a data file from Mathematica to PARI/GP
#
#			xiderivativeatsequalhalf[1] = 0
/ = 0\$/d;
#			blank lines
/^.\?\$/d;
#			2864824522574683538316274864488477`121\
#			.4772*^75
#		==>	2864824522574683538316274864488477`121.4772*^75
#
#			2864824522574683538316274864488477`121\
#			.069
#		==>	2864824522574683538316274864488477`121.069
/`.*\\\$/h;
/`.*\\\$/d;
/`\|\\\$/!H;
/`\|\\\$/!x;
s/\\\n//;
#			688689767666203002349477596872`77.6054*^80
#		==>	688689767666203002349477596872E80`77.6054
s/`\(.*\)\*^\(.*\)/E\2`\1/;
#			14693726641604350416778354683116048`800
#		==>	14693726641604350416778354683116048; \\800
s/`/; \\\\/;

#============================================= end-of-file
EOF

echo \
wget http://faculty.csupueblo.edu/rick.kreminski\
/stieltjesrelated/xiderivativeatsequalhalf.txt
sed -f xid.sed xiderivativeatsequalhalf.txt > xid12.gp
echo \
wget http://faculty.csupueblo.edu/rick.kreminski\
/stieltjesrelated/xiderivativeatsequalone.txt
sed -f xid.sed xiderivativeatsequalone.txt > xid1.gp

##-------------------------------------------- End paste

With the "xid12.gp" and "xid1.gp" data files, one can plot
the derivatives and their decimal precision with GP.

\\------------------------------------------- Start cut

X12 = 2*[1..500];
Y12 = vector(1000);
alias(xiderivativeatsequalhalf,Y12);
Y12 = [Y12[k] | k<-X12];
print("log2 plot of derivatives of xi(s) at s=1/2");
plothraw(X12,apply(exponent,Y12),1)
print("precision of derivatives of xi(s) at s=1/2");
plothraw(X12,apply(precision,Y12),1)

X1 = [1..1400];
Y1 = vector(1400);
alias(xiderivatone,Y1);
print("log2 plot of derivatives of xi(s) at s=1");
plothraw(X1,apply(exponent,Y1),1)
print("precision of derivatives of xi(s) at s=1");
plothraw(X1,apply(precision,Y1),1)

\\-------------------------------------------- End paste

Kreminski's indirect method [9] uses series tranformations starting
with extremely accurate values of the "Stieltjes constants" from
the Laurent series of zeta(s) at s=1 to obtain, firstly a series for
D[log(xi(s))] centered at s=1, secondly at s=1/2, thirdly for
log(xi(s)) centered at s=1/2 by integration, and finally
for xi(s) at s=1/2 by exponentiation.

All the series operations are immediately available in GP, but
the Kreminski Stieltjes constants with more than ten thousand
correct digits computed by Newton-Cotes of very high order [10]
have not all been made available [11].

We will simply compare the resulting transformed series to
the ones obtained directly by the expansion of log(xi(s))
and of xi(s) which shows a disturbing loss of precision
in Kreminski's indirect method.

\\------------------------------------------- Start cut

\\ floating-point `equality' based on suggestion of Karim Belabas
fleq(x,y,r=7/8) = if(x==y,1,exponent(normlp(x-y))/exponent(0.) > r);
exr(x,y) = exponent( if(!x,y,1-y/x) );
IS(e,msg) = print("Test ",msg,if(e, " passed"," failed"));
\\ even and odd part of series
sev(S) = my(p=Vec(Pol(S)),v=variable(S)); sum(k=0,#p\2-1,p[#p-2*k]*v^(2*k));
sod(S) = my(p=Vec(Pol(S)),v=variable(S)); sum(k=0,#p\2-1,p[#p-2*k-1]*v^(2*k+1));

default(seriesprecision,32);
default(realbitprecision,192);

xis(s) = Pi^(-s/2)*gamma(1+s/2)*(s-1)*zeta(s);
lxi(s) = -s/2*log(Pi) + log(gamma(1+s/2)) + log((s-1)*zeta(s));
IS( fleq(lxi(1/2),log(xis(1/2))), "lxi = log(xi)" );

SZ  = z*zeta(1+z); \\ entire function=1-Euler*z+0.0728*z^2-...
S1T = Pol(log(xis(1+z))',z);
S1  = Pol(-1/2*log(Pi) + 1/2*psi(3/2+z/2) + SZ'/SZ,z);
IS( fleq(S1,S1T), "S1 = S1T" );

\\ half of precision is lost after this (partial) `shift'
S2T = sod(log(xis(1/2+x))');
S2  = sod(subst(S1,z,x-1/2));
IS( fleq(S2,S2T,1/2), "S2 ~= S2T (1/2 prec)" );

S3T = sev(log(xis(1/2+x)));
IS( fleq(S3T',S2T), "S3T' = S2T" );
S3 = intformal(S2) + log(xis(1/2));
IS( fleq(S3,S3T,1/2), "S3 ~= S3T (1/2 prec)" );

S4T = sev(xis(1/2+x)); \\ Target: 0.49712 + 0.011486*x^2 + ...
IS( fleq(log(S4T),S3T), "log(S4T) = S3T" );
S4  = sev(exp(S3));
IS( fleq(S4,S4T,1/2), "S4 ~= S4T (1/2 prec)" );

\\ relative errors using Kreminski's trusted values
C4T = Vec(polrecip(S4T));
IS( [ exr( Y12[k]/(2*k)!, C4T[2*k+1] ) | k<-[1..10] ] ==
[-186, -175, -166, -158, -146, -134, -124, -114, -101, -91],\
"of errors in xis(1/2+x) expansion by PARI/GP (S4T)");

C4 = Vec(polrecip(S4));
IS( [ exr( Y12[k]/(2*k)!, C4[2*k+1] ) | k<-[1..10] ] ==
[-138, -123, -109, -96, -84, -72, -60, -49, -38, -28],\
"of errors in xis(1/2+x) expansion by series transformations (S4)");

\\-------------------------------------------- End paste

------------------------------ Features used [3,4]
alias, plothraw, seriesprecision, relbitprecision, psi,
subst, intformal, log(series), exp(series), series/series,
normlp, exponent, variable, log, exp, gamma, Euler, zeta,
polrecip, Vec, Pol(note that the default variable is x)

------------------------------ References
[8] Rick Kreminski (2004,2005) Computations of Taylor coefficient values for xi(s) at s=1/2.
http://faculty.csupueblo.edu/rick.kreminski/stieltjesrelated
[9] Rick Kreminski (2004) Unpublished work method
http://faculty.csupueblo.edu/rick.kreminski/stieltjesrelated/xiweb
[10] Kreminski (2003) Newton-Cotes integration for approximating Stieltjes (generalized Euler) constants
Mathematics of computation 72, 1379-1397
http://www.ams.org/mcom/2003-72-243/S0025-5718-02-01483-7/home.html
[11] Rick  Kreminski (2004) Stieltjes constant values to 6900 digits, gamma_1 to gamma_12
http://faculty.csupueblo.edu/rick.kreminski/stieltjesrelated/gammas1to12/gammas1to12to6900digits/

Jacques Gélinas

```