| Jacques Gélinas on Wed, 28 Nov 2018 19:06:02 +0100 |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
| 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);
read("xid12.gp"); \\ Kreminski (2005)
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);
read("xid1.gp"); \\ Kreminski (2005)
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