% Session p. 148

n = 20000; p_n = 224737; b = sparse(n,1); b(1) = 1;
A = spdiags(primes(p_n)',0,n,n); e = ones(n,2);
for k = 2.^(0:floor(log2(n))), A = A + spdiags(e,[-k k],n,n); end
