% Generates Table 7.6
% number of correct digits counted by hand

startintlab;

format long e

tol = 1.5e-15;
    
n = 2000000; p = 32452843; 

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
bb = sparse(n,1); bb(1) = 1;

table_wo_sol = []; table_sol = []; kmax = 28;
for n = 20000*[1 10 100]
    A = AA(1:n,1:n); b = bb(1:n);
    D=spdiags(diag(A),0,n,n);
    [x,fail,r,k] = pcg(A,b,tol,kmax,D);
    if not(fail),
        tic,
        x1 = midrad(x(1),sup(norm(b-A*intval(x))));
        t = toc;
        % very specific code for counting the number of digits to which
        % the interval bounds agree:
        a = disp2str(inf(x1)); b = disp2str(sup(x1));
        a = a.str; a(a==' ')=[]; a(a=='.')=[]; a=a(1:find(a=='e')-1);
        b = b.str; b(b==' ')=[]; b(b=='.')=[]; b=b(1:find(b=='e')-1);
        no_corr = min(find(a~=b))-1;
        % ------------------------------------------------------------
        table_wo_sol = [table_wo_sol; n no_corr t];
        table_sol = [table_sol; infsup(x1)];
    end
end

table_wo_sol
table_sol