Learning with GP: Taylor coefficients of xi(s) - Lyness & Moler

[Jacques Gélinas <jacquesg00@hotmail.com>]

Using Pari/GP, we will plot the interpolation error as
the number of nodes chosen by Gram in 1895 is increased.
We must first accelerate the interpolation computations.

Gram proceeded by stages. If log Xi(t) = a + b*t^2 + ...
and a is known, then (log(Xi(t)) - a)/t^2 = b + c*t^2 + ...
where b can be found by interpolation, and he repeated
this process to evaluate the other parameters. He only
needed the value at t=0 of each polynomial, which can
be computed efficiently by Neville's algorithm [5]. 
We use a more general p-coefficients method, adapted
from Algol procedures based on Vandermonde systems [6].

Both Gram and [6] suggest the transformation x=t^2
for even functions to halve the degree of interpolation:
if g(x)=f(t^2), then f^(2k)(0)/(2k)! == g^(k)(0)/k! by Taylor.

We finally plot the interpolation error for Xi(t) in fonction
of the number of extra nodes used. It has the allure
of a semi-convergent series: the infinite Newton's
series probably diverges for Xi(t) which has too many zeros
for the growth condition of Carlson's theorem [7].

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

/* Lowest p Coefficients of Interpolating Polynomial 
   for f(x) based on X[1,...,N], N=p+m, m>=0  */

cip(f,X=[0],p=#X,m=0) = {
  local( N=p+m, C=vector(Cdim(N)) );
  if(#X<N, error("cip: ",N-#X," nodes missing"));
  for(k=1, N, update(k,p,X,f(X[k])) );  
  C[N+1-p..N]; \\ constant term last for Pol()
}
Cdim(k=1) = 1 + (k-1)*(k+2)/2;

/*  Neville's interpolation algorithm */
update(k,p,X,fk) = {
  local( xk=X[k], m=Cdim(k), n, xkd );
  C[m] = fk;
  for(d=2, k, xkd = xk - X[k-d+1];
   for(s=1, min(p,d), n = m + d - 2;
     C[m--] = if(s==1, C[n]+(xk*(C[n-k] - C[n]))/xkd,
              if(s==d, (C[n+1] - C[n-k+1])/xkd,
       C[n]+(xk*(C[n-k] - C[n])+C[n+1] - C[n-k+1])/xkd )));
   if(d>p, m -= d-p ));}

/* Lucky seven test with even Legendre polynomials */
pt(n,t,j) = pollegendre(2*n,t);  \\ even P(t) target 
tv(N=15)  = concat(0,[k-1/2|k<-[1..N]]); \\ Gram nodes
xv(N=15)  = apply(sqr, tv(N));       \\ (x=t^2)-nodes
px(n,x)   = pt(n, sqrtint(4*x)/2);   \\ Q(x)=P(sqrt(x))
cpxi(n,m) = cip(x->px(n,x),xv(n+m)); \\ m extra nodes
pti(n,t,m) = subst(Pol(cpxi(n,m),x), x, t^2);
print("Legendre test pti = pt ", if( 7*3 == sum(n=1,7,\
 sum(m=0,2, pti(n,t,m)==pt(n,t) )), "passed", "failed" ));

/* Error plots */
Xi(t)  = if(t^2==-1/4, 1/2, xis(1/2-I*t));
xis(s) = Pi^(-s/2)*gamma(1+s/2)*(s-1)*zeta(s);

ckn(k)   = derivnum(t=0,Xi(I*t),k)/k!; \\ used as exact value
cki(k,m) = cip(x->Xi(I*sqrtint(4*x)/2),xv(k/2+m),k/2+1,m)[1];

rer(x,y) = exponent( if(!x, y, 1-y/x) );   \\ Karim Belabas
M = 16*[0..16];         \\ numbers of extra nodes for plot
ple(k) = my( ck = ckn(k) );\
     plothraw(M, [ rer(ck,cki(k,m))|m<-M ], 1);

print("High-precision plots for k=4,20,36");
print("(coefficient error exponent/number of extra Gram nodes)");

localbitprec(640); forstep(k=4,36,16, ple(k));

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

------------------------------ Features used [3,4]
derivnum, exponent, for, forstep, localbitprec, min, Pol,
plothraw, pollegendre, sqr, sqrtint, subst, sum, (global C)

------------------------------ References
[5] Wikipedia (2018) Neville's algorithm
    http://en.wikipedia.org/wiki/Neville's_algorithm
[6] Lyness & Moler (1966) Van Der Monde systems and numerical differentiation
    http://www.digizeitschriften.de/dms/img/?PID=GDZPPN001166115
[7] Wikipedia (2018) Carlson's_theorem
    http://en.wikipedia.org/wiki/Carlson's_theorem


Jacques Gélinas          Next: Taylor coefficients of xi(s) - Kreminski