% Generates Table 7.2

format short g

p = 224737; n=20000;
AA = spdiags(primes(p)',0,n,n); e = ones(n,2);
for k = 2.^(0:floor(log2(n))), AA = AA + spdiags(e,[-k k],n,n); end

table = [];
for n = [20 200 2000 20000]
   A = AA(1:n,1:n);
   D2 = spdiags(1./sqrt(diag(A)),0,n,n); DA = D2*(A*D2);
   p = A(n,n);
   if n > 200   % use iterative eigenvalue solver ARPACK
       opts.issym = 1; opts.isreal = 1; opts.tol = 1e-6; opts.disp = 0;
       lmin  = eigs(A ,1,'SA',opts); lmax  = eigs(A ,1,'LA',opts);
       lmin2 = eigs(DA,1,'SA',opts); lmax2 = eigs(DA,1,'LA',opts);
   else         % use QR algorithm
       l = eig(A); l2 = eig(full(DA));
       lmin = min(l); lmin2 = min(l2);
       lmax = max(l); lmax2 = max(l2);
   end
   kappa = lmax/lmin; kappa2 = lmax2/lmin2;
   table = [table; n p lmax lmin kappa kappa2];
end

table